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

《列主元高斯消去法和列主元三角分解法解線性方程》由會員分享，可在線閱讀，更多相關(guān)《列主元高斯消去法和列主元三角分解法解線性方程（13頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、計算方法實驗報告1 【課題名稱】 用列主元高斯消去法和列主元三角分解法解線性方程 【目的和意義】 高斯消去法是一個古老的求解線性方程組的方法，但由它改進得到的選主元的高斯消去法則是目前計算機上常用的解低階稠密矩陣方程組的有效方法。 用高斯消去法解線性方程組的基本思想時用矩陣行的初等變換將系數(shù)矩陣A約化為具有簡單形式的矩陣（上三角矩陣、單位矩陣等），而三角形方程組則可以直接回帶求解 用高斯消去法解線性方程組（其中A∈Rnn）的計算量為：乘除法運算步驟為

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

3、不需要任何中間步驟，從而在計算過程中將高斯消去法的步驟進行了進一步的簡略，大大提高了運算速度，這就是三角分解法 采用選主元的方式與列主元高斯消去法一樣，也是為了避免除數(shù)過小，從而保證了計算的精確度 【計算公式】 1、 列主元高斯消去法 設有線性方程組Ax=b，其中設A為非奇異矩陣。方程組的增廣矩陣為 第1步（k=1）：首先在A的第一列中選取絕對值最大的元素，作為第一步的主元素: 然后交換（A，b）的第1行與第l行元素，再進行消元計算。 設列主元素消去法已經(jīng)完成第1步到第k-1步的按列選主元，交換兩行，消元計算得到與原方程組等價的方

4、程組 A(k)x=b(k) 第k步計算如下： 對于k=1，2，…，n-1 （1）按列選主元：即確定t使 （2）如果t≠k，則交換[A，b]第t行與第k行元素。 （3）消元計算 消元乘數(shù)mik滿足: （4）回代求解 2、 列主元三角分解法 對方程組的增廣矩陣 經(jīng)過k-1步分解后，可變成如下形式： 第k步分解，為了避免用絕對值很小的數(shù)作除數(shù)，引進量 ，于是有=。如果 ，則將矩陣的第t行與第k行元素互

5、換，將（i，j）位置的新元素仍記為或，然后再做第k步分解，這時 【列主元高斯消去法程序流程圖】 【列主元高斯消去法Matlab主程序】 function x=gauss1(A,b,c) %列主元法高斯消去法解線性方程Ax=b if (length(A)~=length(b)) %判斷輸入的方程組是否有誤 disp(輸入方程有誤！) return; end disp(原方程為AX=b：) %顯示方程組 A b disp(

6、------------------------) n=length(A); for k=1:n-1 %找列主元 [p,q]=max(abs(A(k:n,k))); %找出第k列中的最大值，其下標為[p,q] q=q+k-1; %q在A(k:n,k)中的行號轉(zhuǎn)換為在A中的行號 if abs(p)

7、 q>k temp1=A(k,:); %列主元所在行不是當前行，將當前行與列主 A(k,:)=A(q,:); 元所在行交換(包括b) A(q,:)=temp1; temp2=b(k,:); b(k,:)=b(q,:); b(q,:)=temp2; end %消元 for i=k+1:n m(i,k)=A(i,k)/A

