1、精選優(yōu)質(zhì)文檔-傾情為你奉上線性代數(shù)方程組的直接解法實(shí)驗(yàn)?zāi)康模壕€性方程組求解的直接法編程實(shí)現(xiàn)，實(shí)驗(yàn)內(nèi)容：線性方程組求解的高斯消去法算法實(shí)現(xiàn) 線性方程組求解的主元素消去法算法實(shí)現(xiàn) 線性方程組求解的LU分解得方法算法實(shí)現(xiàn) 線性方程組求解追趕法算法實(shí)現(xiàn)實(shí)驗(yàn)比較：高斯消去、主元素消去、LU分解都用實(shí)例 這個(gè)進(jìn)行比較。知識(shí)理論解線性方程組的方法大致分為直接法和迭代法。直接法是指假設(shè)計(jì)算過程中不產(chǎn)生舍入誤差，經(jīng)過有限次的運(yùn)算可求的方程組精確解的方法。方程（2-1） 記為：AX=b；一、高斯（Gauss）消元法(1).Gauss消元法是最基本的一種方法。先逐次消去變量，將方程組化成同解的上三角方程組。消元過程

2、：先逐次消去變量，將方程組化成同解的上三角方程組基本思想回代過程：按方程的相反順序求解三角形方程組，得到原方程組的解。(2) Gauss消去法的求解思路為：若先讓第一個(gè)方程組保持不變，利用它消去其余方程組中的,使之變成一個(gè)關(guān)于變?cè)膎-1階方程組。按照（1）中的思路繼續(xù)運(yùn)算得到更為低階的方程組。經(jīng)過n-1步的消元后，得到一個(gè)三角方程。利用求解公式回代得到線性方程組的解。消元過程：第一次消元， （i=1,2,3，n）.將（2-1）中第i個(gè)方程減去第一個(gè)方程乘以（i=1,2,3，n），完成第一次消元， （2-2）其中：；：簡(jiǎn)記為：其中：按上述方法完成n-1次消元后得到。同解的三角方程組： 簡(jiǎn)記為：

3、 回代過程：按逆序逐步回代得到方程的解。 （3）算法：（5）Matlab程序代碼：function RA,RB,n,X=gaus(A,b)B=A,b;n = length(b);RA = rank(A);RB= rank(B);zhicha = RB-RA;if(zhicha=0) disp('因?yàn)镽ARB，所以此方程組無解'); return;endif RA=RB if RA = n disp('因?yàn)镽A=RB=n,所以此方程有唯一解'); X=zeros(n,1); for p = 1:n-1 for k = p+1:n m = B(k,p)/B(p,p)

4、; B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b = B(1:n,n+1); A = B(1:n,1:n); X(n) =b(n)/A(n,n); for q = n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp('因?yàn)镽A=RB<n,所以此方程有無窮解'); endend運(yùn)行檢測(cè)：>>A=1 1 1;12 -3 3;-18 3 -1;b=6;15;-15;RA,RB,n,X=gaus(A,b)因?yàn)镽A=RB=n,所以此方程有唯一解&

5、gt;>RA = 3RB = 3n = 3X = 1.0000 2.00003.0000二、 主元素法1.列主元素法消元(1)基本思想：在每次消元前，在要消去未知數(shù)的系數(shù)中找到絕對(duì)值最大的系數(shù)作為主元，通過對(duì)換行將其換到對(duì)角線上，然后進(jìn)行消元.(2)消元過程： 與Gauss很類似，每次對(duì)對(duì)角線換成最大的值，后面過程與Gauss基本相同。如此經(jīng)過n-1步,（2-1）的增廣矩陣A,b被化為上三角矩陣；回代過程：同Gauss算法一樣回代求解。(3)算法：det<-1對(duì)于k=1,2,3,n-1;|aik,k|=max|aik|(kin)(i)如果aik=0,則停止計(jì)算（det(A)=0）

