lu分解法、列主元高斯法、jacobi迭代法、gaussseidel法的原理及matlab程序

一、實(shí)驗(yàn)?zāi)康募邦}目1.1實(shí)驗(yàn)?zāi)康模簩W(xué)會(huì)用高斯列主元消去法，LU分解法，Jacobi迭代法和Gauss-Seidel迭代法解線性方程組。學(xué)會(huì)用Matlab編寫各種方法求解線性方程組的程序。1.2實(shí)驗(yàn)題目：用列主元消去法解方程組：x+x+3x=42x+x-x+x=1<12342x—x—x+3x=—31234—x+2x+3x—x=4用LU分解法解方程組4x=b,其中'48—240-、12"4)—242412:〔24A=，b=06202—2*—6621J)」2)3.分別用Jacobi迭代法和Gauss-Seidel迭代法求解方程組:10x—x+2x=—118x1—x\3x3=—11<2342x—x+10x=6—x+3x—x+11x=25V1234二、實(shí)驗(yàn)原理、程序框圖、程序代碼等2.1實(shí)驗(yàn)原理2.1.1高斯列主元消去法的原理Gauss消去法的基本思想是一次用前面的方程消去后面的未知數(shù)，從而將方程組化為等價(jià)形式：TOC\o"1-5"\h\z(rrrbx+bx++bx=g1111221nn1bx++bx=g<222???2nn2…以1Ibnnxn=gn這個(gè)過程就是消元，然后再回代就好了。具體過程如下：對(duì)于k=1,2,,n-1，若a(k)豐0,依次計(jì)算kkm.=a（k）/a（k）aj+i）=aj）-ma以）b（k+1）=b（k）-mb（k）,i,j-k+1,,n然后將其回代得到：‘‘‘...IIx=b（n）/a（n）x=（b（k）-況a（k）x）/a（k）,k=n一1,n一2,,1kkkjjkkIj=k+1以上是高斯消去。…但是高斯消去法在消元的過程中有可能會(huì)出現(xiàn)a（k）=0的情況，這時(shí)消元就無法進(jìn)行了，kk即使主元數(shù)a（k）n0,但是很小時(shí)，其做除數(shù)，也會(huì)導(dǎo)致其他元素?cái)?shù)量級(jí)的嚴(yán)重增長(zhǎng)和舍入誤差kk的擴(kuò)散。因此，為了減少誤差，每次消元選取系數(shù)矩陣的某列中絕對(duì)值最大的元素作為主元素。然后換行使之變到主元位置上，再進(jìn)行銷元計(jì)算。即高斯列主元消去法。2.1.2直接三角分解法（LU分解）的原理先將矩陣A直接分解為A=LU則求解方程組的問題就等價(jià)于求解兩個(gè)三角形方程組。直接利用矩陣乘法，得到矩陣的三角分解計(jì)算公式為：[u=a,i=1,2,,n[l=a/u,i=2,,nu=a一區(qū)lu,；,j=k,k+1,,nkjkjkmmj<m=1,k=2,3,nl=（a一蕓lu）/u,i=k+1,k+2,,n且knnTOC\o"1-5"\h\zikikimmkkk???m=1???由上面的式子得到矩陣A的LU分解后，求解Ux=y的計(jì)算公式為卜b1*y=b-^^ly,i=2,3,niiijjIj=1丁yn/unn…x=（y-^Lux）/u,i=n-1,n-2,,1iiijjiilj=i+1???以上為L(zhǎng)U分解法。

2.1.3Jacobi迭代法和Gauss-Seidel迭代法的原理Jcaobi迭代設(shè)線性方程組Ax=b的系數(shù)矩陣A可逆且主對(duì)角元素Oil，"??，…七n均不為零，令D=diag(a,a,...,a)1122nn并將A分解成A=（A-D）+D從而(1)可寫成Dx=（D-A）x+bx=Bx+f其中B=I-D-1&匕=D-ib.以Bi為迭代矩陣的迭代法(公式)x(k+i)=Bx(k)+f11稱為雅可比(Jacobi)迭代法，其分量形式為x（k+1）=—L[b-￡ax（k）]j=1j&iai.、ijj

j=1j=1j&11j壬i（5）i=1,2,...n,k=0,1,2,...（5）其中x(0)=q(o),x2。)’..Xo))為初始向量.（2）Gauss-Seidel迭代由雅可比迭代公式可知，在迭代的每一步計(jì)算過程中是用X&)的全部分量來計(jì)算x(k+1)的r*r*-￡--八|=rAIx、IAA-?人八x(k+1)口―一rrf>I曰-八x(k+1),...,x(k+1)',n—/—"4i+?rC.rI所有分量，顯然在計(jì)算第1個(gè)分量i時(shí)，已經(jīng)計(jì)算出的最新分量1，i-1沒有被利用。把矩陣A分解成A=D-L-U(6)其中D=dag(a,a22,…'ann),-^-分別為A的主對(duì)角元除外的下三角和上三角部分,nn于是,方程組(1)便可以寫成(D-LL=Ux+b其中以B2為迭代矩陣構(gòu)成的迭代法(公式)x(k+i)=Bx(k)+f稱為高斯一塞德爾迭代法,用分量表示的形式為x(k+偵=_L[b-另ax(k+偵-￡ax(k)]iai..ijjjj=ij=+iiii=1,2,n,k=0,1,2,...2.2程序代碼2.2.1高斯列主元的代碼functionGauss(A,b)%A為系數(shù)矩陣，b為右端項(xiàng)矩陣[m,n]=size(A);n=length(b);fork=1:n-1[pt,p]=max(abs(A(k:n,k)));%找出列中絕對(duì)值最大的數(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~=0Ab=[A,b];endendx=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);endfork=1:nfprintf('x[%d]=%f\n',k,x(k));end2.2.2LU分解法的程序functionLU(A,b)%A為系數(shù)矩陣，b為右端項(xiàng)矩陣[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);%開始進(jì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);%開始解方程組Ux=yy(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));end2.2.3Iacobi迭代法的程序functionJacobi(A,b,eps)%A為系數(shù)矩陣，b為后端項(xiàng)矩陣，epe為精度[m,n]=size(A);D=diag(diag(A));%求矩陣DL=tril(A)-D;%求矩陣LU=triu(A)-D;%求矩陣Utemp=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));end

epe為精度2.2.4Gauss-Seidel迭代程序functionGauss_Seidel(A,b,eps)%A為系數(shù)矩陣，b為后端項(xiàng)矩陣，[m,n]=size(A);D=diag(diag(A));%求矩陣DL=D-tril(A);%求矩陣LU=D-triu(A);%求矩陣Utemp=1;x=zeros(m,1);k=0;whileepe為精度temp=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í)驗(yàn)過程中需要記錄的實(shí)驗(yàn)數(shù)據(jù)表格3.1第一題(高斯列主元消去)的數(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.0000003.2第二題(匚。分解法)數(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.2596693.3第三題Jacobi迭代法的數(shù)據(jù)>>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.8423973.4第三題用Gauss_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四、實(shí)驗(yàn)中存在的問題及解決方案4.1存在的問題(1)第一題中在matlab中輸入>>Gauss(A,b)(數(shù)據(jù)省略)得到m=4n=4???Undefinedfunctionorvariable"AbH.Errorin==>Gaussat8[ap,p]=max(abs(Ab(k:n,k)));沒有得到想要的結(jié)果。(2)第二題中在matlab中輸入>>y=LU(A,b)得到y(tǒng)=4.00006.0000-5.0000-3.3571不是方程組的解。(3)第三題中在用高斯賽德爾方法時(shí)在matlab中輸入>>Gauss-Seidel(A,b,eps)結(jié)果程序報(bào)錯(cuò)??？Errorusing==>GaussToomanyoutputarguments.得不到想要的結(jié)果。4.2解決方案針對(duì)第一題中由于程序的第二行漏了一個(gè)分號(hào)導(dǎo)致輸出了m和n的值，第8行中將Ab改為A問題就解決了。由于程序后面出現(xiàn)了矩陣y故輸出的事矩陣y的值，但是我們要的事x的值，故只需要將y改成x，或者直接把y去掉就解決了問題。在function文件中命名不能出現(xiàn)“-”應(yīng)該將其改為下劃線“_”，所以將M文件名“Gauss-Seidel(A,b,eps)”改成“Gauss_Seidel(A,b,eps)”就解決問題了