第一篇：LU分解法、列主元高斯法、Jacobi迭代法、Gauss-Seidel法的原理及Matlab程序

一、實(shí)驗(yàn)?zāi)康募邦}目

1.1 實(shí)驗(yàn)?zāi)康模?/p>

（1）學(xué)會(huì)用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seidel迭代法解線(xiàn)性方程組。

（2）學(xué)會(huì)用Matlab編寫(xiě)各種方法求解線(xiàn)性方程組的程序。1.2 實(shí)驗(yàn)題目：

1.用列主元消去法解方程組：

?x1?x2?3x4?4?2x?x?x?x?1?123

4?

2x?x?x?3x??34?123???x1?2x2?3x3?x4?42.用LU分解法解方程組Ax?b,其中

12?48?240???4??????24241212?，b??4?

A???0???2?6202?????66216????2?3.分別用Jacobi迭代法和Gauss-Seidel迭代法求解方程組：

?10x1?x2?2x3??11?8x?x?3x??11?234

?

2x?x?10x?63?12???x1?3x2?x3?11x4?25

二、實(shí)驗(yàn)原理、程序框圖、程序代碼等

2.1實(shí)驗(yàn)原理

2.1.1高斯列主元消去法的原理

Gauss消去法的基本思想是一次用前面的方程消去后面的未知數(shù)，從而將方程組化為等價(jià)形式：

?b11x1?b12x2???b1nxn?g1?b22x2???b2nxn?g2? ????bnnxn?gn?這個(gè)過(guò)程就是消元，然后再回代就好了。具體過(guò)程如下：

(k)對(duì)于k?1,2,?,n?1，若akk?0,依次計(jì)算

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

(k)(k)mik?aik/akk(k?1)(k)(k)aij?aij?mikakj

bi(k?1)?bi(k)?mikbk(k),i,j?k?1,?,n然后將其回代得到：

(n)(n)?xn?bn/ann?n ?(k)(k)(k)x?(b?ax)/a,k?n?1,n?2,?,1?kkjjkk?kj?k?1?以上是高斯消去。

(k)但是高斯消去法在消元的過(guò)程中有可能會(huì)出現(xiàn)akk?0的情況，這時(shí)消元就無(wú)法進(jìn)行了，(k)即使主元數(shù)akk?0,但是很小時(shí)，其做除數(shù)，也會(huì)導(dǎo)致其他元素?cái)?shù)量級(jí)的嚴(yán)重增長(zhǎng)和舍入誤差的擴(kuò)散。因此，為了減少誤差，每次消元選取系數(shù)矩陣的某列中絕對(duì)值最大的元素作為主元素。然后換行使之變到主元位置上，再進(jìn)行銷(xiāo)元計(jì)算。即高斯列主元消去法。2.1.2直接三角分解法（LU分解）的原理

先將矩陣A直接分解為A?LU則求解方程組的問(wèn)題就等價(jià)于求解兩個(gè)三角形方程組。直接利用矩陣乘法，得到矩陣的三角分解計(jì)算公式為：

?u1i?a1i,i?1,2,?,n??li1?ai1/u11,i?2,?,nk?1? ?ukj?akj??lkmumj，j?k,k?1,?,n?m?1,k?2,3,?n?k?1?l?(a?lu)/u,i?k?1,k?2,?,n且k?n?ikikimmkkk?m?1?由上面的式子得到矩陣A的LU分解后，求解Ux=y的計(jì)算公式為

?y1?b1?i?1??yi?bi??lijyj,i?2,3,?nj?1?

?xn?yn/unn?n??xi?(yi??uijxj)/uii,i?n?1,n?2,?,1j?i?1?以上為L(zhǎng)U分解法。

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理（1）Jcaobi迭代

設(shè)線(xiàn)性方程組

Ax?b(1)的系數(shù)矩陣A可逆且主對(duì)角元素a11,a22,...,ann均不為零,令

D?diag?a11,a22,...,ann并將A分解成

A??A?D??D(2)從而(1)可寫(xiě)成

Dx??D?A?x?b 令

x?B1x?f1

?1?1B?I?DA,f?Db.(3)1其中1?

以B1為迭代矩陣的迭代法(公式)

x?k?1??B1x?k??f1(4)稱(chēng)為雅可比(Jacobi)迭代法,其分量形式為

???x(k?1)i1?aii?bi??aijx(jk)j?1j?in?i?1,2,...n,Tk?0,1,2,...(5)

