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

一、實驗?zāi)康募邦}目實驗?zāi)康模簩W(xué)會用高斯列主元消去法，LU分解法，JacobiGauss-Seidel方程組。Matlab編寫各種方法求解線性方程組的程序。實驗題目：用列主元消去法解方程組：xx3x41 2 42xxxx11 2 3 42xxx3x321 2 3 4x1

2x2

3xx43 4LUAxb其中48 24 0 12 424 24 12 A

4，b 0 6 20 2

26 6 2

2 JacobiGauss-Seidel迭代法求解方程組：10xx2x

111 2 3xx3x1122 3 4xx1

10x63x1

3x2

x11x3

25二、實驗原理、程序框圖、程序代碼等實驗原理高斯列主元消去法的原理bxgbxgbxg1nn 12nn 2bxbx11

122bx222

bxgnnn n這個過程就是消元，然后再回代就好了。具體過程如下：對于k1,2,

n1，若a(k)0,依次計算kk南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實驗報告南昌航空大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院實驗報告第PAGE10第10頁ma(k)/a(k)ik ik kka(k1)a(k)m

a(k),nij ij ik,nb(k1)b(k)m

b(k),i,jk1,然后將其回代得到：

i ixb(n)/a(n)

ikk,1n n nn,1nx(b(k) a(k)x)/a(k),kn1,n2,n以上是高斯消去。

k k kj j kkjk1但是高斯消去法在消元的過程中有可能會出現(xiàn)a(k)0的情況，這時消元就無法進行了，kk即使主元數(shù)a(k)0,但是很小時，其做除數(shù)，也會導(dǎo)致其他元素數(shù)量級的嚴重增長和舍入誤差kk的擴散。因此，為了減少誤差，每次消元選取系數(shù)矩陣的某列中絕對值最大的元素作為主元素。然后換行使之變到主元位置上，再進行銷元計算。即高斯列主元消去法。直接三角分解法（LU）的原理AALU直接利用矩陣乘法，得到矩陣的三角分解計算公式為：,n2, ,nu,n2, ,nla/u

,ii1

11,nk,nukj

akj

m1

lukmmj

,,jk,k1,

n,k2,3,nl(ak1lu )/u,ik1,kik ik immk kkm1

,n且kn由上面的式子得到矩陣A的LU分解后，求解Ux=y的計算公式為yb1 1nybi1ln

y,i2,3,iiiixy

ijjj1/un n nn,1iix(ynux)/,1ii

,in1,n2,LU。

ijj iiji1JacobiGauss-SeidelJcaobiAxb (1)11A11

,a,...,a22

均不為零,令A(yù)

Ddiaga,a11 22

,...,a nn從而(1)可寫成

AADDDxDAxb

(2)令xBxf1 1其中BID1,fDb. (3)1 11以B為迭代矩陣的迭代法(公式)1xk1Bxkf1 1稱為雅可比(Jacobi)迭代法,其分量形式為

(4)1 n x(k1) i aii

bai ij1ji

x(k)j i

k(5) 1 2 其中x0x0,x0,...x 1 2 Gauss-Seidel由雅可比迭代公式可知,在迭代的每一步計算過程中是用xk的全部分量來計算xk1的ix1x1,...,x

k1沒有被利i用。把矩陣A分解成

1 i1ADLU (6)Ddiag11,a22,...,annLUA的主對角元除外的下三角和上三角部分,于是,方程組(1)便可以寫成DLxUxb即xB2其中

xf2BDL1U,22以B為迭代矩陣構(gòu)成的迭代法(公式)2

fDL1b2 (7)xk1Bxkf2 2

(8)稱為高斯—塞德爾迭代法,用分量表示的形式為1 i1 n x(k1) i a

bai

x(k1)aj

x(kjii

jji1 in, k程序代碼高斯列主元的代碼functionGauss(A,b) %A為系數(shù)矩陣，b為右端項矩陣[m,n]=size(A);n=length(b);fork=1:n-1[pt,p]=max(abs(A(k:n,k))); %找出列中絕對值最大的數(shù)p=p+k-1;ifp>kt=A(k,:);A(k,:)=A(p,:);A(p,:)=t; %t=b(k);b(k)=b(p);b(p)=t;endm=A(k+1:n,k)/A(k,k); %開始消元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);ifflag~=0endend

