matlab在科學計算中的應用 線性方程組的直接求解-分析方法ppt課件

1、4.2.3 線性方程組的直接求解分析方法 矩陣的三角分解： LU分111212122212nnFnnnnuuuluuALUIllu2na2111121111111111 iiiikkiiikiiiiikkiiiiiikkiiikkiiiikkiiikaluual uual uual uu格式L是一個單位下三角矩陣，u是一個上三角矩陣 l,u=lu(A) lP-1 Ll = 0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000例：用兩種方法對A進行LU

2、分解123241467A A=1 2 3; 2 4 1; 4 6 7; l,u,p=lu(A)l = 1.0000 0 0 0.5000 1.0000 0 0.2500 0.5000 1.0000u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000p = 0 0 1 0 1 0 1 0 0QR分解 將矩陣A分解成一個正交矩陣與一個上三角矩陣的乘積。格式： Q,R = qr(A) 求得正交矩陣Q和上三角陣R，Q和R滿足A=QR。 例： A = 1 2 3;4 5 6; 7 8 9; 10 11 12; Q,R = qr(A) Q = -0.07

3、76 -0.8331 0.5456 -0.0478 -0.3105 -0.4512 -0.6919 0.4704 -0.5433 -0.0694 -0.2531 -0.7975 -0.7762 0.3124 0.3994 0.3748 R = -12.8841 -14.5916 -16.2992 0 -1.0413 -2.0826 0 0 -0.0000 0 0 0 Cholesky(喬里斯基 )分解 格式： D=chol(A)例：進行Cholesky分解。 A=16 4 8; 4 5 -4; 8 -4 22; D=chol(A)D = 4 1 2 0 2 -3 0 0 31648454842

4、2A（1LU分解： A*X=b 變成 L*U*X=b所以 X=U(Lb) 這樣可以大大提高運算速度。例：求方程組 的一個特解。解： A=4 2 -1;3 -1 2;11 3 -1; B=2 10 8; D=det(A)D = 1012312312342232101138xxxxxxxxx利用矩陣的LU、QR和cholesy分解求方程組的解L,U=lu(A)L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0U = 11.0000 3.0000 -1.0000 0 -1.8182 2.2727 0 0 0.5000 X=U(LB)X = 0.4

5、000 3.2000 6.0000 A*Xans = 2.0000 10.0000 8.0000（2Cholesky分解 若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉置的乘積，方程 A*X=b 變成 R*R*X=b所以 X=R(Rb) （3QR分解 對于任何長方矩陣A，都可以進行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A=QR方程 A*X=b 變形成 QRX=b所以 X=R(Qb) 這三種分解，在求解大型方程組時很有用。其優(yōu)點是運算速度快、可以節(jié)省磁盤空間、節(jié)省內存。 三個變換 在線性方程組的迭代求解中，要用到系數矩陣A的上三角矩陣、對角陣和

6、下三角矩陣。此三個變換在MATLAB中可由以下函數實現。 上三角變換： 格式 triu(A,1) 對角變換： 格式 diag(A) 下三角變換： 格式 tril(A,-1) 例：對此矩陣做三種變換。122111221A A=1 2 -2;1 1 1;2 2 1; triu(A,1)ans = 0 2 -2 0 0 1 0 0 0 tril(A,-1)ans = 0 0 0 1 0 0 2 2 0 b=diag(A); bans = 1 1 1122111221A4.3 線性方程組的迭代解法 迭代法的一般形式 線性方程組：12(0)1)( ) ()( ,) ., 0,1,2,.TijnknkAx

7、bAabb bbxMxgxxkxMg其中為n階非奇異方陣，構造同解方程組： 由此得迭代公式： 其為任取的初始向量中( )kkx當 充分大時，即為原方程的近似解.4.3.1 Jacobi迭代法11 11221121 1222221 122112211221 1 det( )det()0,0(1,2,., ) b nnnnnnnnnnijiinna xa xa xba xa xa xba xa xa xbAaainxxb xfxb x方程組：其中不妨設ijij22iiiiii1 12211b() , nnnnnnnnnaijb xfabfaxb xb xbxf 方程組 Ax=b 可寫成 x=Bx+