?0??0??0??0???為初始向量.x?x,x,...x12n其中（2）Gauss-Seidel迭代

由雅可比迭代公式可知,在迭代的每一步計(jì)算過(guò)程中是用x所有分量,顯然在計(jì)算第i個(gè)分量xi用。

把矩陣A分解成

?k?1??k?的全部分量來(lái)計(jì)算x?k?1??k?1?的時(shí),已經(jīng)計(jì)算出的最新分量x1,...,xi?1?k?1?沒(méi)有被利

A?D?L?U(6)其中D?diag?a11,a22,...,ann?,?L,?U分別為A的主對(duì)角元除外的下三角和上三角部分,第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

于是,方程組(1)便可以寫(xiě)成

?D?L?x?Ux?b 即

x?B2x?f2

其中

B2??D?L?U,?1f2??D?L?b(7)

?1以B2為迭代矩陣構(gòu)成的迭代法(公式)

x?k?1??B2x?k??f2(8)稱(chēng)為高斯—塞德?tīng)柕?用分量表示的形式為

??? x(k?1)ii?1n1(k?1)??bi??aijxj??aijx(jk)j?1j?i?1aii? i?1,2,?n,k?0,1,2,...2.2程序代碼

2.2.1高斯列主元的代碼

function Gauss(A,b)

%A為系數(shù)矩陣，b為右端項(xiàng)矩陣 [m,n]=size(A);n=length(b);for k=1:n-1

[pt,p]=max(abs(A(k:n,k)));

%找出列中絕對(duì)值最大的數(shù)

p=p+k-1;

if p>k

t=A(k,:);A(k,:)=A(p,:);A(p,:)=t;

%交換行使之變到主元位置上

t=b(k);b(k)=b(p);b(p)=t;

end

m=A(k+1:n,k)/A(k,k);

%開(kāi)始消元

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);if flag~=0

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

Ab=[A,b];end end

x=zeros(n,1);

%開(kāi)始回代

x(n)=b(n)/A(n,n);

for k=n-1:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

for k=1:n

fprintf('x[%d]=%fn',k,x(k));

end 2.2.2 LU分解法的程序

function LU(A,b)

%A為系數(shù)矩陣，b為右端項(xiàng)矩陣 [m,n]=size(A);

%初始化矩陣A，b，L和U n=length(b);

L=eye(n,n);U=zeros(n,n);U(1,1:n)=A(1,1:n);

%開(kāi)始進(jìn)行LU分解 L(2:n,1)=A(2:n,1)/U(1,1);for k=2:n

U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);

L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);

end L

%輸出L矩陣 U

%輸出U矩陣

y=zeros(n,1);

%開(kāi)始解方程組Ux=y y(1)=b(1);for k=2:n

y(k)=b(k)-L(k,1:k-1)*y(1:k-1);end x=zeros(n,1);

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

x(n)=y(n)/U(n,n);for k=n-1:-1:1

x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);end for k=1:n

fprintf('x[%d]=%fn',k,x(k));end 2.2.3 Jacobi迭代法的程序

function Jacobi(A,b,eps)

%A為系數(shù)矩陣，b為后端項(xiàng)矩陣，epe為精度 [m,n]=size(A);D=diag(diag(A));

%求矩陣D L=D-tril(A);

%求矩陣L U=D-triu(A);

%求矩陣U temp=1;x=zeros(m,1);k=0;while abs(max(x)-temp)>eps

temp=max(abs(x));

k=k+1;

%記錄循環(huán)次數(shù)

x=inv(D)*(L+U)*x+inv(D)*b;

%雅克比迭代公式 end

for k=1:n

fprintf('x[%d]=%fn',k,x(k));end 2.2.4 Gauss-Seidel迭代程序

function Gauss_Seidel(A,b,eps)

%A為系數(shù)矩陣，b為后端項(xiàng)矩陣，epe為精度 [m,n]=size(A);D=diag(diag(A));

%求矩陣D L=D-tril(A);

%求矩陣L U=D-triu(A);

%求矩陣U temp=1;

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

x=zeros(m,1);k=0;while abs(max(x)-temp)>eps

temp=max(abs(x));

k=k+1;

%記錄循環(huán)次數(shù)

x=inv(D-L)*U*x+inv(D-L)*b;

%Gauss_Seidel的迭代公式 end

for k=1:n