Ab=[A,b];x=zeros(n,1); %x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);endforfprintf('x[%d]=%f\n',k,x(k));endLU分解法的程序functionLU(A,b) %A為系數(shù)矩陣，b為右端項矩[m,n]=size(A); %初始化矩陣A，b，L和Un=length(b);L=eye(n,n);U=zeros(n,n);U(1,1:n)=A(1,1:n); %LU分解L(2:n,1)=A(2:n,1)/U(1,1);fork=2:nU(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);endL %L矩陣U %U矩陣y=zeros(n,1); %y(1)=b(1);fork=2:ny(k)=b(k)-L(k,1:k-1)*y(1:k-1);endx=zeros(n,1);x(n)=y(n)/U(n,n);fork=n-1:-1:1x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);endfork=1:nfprintf('x[%d]=%f\n',k,x(k));endJacobifunctionJacobi(A,b,eps) %A為系數(shù)矩陣，b為后端項矩陣，epe為精度[m,n]=size(A);D=diag(diag(A)); %求矩陣L=D-tril(A); %求矩陣LU=D-triu(A); %temp=1;x=zeros(m,1);k=0;whileabs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %記錄循環(huán)次數(shù)x=inv(D)*(L+U)*x+inv(D)*b; %雅克比迭代公endfork=1:nfprintf('x[%d]=%f\n',k,x(k));endGauss-SeidelfunctionGauss_Seidel(A,b,eps) %A為系數(shù)矩陣，b為后端項矩陣，epe為精度[m,n]=size(A);D=diag(diag(A)); %求矩陣L=D-tril(A); %求矩陣LU=D-triu(A); %temp=1;x=zeros(m,1);k=0;whileabs(max(x)-temp)>epstemp=max(abs(x));k=k+1; %記錄循環(huán)次數(shù)x=inv(D-L)*U*x+inv(D-L)*b; %Gauss_Seidel的迭代公式endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end三、實驗過程中需要記錄的實驗數(shù)據(jù)表格第一題（高斯列主元消去）的數(shù)據(jù)>>A=[1103;21-11;3-1-13;-123-1];>>b=[4;1;-3;4];>>Gauss(A,b)x[1]=-1.333333x[2]=2.333333x[3]=-0.333333x[4]=1.000000第二題（LU分解法）數(shù)據(jù)>>A=[48-240-12;-24241212;06202;-66216];>>b=[4;4;-2;-2];>>LU(A,b)L=1.0000000-0.50001.00000000.50001.00000-0.12500.2500-0.07141.0000U=48.0000-24.00000-12.0000012.000012.00006.00000014.0000-1.000000012.9286x[1]=0.521179x[2]=1.005525x[3]=-0.375691x[4]=-0.259669Jacobi>>A=[10-120;08-13;2-1100;-13-111];b=[-11;-11;6;25];Jacobi(A,b,0.00005)x[1]=-1.467396x[2]=-2.358678x[3]=0.657604x[4]=2.842397Gauss_Seidel迭代的數(shù)據(jù)>>A=[10-120;08-13;2-1100;-13-111];>>b=[-11;-11;6;25];>>Gauss_Seidel(A,b,0.00005)x[1]=-1.467357x[2]=-2.358740x[3]=0.657597x[4]=2.842405四、實驗中存在的問題及解決方案存在的問題第一題中在matlab中輸入>>數(shù)據(jù)省略得到m=4 n=4 Undefinedfunctionorvariable"Ab".Errorin==>Gaussat8[ap,p]=max(abs(Ab(k:n,k)));沒有得到想要的果。（2）第二題中在 matlab中輸入>>y=LU(A,b)得到y(tǒng)=4.0000 6.0000 -5.0000-3.3571不是方程組的解。（3）第三題中在用高斯賽德爾方法時在matlab中輸入>>Gauss-Seidel(A,b,eps)結(jié)果程序錯Errorusing==>Gauss Toomanyoutput得不到想要的結(jié)果。解決方案mn8行中將Ab改為A問題就解決了。由于程序后面出現(xiàn)了矩陣yy的值，但是我們要的事xyxy去掉就解決了問題。function文件中命名不能出現(xiàn)“”應(yīng)該將其改為下劃線“_M“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel