列主元高斯消去法和列主元三角分解法解線性方程

1、計(jì)算方法實(shí)驗(yàn)報(bào)告1【課題名稱】用列主元高斯消去法和列主元三角分解法解線性方程【目的和意義】高斯消去法是一個(gè)古老的求解線性方程組的方法,但由它改進(jìn)得到的選主元的高斯消去法則是目前計(jì)算機(jī)上常用的解低階稠密矩陣方程組的有效方法。用高斯消去法解線性方程組的基本思想時(shí)用矩陣行的初等變換將系數(shù)矩陣A約化為具有簡(jiǎn)單形式的矩陣（上三角矩陣、單位矩陣等），而三角形方程組則可以直接回帶求解用高斯消去法解線性方程組 Ax b （其中A Rn xn）的計(jì)算量為：乘除法運(yùn)算步驟為MD n (n 1) (n 1) n (2 n 1) n (n 1) n(n 1) n 3* 22 6 2 2n3，加減運(yùn)算步驟為AS (n

2、1) n (2 n 1 ) n( n 1 ) n (n 1)(n 1) n (2 n 2 25）_6。相比之下，傳統(tǒng)的克萊姆U “ c19法則則較為繁瑣，如求解 20階線性方程組，克萊姆法則大約要5 10次乘法，而用高斯消去法只需要3060次乘除法。在高斯消去法運(yùn)算的過(guò)程中，如果出現(xiàn)abs（A（i,i）等于零或過(guò)小的情況，則會(huì)導(dǎo)致矩陣元素?cái)?shù)量級(jí)嚴(yán)重增長(zhǎng)和舍入誤差的擴(kuò)散，使得最后的計(jì)算結(jié)果不可靠，所以目前計(jì)算機(jī)上常用的解低階稠密矩陣方程的快速有效的方法時(shí)列主元高斯消去法,從而使計(jì)算結(jié)果更加精確。2、列主元三角分解法高斯消去法的消去過(guò)程，實(shí)質(zhì)上是將A分解為兩個(gè)三角矩陣的乘積A=LU，并求解Ly=b

3、的過(guò)程?；貛н^(guò)程就是求解上三角方程組Ux=y。所以在實(shí)際的運(yùn)算中，矩陣L和U可以直接計(jì)算出，而不需要任何中間步驟， 從而在計(jì)算過(guò)程中將高斯消去法的步驟進(jìn)行了進(jìn)一步的簡(jiǎn)略，大大提咼了運(yùn)算速度，這就是三角分解法采用選主元的方式與列主元高斯消去法一樣,也是為了避免除數(shù)過(guò)小，從而保證了計(jì)算的精確度【計(jì)算公式】1、列主元高斯消去法設(shè)有線性方程組Ax=b，其中設(shè)A為非奇異矩陣。方程組的增廣矩陣為aii3,2a1na,ba21Ma22a2nbl b2 M第1步（k=1）:首先在然后交換（A，b ）的第3|1Man1an2annbnA的第一列中選取絕對(duì)值最大的元素3|1，作為第一步的主元素：311max3&#