8、f 由此可構造迭代法： x(k+1)=Bx(k)+f 其中：B=I-D-1A=D-1(L+U) , f=D-1b. D=diag(A), L,U分別為-A的嚴格下三角和嚴格上三角矩陣121311121232121131000nnnnnnnnnbbbbbbbbBbbbbfunction y=jacobi(a,b,x0)D=diag(diag(a); U=-triu(a,1); L=-tril(a,-1);B=D(L+U); f=Db;y=B*x0+f;n=1;while norm(y-x0)=1.0e-6 x0=y; y=B*x0+f; n=n+1;endn 例：用Jacobi方法求解， 設x(

9、0)=0,精度為10-6。 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; jacobi(a,b,0;0;0) n = 11 ans = 0.9958 0.9579 0.7916121232310910272106xxxxxxx4.3.2 Gauss-Seidel迭代法(1)( )kkJacobixBxf迭代：1( )( )12211( )( )221 122( )( )( )1 12211 b , nkknnkknnkkknnnnnnxxb xfxb xb xfxb xb xbxf(k+1)(k+1)(k+1)(1)kx循環(huán)計算過程中存有兩個向量：( )kx

10、1( )( )12211( )221 122(1)(1)(11)(1)12211 b nkknnkkkknnnnnnnknxxb xfxb xb xfxb xb xbxf(k+1)(k+1)(k+1)迭代公式改寫為：(1)(1)( )-1-1 LULUB. L=D L, U=D Ukkkxxxf其中 和 分別為矩陣 的嚴格下三角和嚴格上三角矩陣(1)( )(1)-1( )-1-1 (I-L)U(I-L) U(I-L)%det(I-L)1,Dkkkkxxfxxffb-1-1-1-1-1-1-1(1)( )-1-1-11 (I-L) U=(I-D L) D U=(I-L) D(I-DGf(D-L)

11、 Uf=(D-L) bL) DkkxxGbb由此得原方程 Axb 求解的迭代方程：Ddiag(A)L和U分別為矩陣-A的嚴格下三角和嚴格上三角矩陣其中：function y=seidel(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G=(D-L)U ;f=(D-L)b;y=G*x0+f; n=1;while norm(y-x0)=1.0e-6 x0=y; y=G*x0+f; n=n+1;endn 例：對上例用Gauss-Seidel迭代法求解 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; seidel(a

12、,b,0;0;0) n = 7 ans = 0.9958 0.9579 0.7916121232310910272106xxxxxxxJacobinans.迭代： 11 09958 09579 07916例：分別用Jacobi和G-S法迭代求解,看是否收斂。123122911172216xxx Jacobi迭代法： a=1 2 -2; 1 1 1; 2 2 1; b=9; 7; 6; jacobi(a,b,0;0;0)n = 4ans = -27 26 8 seidel(a,b,0;0;0)n = 1011ans = 1.0e+305 * -Inf Inf -1.7556Seidel迭代法：對

13、有些問題Gauss-Seidel迭代比Jacobi迭代收斂快對有些問題Gauss-Seidel迭代比Jacobi迭代收斂的慢對有些問題Jacobi迭代收斂但Gauss-Seidel迭代發(fā)散4.3.3 SOR迭代法 在很多情況下，J法和G-S法收斂較慢，所以考慮對G-S法進行改進。于是引入一種新的迭代法逐次超松弛迭代法Succesise Over-Relaxation),記為SQR法。 迭代公式為： X(k+1)= (D-wL)-1(1-w)D+wU)x(k) + w(D-wL)-1 b 其中：w最佳值在1, 2之間，不易計算得到，因此 w通常有經驗給出。function y=sor(a,b,w

14、,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);M=(D-w*L)(1-w)*D+w*U); f=(D-w*L)b*w;y=M*x0+f; n=1;while norm(y-x0)=1.0e-6 x0=y; y=M*x0+f; n=n+1;endn例：上例中，當w=1.103時，用SOR法求解原方程。 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; sor(a,b,1.103,0;0;0)n = 8ans = 0.9958 0.9579 0.79164.3.4 兩步迭代法 當線性方程系數矩陣為對稱正定時，可用一種特殊