6、(ii)如果aik=k,則zhuan（i）換行 akj<-> aik,j(j=k,k+1,n) bk<-> bik后邊就是Gauss消元(4).Matlab程序function X=liezhu(A,b)B=A,b;n=length(b);RA=rank(A);RB=rank(B);if RA=RB if RA=n disp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且為n，此方程有唯一解') X=zeros(n,1); for p=1:n-1 y,j=max(abs(B(p:n,p); C=B(p,:); B(p,:)= B(j+p-1,:); B(j+p-1,

7、:)=C; end for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p); B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1); A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); end else disp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且小于n，方程組有無窮解') endelse disp('因?yàn)橄禂?shù)矩陣與增廣矩陣的秩不相同，所以此方程組無解'

8、)endend測(cè)試結(jié)果： >>clear allclcA=1 1 1;12 -3 3;-18 3 -1;b=6 15 -15'liezhu(A,b)>>因?yàn)橄禂?shù)矩陣與增廣矩陣的秩相同且為n，此方程有唯一解ans = 1.0000 2.0000 3.0000例2測(cè)試結(jié)果>> clear allclcA=0.50 1.1 3.1;2.0 4.5 0.36;5.0 0.96 6.5;b=6.0 0.020 0.96'liezhu(A,b)>>ans = -2.6000 1.00002.0000三、直接三角分解法1.高斯消元的矩陣表示 第

9、一次消元 經(jīng)過K次消元后得到： 消元過程：則 （1）LU分解：，因?yàn)锳的各元素已知，可求出L和U的各元素。（2） 看A的第k行，有，再看A的第k列，有。這樣就可以計(jì)算計(jì)算出U的第k行和L的第k列的元素。（3） 從k=1到k=n交替計(jì)算U和L，就能逐步計(jì)算出U和L的全部元素。（4） 分兩步求解LUx=b，先解方程組Ly=b，因?yàn)長(zhǎng)是下三角矩陣，求解只要逐步向前回代。第二步是解方程組Ux=y，因?yàn)閁是上三角矩陣，用逐次向后回代的方法求出x。3.Matlab程序function L,U,Y,X = Doolittle(A,b)n = length(b);L = zeros(n,n);U = zero

10、s(n,n);%先將A進(jìn)行LU分解for i = 1:n %把第一列第一行求出 for j = i:n if i = 1 A(i,j) = A(i,j); if j>i A(j,i) = A(j,i)/A(i,i); end end end %下面計(jì)算U的第i行值 for j = i:n if i= 1 m = 0; for k = 1:i-1 m = m+A(i,k)*A(k,j); end A(i,j) = A(i,j)-m; end end %下面計(jì)算L的第i列值 for j = i:n if i = 1&&j>i l = 0; for k = 1:i-1 l

11、 = l+A(j,k)*A(k,i); end A(j,i) = (A(j,i)-l)/A(i,i); end endend%將LU的值分別加入for i = 1:n for j = i:n U(i,j) = A(i,j); end for j = 1:i if i = j L(i,j) = 1; else L(i,j) = A(i,j); end endend% LY=B，UX=Y，進(jìn)行求解Ab = L,bY = zeros(n,1);Y(1) = Ab(1,n+1)/Ab(1,1);for i = 2:n for j = 1:i-1 Y(i) = Y(i)+Ab(i,j)*Y(j); en

