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

,.計(jì)算方法實(shí)驗(yàn)報(bào)告1【課題名稱】用列主元高斯消去法和列主元三角分解法解線性方程【目的和意義】高斯消去法是一個(gè)古老的求解線性方程組的方法，但由它改進(jìn)得到的選主元的高斯消去精品文檔放心下載法則是目前計(jì)算機(jī)上常用的解低階稠密矩陣方程組的有效方法。精品文檔放心下載用高斯消去法解線性方程組的基本思想時(shí)用矩陣行的初等變換將系數(shù)矩陣A約化為具精品文檔放心下載有簡單形式的矩陣（上三角矩陣、單位矩陣等），而三角形方程組則可以直接回帶求解謝謝閱讀用高斯消去法解線性方程組Axb（其中A∈Rn×n）的計(jì)算量為：乘除法運(yùn)算步驟為精品文檔放心下載n(n1)(n1)n(2n1)n(n1)n(n1)n3nMD26223n23，加減運(yùn)算步驟為AS(n1)n(2n1)n(n1)n(n1)(n1)n(2n5)6226。相比之下，傳統(tǒng)的克萊姆法則則較為繁瑣，如求解20階線性方程組，克萊姆法則大約要51019次乘法，而用高斯消去法只需要3060次乘除法。謝謝閱讀在高斯消去法運(yùn)算的過程中，如果出現(xiàn)abs(A(i,i))等于零或過小的情況，則會(huì)導(dǎo)致矩陣謝謝閱讀元素?cái)?shù)量級嚴(yán)重增長和舍入誤差的擴(kuò)散，使得最后的計(jì)算結(jié)果不可靠，所以目前計(jì)算機(jī)上常精品文檔放心下載用的解低階稠密矩陣方程的快速有效的方法時(shí)列主元高斯消去法，從而使計(jì)算結(jié)果更加精確。精品文檔放心下載2、列主元三角分解法高斯消去法的消去過程，實(shí)質(zhì)上是將A分解為兩個(gè)三角矩陣的乘積A=LU，并求解Ly=b謝謝閱讀的過程?；貛н^程就是求解上三角方程組Ux=y。所以在實(shí)際的運(yùn)算中，矩陣L和U可以直精品文檔放心下載接計(jì)算出，而不需要任何中間步驟，從而在計(jì)算過程中將高斯消去法的步驟進(jìn)行了進(jìn)一步的謝謝閱讀,.簡略，大大提高了運(yùn)算速度，這就是三角分解法采用選主元的方式與列主元高斯消去法一樣，也是為了避免除數(shù)過小，從而保證了計(jì)算的精謝謝閱讀確度【計(jì)算公式】1、列主元高斯消去法設(shè)有線性方程組Ax=b，其中設(shè)A為非奇異矩陣。方程組的增廣矩陣為感謝閱讀aaLLaMb11121n1aaLLaMb21222n2MM[a,b]al1MMMMMMaaLLaMbn1n2nnn第1步（k=1）：首先在A的第一列中選取絕對值最大的元素al1，作為第一步的主元素:精品文檔放心下載amaxa0l11ini1然后交換（A，b）的第1行與第l行元素，再進(jìn)行消元計(jì)算。精品文檔放心下載設(shè)列主元素消去法已經(jīng)完成第1步到第k-1步的按列選主元，交換兩行，消元計(jì)算得到與感謝閱讀原方程組等價(jià)的方程組 A(k)x=b(k)a(1) 11[A,b][A(k),b(k)]精品文檔放心下載

a(1)La(1)a(1)Mb(1)121k1n1a(2)La(2)a(2)Mb(2)222k2n2OMMMa(k)La(k)Mb(k)kkknka(k)a(k)Mb(k)k1,kk1,nk1MMMM,.第k步計(jì)算如下：對于k=1，2，…，n-1a(k)maxa(k)0tkkinik（1）按列選主元：即確定t使（2）如果t≠k，則交換[A，b]第t行與第k行元素。感謝閱讀（3）消元計(jì)算ama,(ik1,L,n)ikikikakkaama, ij ij ik kj bbmb,i i ik k

