matlab解線性方程組線性方程組及MATLAB應用
1matlab 解線性方程組 線性方程組及 MATLAB 應用數值實驗 線性方程組與 MATLAB 應用王1．實驗目的：了解矩陣的範數與條件數。 實驗内容：已知矩陣111111A11111111求 A1，A2，A 和 cond2(A)。 11解：編寫了一個 M 檔案來求矩陣 A 的範數與條件數：test3_1.m 如下：A=[1 1 1 1;-1 1 -1 1;-1 -1 1 1;1 -1 -1 1]; norm(A,1) norm(A,2) norm(A,inf) cond(A,2)計算結果依次是： 4 2 4 1.00002．實驗目的：研究高斯消去法的數值穩定性(出現小主2元) 。實驗内容：設方程組 Axb，其中兩個矩陣如下，分别對以上兩個方程組0.310155.291(1)A111.21159.1746.786.13012 ，b11952221159.14371032.09999999999999(2)A251100185.9000000000000162 ，b25511023(1)計算矩陣的條件數，判斷系數矩陣是良态的還是病态的？解： 本題編寫了一個 test3_21 的 M 檔案如下：A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2]; cond(A1) cond(A2)求得兩個矩陣的條件數分别為 68.4296 和 8.9939，易知這矩陣 A1 的條件數遠遠大于 1，而矩陣 A2 的條件數剛大于1，故這,矩陣 A1 為病态矩陣，矩陣 A2 為良态矩陣。 (2)用列主元消去法求得 L 和 U 及解向量 x1,x2R；解：本題利用 Matlab 的列主元三角分解函數 lu()；具體求解如下： >> A1=[0.3*1e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; >> A2=[10 -7 0 1;-3 2.099999999999999 6 2;5 -1 5 -1;0 1 0 2];4>> b1=[59.17;46.78;1;2];>> b2=[8;5.0000000000001;5;1]; >> [L1,U1]=lu(A1)L1 = 0.0000 1.0000 0 0 0.4724 -0.1755 1.0000 0 1.0000 0 0 0 0.0893 0.0202 -0.1738 1.0000 U1 = 11.2000 9.0000 5.0000 2.0000 0 59.1400 3.0000 1.0000 0 0 4-2.8354 1.2307 0 0 0 1.0151 >> [L2,U2]=lu(A2)L2 =1.0000 0 0 0 -0.3000 -0.0000 1.0000 0 0.5000 1.0000 0 0 0 0.4000 -0.3333 1.0000 U2 =10.0000 -7.0000 0 1.0000 0 2.5000 5.0000 -1.5000 0 0 6.0000 2.3000 0 0 0 3.3667 >> y1=L1\b1; >> x1=U1\y1 x1 =3.8457 1.6095 -15.4761 10.4113 >> y2=L2\b2; >> x2=U2\y2 x2 =0.1337 -0.8218 0.8842 0.91091,x2R4x 用不選主元的高斯消去法求得 L 和 U 及解向量；解：編寫一個 LU_Fact 的 M 檔案儲存不選主元的 LU 分解法然後調用求解。具體實作如下LU_Fact.m 如下所示：function [x,L,U,index]=LU_Fact(A,b)%普通 LU 分解。 A 為要分解的矩陣， b 為方程組右端項；%x 為方程組的解，L 為機關下三角陣，U 為上三角陣，index 為訓示變量 [n,m]=size(A); nb=length(b); if n ~= m5error(‘矩陣 A 的行列必須相等！’);return; endif m ~= nberror(‘矩陣 A 的行數必須和 b 的行數相同！’); return; endL=eye(n);U=zeros(n);index=1;x=zeros(n,1);y=zeros(n,1); for k=1:n for j=k:n z=0;for q=1:k-1z=z+L(k,q)*U(q,j); endU(k,j)=A(k,j)-z; endif abs(U(k,k))for i=k+1:n z=0;for q=1:k-1z=z+L(i,q)*U(q,k); endL(i,k)=(A(i,k)-z)/U(k,k); end endy(1)=b(1); for k=2:n z=0;for j=1:k-1z=z+L(k,j)*y(j); endy(k)=b(k)-z; endx(n)=y(n)/U(n,n); for k=n-1:-1:1 z=0;for j=k+1:nz=z+U(k,j)*x(j); endx(k)=(y(k)-z)/U(k,k); end求解過程及結果如下： 對于方程組 1 的求解過程： >> 6[x,L,U,index]=LU_Fact(A1,b1)x = 0 0 0 0L = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1U = 0.0000 59.1400 3.0000 1.0000 0 0 0