8、(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消元處理 end end disp(消元后所得到的上三角陣是) A %顯示消元后的系數(shù)矩陣 b(n)=b(n)/A(n,n); %回代求解 for i=n-

9、1:-1:1 b(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ù)解題】 【列主元三角分解法程序流程圖】 【列主元三角分解法Matlab主程序】 ①自己編的程序： function x=PLU(A,b,eps) %定義函數(shù)列主元三角分解法函數(shù) if (length(A)~=length(b)) %判斷輸入的方程組是否有誤 disp(輸入方程有誤！)

10、 return; end disp(原方程為AX=b：) %顯示方程組 A b disp(------------------------) n=length(A); A=[A b]; %將A與b合并，得到增廣矩陣 for r=1:n if r==1 for i=1:n [c d]=max(abs(A(:,1))); %選取最大列向量，并做行交換 if c<=eps

11、 %最大值小于e，主元太小，程序結(jié)束 break; else end d=d+1-1; p=A(1,:); A(1,:)=A(d,:); A(d,:)=p; A(1,i)=A(1,i); end A(1,2:n)=A(1,2:n); A(2:n,1)=A(2:n,1)/A(1,1);

12、 %求u(1,i) else u(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r); %按照方程求取u(r,i) if abs(u(r,r))<=eps %如果u（r,r）小于e，則交換行 p=A(r,:); A(r,:)=A(r+1,:); A(r+1,:)=p; else end for i=r:n A(

13、r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根據(jù)公式求解，并把結(jié)果存在矩陣A中 end for i=r+1:n A(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r))/A(r,r); %根據(jù)公式求解，并把結(jié)果存在矩陣A中 end end end y(1)=A(1,n+1); for i=2:n h=0; for k=1:i-1 h=h+A(i,k)*y(k);

14、 end y(i)=A(i,n+1)-h; %根據(jù)公式求解y(i) end x(n)=y(n)/A(n,n); for i=n-1:-1:1 h=0; for k=i+1:n h=h+A(i,k)*x(k); end x(i)=(y(i)-h)/A(i,i); %根據(jù)公式求解x(i) end A disp(AX=b的解x是) x=x;

15、 %輸出方程的解 ②可直接得到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ù)矩陣是否為方陣 if m~=n error(矩陣不是方陣) return

16、 end if det(A)==0 %判斷系數(shù)矩陣能否被三角分解 error(矩陣不能被三角分解) end u=A;p=eye(m);l=eye(m); %將系數(shù)矩陣三角分解，分別求出P,L,U for i=1:m for j=i:m t(j)=u(j,i); for k=1:i

17、-1 t(j)=t(j)-u(j,k)*u(k,i); end end a=i;b=abs(t(i)); for j=i+1:m if b

18、 for j=1:m c=u(i,j); u(i,j)=u(a,j); u(a,j)=c; end for j=1:m c=p(i,j); p(i,j)=p(a,j); p(a,j)=c; end

19、 c=t(a); t(a)=t(i); t(i)=c; end u(i,i)=t(i); for j=i+1:m u(j,i)=t(j)/t(i); end for j=i+1:m for k=1:i-1 u(i,j)=u(i,j)-u(i

20、,k)*u(k,j); end end end l=tril(u,-1)+eye(m); u=triu(u,0) (Ⅱ). function x=PLU2(A,b) %定義列主元三角分解法的函數(shù) [l,u,p]=PLU1(A) %調(diào)用PLU分解系數(shù)矩陣A m=length(A);

21、 %由于A左乘p，故b也要左乘p v=b; for q=1:m b(q)=sum(p(q,1:m)*v(1:m,1)); end b(1)=b(1) %求解方程Ly=b for i=2:1:m b(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:1 b(i

22、)=(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ù)解題】 ① ② 【編程疑難】 這是第一次用matlab編程，對matlab的語句還不是非常熟悉，因此在編程過程中，出現(xiàn)了許多錯誤提示。并且此次編程的兩種方法對矩陣的運算也比較復雜。問題主要集中在循環(huán)控制中，循環(huán)次數(shù)多了一次或者缺少了一次，導致數(shù)據(jù)錯誤，一些基本的編程語句在語法上也會由于生疏而產(chǎn)生許多問題，但是語句的錯誤由于系統(tǒng)會提示，比較容易進行修改，數(shù)據(jù)計算過程中的一些邏輯錯誤，比如循環(huán)變量的控制，這些系統(tǒng)不會提示錯誤，需要我們細心去發(fā)現(xiàn)錯誤，不斷修正，調(diào)試。