(i,jk1,L,n)(ik1,L,n)消元乘數(shù)mik滿足:（4）回代求解xbnnann(bnax)ijjixji1,(in1,n2,,1)iaii,.2、列主元三角分解法A[A,b]對方程組的增廣矩陣經(jīng)過k-1步分解后，可變成如下形式：uuLuuLuLuy11121,k11k1j1n1l21u22Lu2,k1u2kLu2jLu2ny2MMMMMMMlk1,1lLuuLuLuyk1,2k1,k1k1,kk1,jk1,nk1AllLlaLaLabk1k2k,k1kkkjknkMMMMMMMllLlaLaLabi1i2i,k1ikijiniMMMMMMMllLlaLaLabn1n2n,k1nknjnnn第k步分解，為了避免用絕對值很小的數(shù)ukk作除數(shù)，引進(jìn)量謝謝閱讀ualu(jk,k1,L,n;k1,2,L,n)k1kjkjkmmjm1l(ak1u)/u(ik1,k2,L,n;k1,2,L,n)likikimmkkkm1k1lu(ik,k1,L,n)saiikimmk，于是有ukk=sk。如果smaxs，則將矩陣的第m1tkini行與第k行元素互換，將（i，j）位置的新元素仍記為ljj或ajj，然后再做第k步分解，這感謝閱讀時(shí),.u s (s即交換前的s)，kk k k ts/s（ik1,k2,L,n)感謝閱讀ik i kl1(ik1,k2,L,n),感謝閱讀ik【列主元高斯消去法程序流程圖】【列主元高斯消去法Matlab主程序】functionx=gauss1(A,b,c) %列主元法高斯消去法解線性方程Ax=b感謝閱讀,.if(length(A)~=length(b)) %判斷輸入的方程組是否有誤謝謝閱讀disp('輸入方程有誤！')return;enddisp('原方程為AX=b：') %顯示方程組Abdisp('------------------------')精品文檔放心下載n=length(A);fork=1:n-1 %找列主元[p,q]=max(abs(A(k:n,k))); %找出第k列中的最大值，其下標(biāo)為[p,q]精品文檔放心下載q=q+k-1; %q在A(k:n,k)中的行號轉(zhuǎn)換為在A中的行感謝閱讀號ifabs(p)<cdisp('列元素太小，det(A)≈0');break;elseifq>ktemp1=A(k,:); %列主元所在行不是當(dāng)前行，將當(dāng)前行與列謝謝閱讀主,.A(k,:)=A(q,:); 元所在行交換(包括b)精品文檔放心下載A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=temp2;end%消元fori=k+1:nm(i,k)=A(i,k)/A(k,k); %A(k,k)將A(i,k)消為0所乘系數(shù)感謝閱讀A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); %第i行消元處理感謝閱讀b(i)=b(i)-m(i,k)*b(k); %b消元處理謝謝閱讀endenddisp('消元后所得到的上三角陣是')A %顯示消元后的系數(shù)矩陣b(n)=b(n)/A(n,n); %回代求解fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n)*b(i+1:n)))/A(i,i);精品文檔放心下載endclearx;,.disp('AX=b的解x是')x=b;【調(diào)用函數(shù)解題】,.【列主元三角分解法程序流程圖】,.【列主元三角分解法Matlab主程序】①自己編的程序：functionx=PLU(A,b,eps) %定義函數(shù)列主元三角分解法函謝謝閱讀數(shù)if(length(A)~=length(b)) %判斷輸入的方程組是否有誤精品文檔放心下載,.disp('輸入方程有誤！')return;enddisp('原方程為AX=b：') %顯示方程組Abdisp('------------------------')感謝閱讀n=length(A);A=[Ab]; %將A與b合并，得到增廣矩陣謝謝閱讀forr=1:nifr==1fori=1:n[cd]=max(abs(A(:,1))); %選取最大列向量，并做行交換精品文檔放心下載ifc<=eps %最大值小于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);A(2:n,1)=A(2:n,1)/A(1,1); %求u(1,i)精品文檔放心下載elseu(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r); %按照方程求取u(r,i)謝謝閱讀ifabs(u(r,r))<=eps %如果u（r,r）小于e，則交換行感謝閱讀p=A(r,:);A(r,:)=A(r+1,:);A(r+1,:)=p;elseendfori=r:nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根據(jù)公式求解，并把結(jié)果存在矩陣A謝謝閱讀中endfori=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r))/A(r,r);%根據(jù)公式求解，并把結(jié)果存在矩陣A中精品文檔放心下載endendend,.y(1)=A(1,n+1);fori=2:nh=0;fork=1:i-1h=h+A(i,k)*y(k);endy(i)=A(i,n+1)-h; %根據(jù)公式求解y(i)謝謝閱讀endx(n)=y(n)/A(n,n);fori=n-1:-1:1h=0;fork=i+1:nh=h+A(i,k)*x(k);endx(i)=(y(i)-h)/A(i,i); %根據(jù)公式求解x(i)感謝閱讀endAdisp('AX=b的解x是')x=x'; %輸出方程的解,.②可直接得到P,L,U并解出方程解的的程序（查閱資料得子函數(shù)PLU1，其作用是將矩陣A分解成L乘以U的形式。PLU2為調(diào)用PLU1解題的程序，是自己編的）精品文檔放心下載(Ⅰ).function[l,u,p]=PLU1(A) %定義子函數(shù)，其功能為列主元三角分解系數(shù)矩感謝閱讀陣A[m,n]=size(A); %判斷系數(shù)矩陣是否為方感謝閱讀陣ifm~=nerror('矩陣不是方陣')returnendifdet(A)==0%判斷系數(shù)矩陣能否被三角分解謝謝閱讀error('矩陣不能被三角分解')endu=A;p=eye(m);l=eye(m); %將系數(shù)矩陣三角分解，分別求出謝謝閱讀P,L,Ufori=1:mforj=i:mt(j)=u(j,i);fork=1:i-1,.t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i));forj=i+1:mifb<abs(t(j))b=abs(t(j));a=j;endendifa~=iforj=1:mc=u(i,j);u(i,j)=u(a,j);u(a,j)=c;endforj=1:mc=p(i,j);p(i,j)=p(a,j);p(a,j)=c;endc=t(a);,.t(a)=t(i);t(i)=c;endu(i,i)=t(i);forj=i+1:mu(j,i)=t(j)/t(i);endforj=i+1:mfork=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);謝謝閱讀endendendl=tril(u,-1)+eye(m);u=triu(u,0)(Ⅱ).functionx=PLU2(A,b) %定義列主元三角分解法的精品文檔放心下載函數(shù)[l,u,p]=PLU1(A) %調(diào)用PLU分解系數(shù)矩謝謝閱讀陣A,.m=length(A); %由于A左乘p，故b也要左感謝閱讀乘pv=b;forq=1:mb(q)=sum(p(q,1:m)*v(1:m,1));精品文檔放心下載endb(1)=b(1) %求解方程Ly=bfori=2:1:mb(i)=(b(i)-sum(l(i,1:i-1)*b(1:i-1)));精品文檔放心下載endb(m)=b(m)/u(m,m); %求解方程Ux=yfori=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)))/u(i,i);謝謝閱讀endclearx;disp('AX=b的解x