4、39;ii1行與第I行元素，再進(jìn)行消元計(jì)算。設(shè)列主元素消去法已經(jīng)完成第1步到第k-1步的按列選主元，交換兩行，消元計(jì)算得到與原方程組等價(jià)的方程組A（k）x=b（k）A,b A(k),b(k)L瞪31(n)Mb32?L3；?322)Mb22)OMMM3kk)L3kn)Mbkk)(k)3k 1,k(k)3k 1,nMbkk)MMMM31（1）第k步計(jì)算如下:對(duì)于k=1 , 2，n-1(k)(k)atkmaxaik(1 )按列選主元：即確定 t(2)如果t Mk，則交換A ,b第t行與第k行元素。消元計(jì)算aikmikaik,(iakk1,L , n)abimk akj,mkbk,(i,j k 1L

5、,n)(ik 1L ,n)£ l,(z = A -h 1,消元乘數(shù)mik滿足:(4 )回代求解XnbnannXi(bij iaijXj)1,(i n 1, n 2,1)aii步分解，為了避免用絕對(duì)值很小的數(shù)Ukk作除數(shù),Ukjakjlik(aikk 1lkmumjm 1k 1limumk) / ukkm 1(j(iaik1limUmk(ik,km 11,Lk, k 1,L , n; k 1,2丄,n)k 1,k2,L ,n;k 1,2丄，n),n），于是有Ukk= Sk。如果Stmaxs，則將矩陣的第t行與第k行元素互換，將（i. a -i，j）位置的新元素仍記為 jj或jj，然后再

6、做第k步分解，這2、列主元三角分解法AA,b對(duì)方程組的增廣矩陣經(jīng)過(guò)k-1步分解后，可變成如下形式:u11U12LU1,k 1u1kLu1jL%y1I21U22LU2,k 1U2kLU2jLU2ny2MMMMMMMlk 1,1lk 1,2LUk 1,k 1Uk 1,kLUk 1,jLUk 1,nyk1Alk1lk2Llk,k 1akkLakjLaknbkMMMMMMMli1li2Lli ,k 1aikLajLanbMMMMMMMln1ln2Lln,k 1ankLanjLannbnukksklikSi /Sk且lik1 (i(Sk即交換前的St), (i k 1,k 2,L ,n)k 1,k2,L

7、 , n),【列主元高斯消去法程序流程圖】輸入方程：系數(shù)矩陪A，t以辰主元允許S小值c程組是舌肖k-1,2n-1按別選刖主云亙T標(biāo)p, q=in 3Dt 3 s（；A（tt:切J并將tTAClcnh:）巾的行號(hào)4換苗在丸中的行、即（=q+k-li=k+-lJ . . -nm0kj=AG*j/A<k/E3A«k :環(huán)-凸(yk -in G10 * 4畑 k: *回代求紹i=n- J. 1%Tij=hCrH咼 Qxn)城：=0(手-BLHtiJ；攻:ii+l ：可*%i+l 榊沖 ©p J 諭出i+算幫陪機(jī)IKCJajQZJ (p) VCVES4輸岀“方程處有誤"

8、;停機(jī)VE3判斷列午訐旱否存當(dāng)前行k輔岀“主元太小“，弘L < A> =0t?Knpltmp2=bCSr；工bCq,)=temps.【列主元高斯消去法Matlab主程序】function x=gauss1(A,b,c)%列主元法高斯消去法解線性方程Ax=bif (le ngth(A)=le ngth(b)%判斷輸入的方程組是否有誤enddisp('輸入方程有誤!return;disp('原方程為 AX=b :')dis PCn=len gth(A);for k=1: n-1')p ,q=max(abs(A (k:n, k);q=q+k-1;if ab

9、s (p)<cdisp('列元素太小，det(A) y);break;elseif q>ktem P仁A(k,:);%顯示方程組%找列主元%找出第k列中的最大值，其下標(biāo)為p,q%q在A(k:n,k)中的行號(hào)轉(zhuǎn)換為在A中的行%列主元所在行不是當(dāng)前行，將當(dāng)前行與列A(k,:)=A(q,:);元所在行交換(包括b)endA(q,:)=tem p1;tem p2=b(k,:);b(k,:)=b(q,:);b(q,:)=te mp2;%消元for i=k+1: nm(i,k)=A(i,k)/A(k,k);%A(k,k) 將 A(i,k)消為0所乘系數(shù)%第i行消元處理%b消元處理A(i

10、,k: n)=A(i,k: n)-m(i,k)*A(k,k: n);b(i)=b(i)-m(i,k)*b(k);endenddisp('消元后所得到的上三角陣是')%顯示消元后的系數(shù)矩陣%回代求解b(n )=b( n)/A(n,n);for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1: n)*b(i+1: n)/A(i,i);end clear x;disp('AX=b的解 x 是')x=b;【調(diào)用函數(shù)解題】>>>>>>>>Command WindowA= 3 4;1 -1 1.212： b=

11、i ;2；3;c=0,0000 L .?r= £au331 (Aj b, c)原方程為AS二b：-11消元后師得到的上三角眩是-0.3333Q, 5000【列主元三角分解法程序流程圖】【列主元三角分解法Matlab主程序】自己編的程序:%定義函數(shù)列主元三角分解法函fun ction x=P LU(A,b,e ps)%判斷輸入的方程組是否有誤if (le ngth(A)=le ngth(b)dispC輸入方程有誤！')return;enddisp('原方程為 AX=b :')%顯示方程組dis PCn=le ngth(A);A=A b;%將A與b合并，得到增廣矩

12、陣for r=1: nif r=1for i=1:nc d=max(abs(A(:,1);%選取最大列向量，并做行交換if c<=e ps%最大值小于e，主元太小，程序結(jié)束break;elseendd=d+1-1;p=A(1,:);A(1,:)=A(d,:);A(d,:)=p;A(1,i)=A(1,i);endA(1,2: n)=A(1,2: n);% 求 u(1,i)A(2: n,1)=A(2: n,1)/A(1,1);else%按照方程求取u(r,i)u(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r);if abs(u(r,r)v=e ps%如果u (r,r)小于

13、e,則交換行P=A(r,:);A(r,:)=A(r+1,:);A(r+1,:)=p;elseendfor i=r: nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i);%根據(jù)公式求解，并把結(jié)果存在矩陣Aendfor i=r+1: nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r)/A(r,r);%根據(jù)公式求解，并把結(jié)果存在矩陣endendendy(1)=A(1, n+1);for i=2:nh=0;for k=1:i-1h=h+A(i,k)*y(k);endendy(i)=A(i, n+1)-h;%根據(jù)公式求解y(i)x(n )=y( n)/A(n

14、,n);for i=n-1:-1:1h=0;for k=i+1: nh=h+A(i,k)*x(k);endx(i)=(y(i)-h)/A(i,i);%根據(jù)公式求解x(i)end disp('AX=b 的解 x 是')x=x'%輸出方程的解可直接得到P,L,U并解出方程解的的程序（查閱資料得子函數(shù) PLU1，其作用是將矩陣A分解成L乘以U的形式。PLU2為調(diào)用PLU1解題的程序，是自己編的）(I ).fun ctio n l,u, p=PLU1(A)%定義子函數(shù)，其功能為列主元三角分解系數(shù)矩m, n=size(A);%判斷系數(shù)矩陣是否為方if m=nerror（'

15、矩陣不是方陣'）return%判斷系數(shù)矩陣endif det(A)=0能否被三角分解error（'矩陣不能被三角分解'）endPL,Uu=A;p=eye(m);匸eye(m);%將系數(shù)矩陣三角分解，分別求出for i=1:mfor j=i:mt(j)=u(j,i);for k=1:i-1endendt(j)=t(j)-u(j,k)*u(k,i);a=i;b=abs(t(i);for j=i+1:mif b<abs(t(j)b=abs(t(j);a=j；endendif a=ifor j=1:mc=u(i,j);u(i,j)=u(a,j);u(a,j)=c;endf

16、or j=1:mP(i，j)=P(a,j)；P (a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endu(i,i)=t(i);for j=i+1:mu(j,i)=t(j)/t(i);endfor j=i+1:mfor k=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);endendend l=tnl(u,-1)+eye(m);u=triu(u,0)%定義列主元三角分解法的fun ction x=P LU2(A,b)函數(shù)%調(diào)用PLU分解系數(shù)矩l,u, p=PLU1(A)m=le ngth(A);%由于A左乘p,故b也要左v=b;for q=1:mb(q)=

17、sum( p(q,1:m)*v(1:m,1);end%求解方程b(1)=b(1)Ly=b for i=2:1:mb(i)=(b(i)-sum(l(i,1:i-1)*b(1:i-1);end%求解方程b(m)=b(m)/u(m,m);Ux=y for i=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)/u(i,i);end clear x;disp('AX=b 的解 x 是')x=b;【調(diào)用函數(shù)解題】Command Window» a= 03-1» b=l：2：31» FLU(ab, th 0001)I.leeT-0.3333Q. EOOOCommand Window1 ;2 1>>£=<! 3 4； 1 -» b=El：2；3:>> s=PLU2CA, b)1,00001-00000 5001