(完整word版)高斯消元法和三角分解法

1、實驗四 線性方程組的直接數(shù)值解法信息與計算科學(xué)金融 崔振威201002034031一、實驗?zāi)康模壕幊虒崿F(xiàn)高斯消元法和三角分解法二、實驗內(nèi)容：P108.15、P120.1三、實驗要求：1、分別對所給算例利用高斯消元法和三角分解法進(jìn)行數(shù)值求解2、分析系數(shù)矩陣P 108.15的算法穩(wěn)定性主程序高斯消元法:function x,det,flag=Gauss(A,b)%求線性方程組的列主元Gauss消去法%A為方程組的系數(shù)矩陣%b為方程組的右端項%x為方程組的解%det為系數(shù)矩陣A的行列式的值%flag為指標(biāo)向量，flag為failure表示計算失敗，為 ok時表示計算成功n ,m=size(A); n

2、b=le ngth(b);%當(dāng)方程組與列的維數(shù)不相等時停止計算，并輸出出錯信息if n=m');error('方程組與右端項列的維數(shù)不相等，錯誤。return;end%賦初始值，計算flag='OK'det=1;x=zeros( n,1);for k=1: n-1%選主元max1=0;for i=k: nif abs(A(i,k)>max1max1=abs(A(i,k);r=i;endendif max1<1e-10flag='failure'retur n;end%交換兩行if r>kfor j=k: nz=A(k,j);A(

3、k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;end%消元過程for i=k+1: nm=A(i,k)/A(k,k);for j=k+1: nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);end det=det*A( n,n);%回代過程if abs(A( n,n) )<1e-10 flag='failure'retur n;endfor k=n :-1:1for j=k+1: n b(k)=b(k)-A(k,j)*x(j);en

4、dx(k)=b(k)/A(k,k);end x(k)=b(k)/A(k,k); end三角分解法:fun ctio n L,U,flag=LU_Decom(A)%求矩陣A的LU分解%A為被分解的矩陣%L為單位下三角陣%U為單位上三角陣%flag為指標(biāo)向量，flag為failure表示計算失敗，為 ok時表示計算成功n, m=size(A);%要求所分解的矩陣是方陣，否則停止計算，并輸出錯誤信息if n=merror('所分解的矩陣不是一個方陣');return;endL=eye( n);U=zeros (n );flag='OK'for k=1: nfor j=

5、k: nz=0;for q=1:k-1z=z+L(k,q)*U(q,j);endU(k,j)=A(k,j)-z;x =1.0e+002 *endif abs(U(k,k)<e ps flag='failure'retur n;endfor i=k+1: nz=0;for q=1:k-1 z=z+L(i,q)*U(q,k);endL(i,k)=(A(i,k)-z)/U(k,k);end endP108 15解：1、利用高斯消元，在 matlab的命令窗口中輸入命令，并得出結(jié)果:>> A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4

6、1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 0,x,det,flag=Gauss(hilb(4),b)1.000000000000000.500000000000000.333333333333330.250000000000000.500000000000000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000000.166666666666670.250000000000000.200000000000000.1666666666666

7、70.142857142857140.15999999999999-1.199999999999922.39999999999980-1.39999999999987det =1.653439153439300e-007 flag =OK 2、利用三角消元，在 matlab的命令窗口中輸入命令，并得出結(jié)果: >> A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 00',L,U=lu(A),x=U(Lb)1.000000000000000.500000000000000.3333333

8、33333330.250000000000000.500000000000000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000000.166666666666670.250000000000000.200000000000000.166666666666670.142857142857141.000000000000000.500000000000000.333333333333330.250000000000001.000000000000000.500

9、000000000001.000000000000001.000000000000000.333333333333331.000000000000000.250000000000000.90000000000000-0.600000000000001.000000000000000.083333333333330.088888888888890.08333333333333-0.00555555555556-0.008333333333331.0e+002 *0.15999999999999-1.199999999999922.39999999999980-1.39999999999987由上

10、面的輸出結(jié)果可以看出:00.00035714285714通過高斯消元和三角分解所的出的值二者之間的相差很小，甚至從結(jié)果無法判斷二者的相差多少。1、利用高斯消元，在 matlab的命令窗口中輸入命令，并得出結(jié)果: >> A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 0,x,det,flag=Gauss(A,b)1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000.16670.25000.20000.1667

11、0.142916.0000-120.0000 det =1.6534e-007240.0000-140.0000flag =OK 2、利用三角消元，在 matlab的命令窗口中輸入命令，并得出結(jié)果:>> format >> A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 0',L,U=lu(A),x=U(Lb)1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000.16670.25000.

12、20000.16670.14291.00000.50001.00001.00000.33331.00000.25000.9000-0.60001.00000.08330.08890.0833-0.0056-0.0083U =1.00000.50000.33330.25000.000435.000016.0000-120.0000240.0000-140.0000分析：通過三角分解求解線性方程組，當(dāng)矩陣為四位有效數(shù)字表示時和用分?jǐn)?shù)表示時，其輸出的解結(jié)果明顯不同，當(dāng)二者的差值在10-4到10-5時輸出結(jié)果相差值有 1點(diǎn)多（a、b中的三角分解的輸出結(jié)果）。因此可以看出當(dāng)矩陣系數(shù)發(fā)生明顯變化其輸出結(jié)果變化也會明 顯。其穩(wěn)定性不好。P120.1解：1、利用高斯消元，在matlab的命令窗口中輸入命令，并得出結(jié)果:3 5;0 0 2 5;-2 -6 -3 1,b=1 2 3 4,x,det,flag=Gauss(A,b)>> A=1 3 5 7;2 -13-1-2-6-31.34290.6857-3.00001.8000det =flag =OK 2、利用三角消元，在 matlab的命令窗口中輸入命令，并得出結(jié)果: >> A=