fprintf('x[%d]=%fn',k,x(k));end

三、實(shí)驗(yàn)過(guò)程中需要記錄的實(shí)驗(yàn)數(shù)據(jù)表格

3.1第一題（高斯列主元消去）的數(shù)據(jù)

>> A=[1 1 0 3;2 1-1 1;3-1-1 3;-1 2 3-1];>> b=[4;1;-3;4];>> Gauss(A,b)x[1]=-1.333333 x[2]=2.333333 x[3]=-0.333333 x[4]=1.000000

3.2 第二題（LU分解法）數(shù)據(jù)

>> A=[48-24 0-12;-24 24 12 12;0 6 20 2;-6 6 2 16];>> b=[4;4;-2;-2];>> LU(A,b)L =

1.0000

0

0

0

-0.5000

1.0000

0

0

0

0.5000

1.0000

0

-0.1250

0.2500

-0.0714

1.0000 U =

48.0000-24.0000

0-12.0000

0

12.0000

12.0000

6.0000

0

0

14.0000

-1.0000

0

0

0

12.9286

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

x[1]=0.521179 x[2]=1.005525 x[3]=-0.375691 x[4]=-0.259669

3.3第三題Jacobi迭代法的數(shù)據(jù)

>> A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];b=[-11;-11;6;25];Jacobi(A,b,0.00005)x[1]=-1.467396 x[2]=-2.358678 x[3]=0.657604 x[4]=2.842397

3.4第三題用Gauss_Seidel迭代的數(shù)據(jù)

>> A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11];>> b=[-11;-11;6;25];>> Gauss_Seidel(A,b,0.00005)x[1]=-1.467357 x[2]=-2.358740 x[3]=0.657597 x[4]=2.842405

四、實(shí)驗(yàn)中存在的問(wèn)題及解決方案

4.1存在的問(wèn)題

（1）第一題中在matlab中輸入>> Gauss(A,b)（數(shù)據(jù)省略）得到m =4 n =4 ??? Undefined function or variable “Ab”.Error in ==> Gauss at 8[ap,p]=max(abs(Ab(k:n,k)));沒(méi)有得到想要的結(jié)果。

（2）第二題中在matlab中輸入>> y=LU(A,b)得到y(tǒng) =4.0000

6.0000

-5.0000

-3.3571不是方程組的解。

（3）第三題中在用高斯賽德?tīng)柗椒〞r(shí)在matlab中輸入>> Gauss-Seidel(A,b,eps)結(jié)果程序報(bào)錯(cuò)??? Error using ==> Gauss

Too many output arguments.得不到想要的結(jié)果。

第 頁(yè)

南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實(shí)驗(yàn)報(bào)告

4.2解決方案

（1）針對(duì)第一題中由于程序的第二行漏了一個(gè)分號(hào)導(dǎo)致輸出了m和n的值，第8行中將Ab改為A問(wèn)題就解決了。

（2）由于程序后面出現(xiàn)了矩陣y故輸出的事矩陣y的值，但是我們要的事x的值，故只需要將y改成x，或者直接把y去掉就解決了問(wèn)題。

（3）在function文件中命名不能出現(xiàn)“-”應(yīng)該將其改為下劃線(xiàn)“_”，所以將M文件名“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解決問(wèn)題了。

五、心得體會(huì)

本次試驗(yàn)涉及到了用高斯列主元消去法，LU分解法，Jacobi迭代法以及Gauss-seidel迭代法等四種方法。需要對(duì)這些方法的原理都要掌握才能寫(xiě)出程序，由于理論知識(shí)的欠缺，我花了很大一部分時(shí)間在看懂實(shí)驗(yàn)的原理上，看懂了實(shí)驗(yàn)原理之后就開(kāi)始根據(jù)原理編寫(xiě)程序，程序中還是出現(xiàn)了很多的低級(jí)錯(cuò)誤導(dǎo)致調(diào)試很久才能運(yùn)行。

通過(guò)這次試驗(yàn)使我深刻的體會(huì)到理論知識(shí)的重要性，沒(méi)有理論知識(shí)的支撐是寫(xiě)不出程序來(lái)的。寫(xiě)程序時(shí)還會(huì)犯很多低級(jí)的錯(cuò)誤，以后一定要加強(qiáng)理論知識(shí)的學(xué)習(xí)，減少編程時(shí)低級(jí)錯(cuò)誤的產(chǎn)生。

第 頁(yè)