12、d Y(i) = (Ab(i,n+1)-Y(i)/Ab(i,i);endAb = U,Y;X = zeros(n,1);X(n) = Ab(n,n+1)/Ab(n,n);for i = n-1:-1:1 for j = i+1:n X(i) = X(i)+Ab(i,j)*X(j); end X(i) = (Ab(i,n+1)-X(i)/Ab(i,i);endend運(yùn)行檢測(cè)：>>A=1 1 1;12 -3 3;-18 3 -1;b=6;15;-15; L,U,Y,X = Doolittle(A,b)運(yùn)行結(jié)果：>>Ab = 1.0000 0 0 6.0000 12.0000

13、 1.0000 0 15.0000 -18.0000 -1.4000 1.0000 -15.0000L = 1.0000 0 0 12.0000 1.0000 0 -18.0000 -1.4000 1.0000U = 1.0000 1.0000 1.0000 0 -15.0000 -9.0000 0 0 4.4000Y = 6.0000 -57.0000 13.2000X = 1.0000 2.00003.0000四、解三角方程組的追趕法1.設(shè)系數(shù)矩陣為三對(duì)角矩陣則方程組Ax=f稱為三對(duì)角方程組。2.設(shè)矩陣A非奇異，A有Crout分解A=LU，其中L為下三角矩陣，U為單位上三角矩陣，記先依次求

14、出L，U中的元素后，令Ux=y，先求解下三角方程組Ly=f得出y，再求解上三角方程組Ux=y。求解三對(duì)角方程組的2追趕法將矩陣三角分解的計(jì)算與求解兩個(gè)三角方程組的計(jì)算放在一起，使算法吏為緊湊。其計(jì)算公式為：3.算法：1.輸入a=(a2,，an), b=(b1,，bn), c=(c1,，cn-1), d=(d1,，dn),n2.對(duì)i=2,3，nai/bi-1=>ai (li)bi-ci-1ai=>bi (ui)di-aidi-1=>di (yi)3.dn/bn=>dn（xn）4.對(duì)i=n-1,n-2，1(di-cidi+1)/bi=>di (xi)5.輸出d=x,停

15、機(jī)4.Matlab程序：clc;clear;%如果需要驗(yàn)證其他矩陣的解，只需更改A、B和n的值即可A = 2,1,0,0;1,3,1,0;0,1,4,-1;0,0,2,3B = 3;5;4;5n = 4;%獲得A中不為0的數(shù)存到相應(yīng)數(shù)組中for i = 1:n for j = 1:n if i = j b(i) = A(i,j); end if i = j-1 c(j-1) = A(i,j); end if j = i-1 a(i) = A(i,j); end endend%將變換后的值存入相應(yīng)矩陣中u(1) = b(1);for i = 2:n; l(i) = a(i)/u(i-1); u(

16、i) = b(i)-l(i)*c(i-1);endfor i = 1:n for j = 1:n if i = j L(i,j) = 1; U(i,j) = u(i); end if i = j-1 U(i,j) = c(i); end if j = i-1 L(i,j) = l(i); end endend% LY=B，UX=Y，進(jìn)行求解disp('增廣矩陣');AB = L,BB = zeros(n,1);B(1) = AB(1,n+1)/AB(1,1);for i = 2:n for j = 1:i-1 B(i) = B(i)+AB(i,j)*B(j); end B(i)

17、 = (AB(i,n+1)-B(i)/AB(i,i);enddisp('解方程LY=B得');Y = BAB = U,B;B = zeros(n,1);B(n) = AB(n,n+1)/AB(n,n);for i = n-1:-1:1 for j = i+1:n B(i) = B(i)+AB(i,j)*B(j); end B(i) = (AB(i,n+1)-B(i)/AB(i,i);enddisp('方程的解為：');X = B運(yùn)行結(jié)果：A = 2 1 0 0 1 3 1 0 0 1 4 -1 0 0 2 3B = 3 5 4 5增廣矩陣AB = 1.0000 0 0 0 3.0000 0.5000 1.0000 0 0 5.0000 0 0.4000 1.0000 0 4.0000 0 0 0.5556 1.0000 5.0000解方程LY=B得Y = 3.0000 3.5000 2.6000 3.5556方程的解為：X = 1.0000 1.0000 1.00001.0000五、算法比較：Gauss消去法是針對(duì)一般的線性方程組，與線性代數(shù)中的初等變換解線性方