15、的迭代法來解決，其迭代公式為： （D-Lx(k+1/2) =U x(k) +b （D-Ux(k+1)=Lx(k+1/2) +b= x(k+1/2) =(D-L)-1 U x(k) + (D-L)-1 b x(k+1)= (D-U)-1 Lx(k+1/2) + (D-U)-1 bfunction y=twostp(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G1=(D-L)U; f1=(D-L)b;G2=(D-U)L; f1=(D-U)b;y=G1*x0+f1; y=G2*y+f2; n=1;while norm(y-x0)=1.0e-6 x

16、0=y; y=G1*x0+f1; y=G2*y+f2; n=n+1;endn 例：求解方程組 a=10 -1 2 0; -1 11 -1 3; 2 -1 10 3; 0 3 -1 8; b=6; 25; -11; 15; twostp(a, b, 0; 0; 0; 0) n = 7 ans = 1.0791 1.9824 -1.4044 0.9560123410120611113252110111031815xxxx4.4 線性方程組的符號解法 在MATLAB的Symbolic Toolbox中提供了線性方程的符號求解函數，如 linsolve(A,b) 等同于 X = sym(A)sym(b

17、). solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN ) 例： A=sym(10,-1,0;-1,10,-2;0,-2,10)； b=(9; 7; 6)； linsolve(A,b) ans = 473/475 91/95 376/475 vpa(ans) ans = .99578947368421052631578947368421 .95789473684210526315789473684211 .79157894736842105263157894736842例： x,y = solve(x2 + x*y + y = 3,x2 - 4*x + 3 = 0

18、,x,y) x = 1 3 y = 1 -3/2 4.5 稀疏矩陣技術 稀疏矩陣的建立： 格式 S=sparse(i,j,s,m,n) 生成一mxn階的稀疏矩陣，以向量i和j為坐標的位置上對應元素值為s。 例： n=5; a1=sparse(1:n, 1:n, 4*ones(1,n), n, n) a1 = (1,1) 4 (2,2) 4 (3,3) 4 (4,4) 4 (5,5) 4 例： a2=sparse(2:n, 1:n-1,ones(1,n-1),n,n) a2 = (2,1) 1 (3,2) 1 (4,3) 1 (5,4) 1 full(a2) ans = 0 0 0 0 0 1

19、0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 例：n=5,建立主對角線上元素為4，兩條次對角線為1的三對角陣。 n=5; a1=sparse(1:n,1:n,4*ones(1,n),n,n); a2=sparse(2:n,1:n-1,ones(1,n-1),n,n); a=a1+a2+a2 a = (1,1) 4 (2,1) 1 (1,2) 1 (2,2) 4 (3,2) 1 (2,3) 1 (3,3) 4 (4,3) 1 (3,4) 1 (4,4) 4 (5,4) 1 (4,5) 1 (5,5) 4 full(a)ans = 4 1 0 0 0 1 4 1 0

20、0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 格式 A=spdiags(B,d,m,n) 生成一mxn階的稀疏矩陣，使得B的列放在由d指定的位置。 例： n=5 b=spdiags(ones(n,1),4*ones(n,1),ones(n,1), -1,0,1,n,n); full(b) ans = 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4 格式： spconvert(dd) 對于無規(guī)律的稀疏矩陣，可使用此命令由外部數據轉化為稀疏矩陣。 調用形式為：先用load函數加載以行表示對應位置和元素值的.dat文本文件，再用此

21、命令轉化為稀疏矩陣。 例：無規(guī)律稀疏矩陣的建立。 首先編制文本文件sp.dat如下： 5 1 5.00 3 5 8.00 4 4 2.00 5 5 0 load sp.dat spconvert(sp)ans = (5,1) 5 (4,4) 2 (3,5) 8 full(ans)ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0 0 2 0 5 0 0 0 0 稀疏矩陣的計算： 同滿矩陣比較，稀疏矩陣在算法上有很大的不同。具體表現在存儲空間減少，計算時間減少。 例：比較求解下面方程組n1000時兩種方法的差別。124111411141nxxx n=1000; a1=

22、sparse(1:n,1:n,4*ones(1,n),n,n); a2=sparse(2:n,1:n-1,ones(1,n-1),n,n); a=a1+a2+a2; b=ones(1000,1); tic; x=ab; t1=toct1 = 0.4800 a=full(a); tic; x=ab; t2=toct2 = 1.32204.6 矩陣的特征值問題 4.6.1一般矩陣的特征值與特征向量 格式： d=eig (A) 只求解特征值。 格式： V, D=eig (A) 求解特征值和特征向量。 例： 直接求解： A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15

23、 1; eig(A) ans = 34.0000 8.9443 -8.9443 0.0000精確解： eig(sym(A)ans = 0 34 4*5(1/2) -4*5(1/2)高精度數值解： vpa(ans,70)ans = 0 34. 8.944271909999158785636694674925104941762473438446102897083588981642084 -8.9442719099991587856366946749251049417624734384461028 97083588981642084 同時求出特征值與特征向量： 直接求解： v, d = eig(A) v = -0.5000 -0.8236 0.3764 -0.2236 -0.5000 0.4236 0.0236 -0.6708 -0.5000 0.0236 0.4236 0.6708 -0.5000 0.3764 -0.8236 0.2236 d = 34.0000 0 0 0 0 8.