Matlab實訓8-線性代數(shù)問題的求解.ppt

1、實訓九 線性代數(shù)問題求解,矩陣 線性方程組的直接解法 線性方程組的迭代法 線性方程組的符號解法 稀疏矩陣技術(shù) 特征值與特征向量,9.1 矩陣9.1.1特殊矩陣的輸入,數(shù)值矩陣的輸入 零矩陣、幺矩陣及單位矩陣 生成nn方陣： A=zeros(n), B=ones(n), C=eye(n) 生成mn矩陣： A=zeros(m,n), B=ones(m,n), C=eye(m,n) 生成和矩陣B同樣位數(shù)的矩陣： A=zeros(size(B),隨機元素矩陣 若矩陣隨機元素滿足0,1區(qū)間上的均勻分布 生成nm階標準均勻分布為隨機數(shù)矩陣： A=rand(n,m) 生成nn階標準均勻分布為隨機數(shù)方陣： A

2、=rand(n),對角元素矩陣 已知向量生成對角矩陣： A=diag(V) 已知矩陣提取對角元素列向量： Vdiag(A) 生成主對角線上第k條對角線為V的矩陣： A=diag(V,k),例：diag( )函數(shù)的不同調(diào)用格式 C=1 2 3; V=diag(C) % 生成對角矩陣 V = 1 0 0 0 2 0 0 0 3 V1=diag(V) % 將列向量通過轉(zhuǎn)置變換成行向量 V1 = 1 2 3 C=1 2 3; V=diag(C,2) % 主對角線上第 k條對角線為C的矩陣 V = 0 0 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0,生成三

3、對角矩陣: V=diag(1 2 3 4)+diag(2 3 4,1)+diag(5 4 3,-1) V = 1 2 0 0 5 2 3 0 0 4 3 4 0 0 3 4,Hilbert矩陣及逆Hilbert矩陣 生成n階的Hilbert矩陣： A=hilb(n) 求取逆Hilbert矩陣： B=invhilb(n),Hankel(漢克 ) 矩陣 其中：第一列的各個元素定義為C向量，最后一行各個元素定義為R。H為對稱陣。 H1=hankel(C) 由 Hankel 矩陣反對角線上元素相等得出一下三角陣均為零的Hankel 矩陣,Vandermonde(范德蒙)矩陣,伴隨矩陣 其中：P(s)為

4、首項系數(shù)為1的多項式。,符號矩陣的輸入 數(shù)值矩陣A轉(zhuǎn)換成符號矩陣： B=sym(A) 例： A=hilb(3) A = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 B=sym(A) B = 1, 1/2, 1/3 1/2, 1/3, 1/4 1/3, 1/4, 1/5,9.1.2 矩陣基本概念與性質(zhì),行列式 格式 ：d=det(A) 例：求行列式 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; det(A) ans = 0,例： tic, A=sym(hilb(20);

5、det(A), toc ans = 1/2377454716768534509091644243427616440175419837753486493033185331234419759310644585187585766816573773440565759867265558971765638419710793303386582324149811241023554489166154717809635257797836800000000000000000000000000000000000 elapsed_time = 2.3140 高階的Hilbert矩陣是接近奇異的矩陣。,矩陣的跡 格式：

6、t=trace(A) 矩陣的秩 格式：r=rank(A) 用默認的精度求數(shù)值秩 r=rank(A， ) 給定精度下求數(shù)值秩 矩陣的秩也表示該矩陣中行列式不等于0的子式的最大階次?？勺C行秩和列秩（線性無關(guān)的）應相等。,例 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; rank(A) ans = 3 該矩陣的秩為3，小于矩陣的階次，故為非滿秩矩陣。 例 H=hilb(20); rank(H) 數(shù)值方法 ans = 13 H=sym(hilb(20); rank(H) % 解析方法，原矩陣為非奇異矩陣 ans = 20,矩陣范數(shù),矩陣的范數(shù)定義： 格式：

7、N=norm(A) 求解默認的2范數(shù) N=norm(A，選項) 選項可為1,2,inf等,例：求一向量、矩陣的范數(shù) a=16 2 3 13; norm(a), norm(a,2), norm(a,1), norm(a,Inf) ans = 2.092844953645635e+001 2.092844953645635e+001 3.400000000000000e+001 1.600000000000000e+001 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; norm(A), norm(A,2), norm(A,1), norm(A,Inf)

8、ans = 34 34 34 34 符號運算工具箱未提供norm( )函數(shù)，需先用double( )函數(shù)轉(zhuǎn)換成雙精度數(shù)值矩陣，再調(diào)用norm( )函數(shù)。,特征多項式 格式： C=poly(A) 例： A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; poly(A) 直接求取 ans = 1.000000000000000e+000 -3.399999999999999e+001 -7.999999999999986e+001 2.719999999999999e+003 -2.819840539024018e-012 A=sym(A); poly(A)

9、運用符號工具箱 ans = x4-34*x3-80*x2+2720*x,矩陣多項式的求解,符號多項式與數(shù)值多項式的轉(zhuǎn)換 格式： f=poly2sym(P) 或 f=poly2sym(P,x) 格式： P=sym2poly(f),例： P=1 2 3 4 5 6; % 先由系數(shù)按降冪順序排列表示多項式 f=poly2sym(P,v) % 以 v 為算子表示多項式 f = v5+2*v4+3*v3+4*v2+5*v+6 P=sym2poly(f) P = 1 2 3 4 5 6,矩陣的逆矩陣 格式： C=inv(A) 例： format long; H=hilb(4); H1=inv(H) H1

10、= 1.0e+003 * 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999 -0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984 0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961 -0.13999999999999 1.67999999999984 -4.19999999999961 2.79999999999974,檢驗： H*H

11、1 ans = 1.00000000000001 0.00000000000023 -0.00000000000045 0.00000000000023 0.00000000000001 1.00000000000011 -0.00000000000011 0.00000000000011 0.00000000000001 0 1.00000000000011 0 0.00000000000000 0.00000000000011 -0.00000000000011 1.00000000000011 計算誤差范數(shù)： norm(H*inv(H)-eye(size(H) ans = 6.23579

12、8190375727e-013 H2=invhilb(4); norm(H*H2-eye(size(H) ans = 5.684341886080802e-014, H=hilb(10); H1=inv(H); norm(H*H1-eye(size(H) ans = 0.00264500826202 H2=invhilb(10); norm(H*H2-eye(size(H) ans = 1.612897415528547e-005 H=hilb(13); H1=inv(H); norm(H*H1-eye(size(H) Warning: Matrix is close to singular

13、or badly scaled. Results may be inaccurate. RCOND = 2.339949e-018. ans = 53.23696008570294 H2=invhilb(13); norm(H*H2-eye(size(H) ans = 11.37062973181391 對接近于奇異矩陣，高階一般不建議用inv( )，可用符號工具箱。, H=sym(hilb(7); inv(H) ans = 49, -1176, 8820, -29400, 48510, -38808, 12012 -1176, 37632, -317520, 1128960, -194040

14、0, 1596672, -504504 8820, -317520, 2857680, -10584000, 18711000, -15717240, 5045040 -29400, 1128960, -10584000, 40320000, -72765000, 62092800, -20180160 48510, -1940400, 18711000, -72765000, 133402500, -115259760, 37837800 -38808, 1596672, -15717240, 62092800, -115259760, 100590336, -33297264 12012,

15、 -504504, 5045040, -20180160, 37837800, -33297264, 11099088 H=sym(hilb(30); norm(double(H*inv(H)-eye(size(H) ans = 0,例：奇異陣求逆 A=16 2 3 13; 5 11 10 8; 9 7 6 12; 4 14 15 1; format long; B = inv(A) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. B

16、= 1.0e+014 * 0.93824992236885 2.81474976710656 -2.81474976710656 -0.93824992236885 2.81474976710656 8.44424930131968 -8.44424930131968 -2.81474976710656 -2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656 -0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885 norm(A*

17、B-eye(size(A) 檢驗 ans = 1.64081513306419 A=sym(A); inv(A) 奇異矩陣不存在一個相應的逆矩陣，用符號工具箱的函數(shù)也不行 ? Error using = sym/inv Error, (in inverse) singular matrix,同樣適用于含有變量的矩陣求逆。 例： syms a1 a2 a3 a4; C=a1 a2;a3 a4； inv(C) ans = -a4/(-a1*a4+a2*a3), a2/(-a1*a4+a2*a3) a3/(-a1*a4+a2*a3), -a1/(-a1*a4+a2*a3),矩陣的相似變換與正交矩陣

18、其中：A為一方陣，B矩陣非奇異。 相似變換后，X矩陣的秩、跡、行列式與特征值等均不發(fā)生變化，其值與A矩陣完全一致。 對于一類特殊的相似變換滿足如下條件，稱為正交基矩陣。,例： A=5,9,8,3; 0,3,2,4; 2,3,5,9; 3,4,5,8; Q=orth(A) Q = -0.6197 0.7738 -0.0262 -0.1286 -0.2548 -0.1551 0.9490 0.1017 -0.5198 -0.5298 -0.1563 -0.6517 -0.5300 -0.3106 -0.2725 0.7406 norm(Q*Q-eye(4) ans = 4.6395e-016 no

19、rm(Q*Q-eye(4) ans = 4.9270e-016,例： A=16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1; Q=orth(A) A為奇異矩陣，故得出的Q為長方形矩陣 Q = -0.5000 0.6708 0.5000 -0.5000 -0.2236 -0.5000 -0.5000 0.2236 -0.5000 -0.5000 -0.6708 0.5000 norm(Q*Q-eye(3) ans = 1.0140e-015,9.2 線性方程組直接解法9.2.1線性方程組直接求解矩陣除法,關(guān)于線性方程組的直接解法，如Gauss消去法、選主元消去法

20、、平方根法、追趕法等等，在MATLAB中，只需用“”或“”就解決問題。它內(nèi)部實際包含著許許多多的自適應算法，如對超定方程用最小二乘法，對欠定方程時它將給出范數(shù)最小的一個解，解三對角陣方程組時用追趕法等等。 格式： x=Ab,例：解方程組 A=.4096,.1234,.3678,.2943;.2246,.3872,.4015,.1129; .3645,.1920,.3781,.0643;.1784,.4002,.2786,.3927; b=0.4043 0.1550 0.4240 -0.2557; x=Ab; x ans = -0.1819 -1.6630 2.2172 -0.4467,9.2.

21、2線性方程組直接求解判定求解,例： A=1 2 3 4; 4 3 2 1; 1 3 2 4; 4 1 3 2; B=5 1; 4 2; 3 3; 2 4; C=A B; rank(A), rank(C) ans = 4 ans = 4 x=inv(A)*B x = -1.8000 2.4000 1.8667 -1.2667 3.8667 -3.2667 -2.1333 2.7333,檢驗 norm(A*x-B) ans = 7.4738e-015 精確解 x1=inv(sym(A)*B x1 = -9/5, 12/5 28/15, -19/15 58/15, -49/15 -32/15, 41

22、/15 檢驗 norm(double(A*x1-B) ans = 0,原方程組對應的齊次方程組的解 求取A矩陣的化零矩陣： 格式： Z=null(A) 求取A矩陣的化零矩陣的規(guī)范形式： 格式： Z=null(A, r ),例： 判斷可解性 A=1 2 3 4; 2 2 1 1; 2 4 6 8; 4 4 2 2; B=1;3;2;6; C=A B; rank(A), rank(C) ans = 2 2 Z=null(A,r) % 解出規(guī)范化的化零空間 Z = 2.0000 3.0000 -2.5000 -3.5000 1.0000 0 0 1.0000, x0=pinv(A)*B % 得出一個

23、特解 x0 = 0.9542 0.7328 %全部解 -0.0763 -0.2977 驗證得出的解 a1=randn(1); a2=rand(1); % 取不同分布的隨機數(shù) x=a1*Z(:,1)+a2*Z(:,2)+x0; norm(A*x-B) ans = 4.4409e-015,解析解 Z=null(sym(A) Z = 2, 3 -5/2, -7/2 1, 0 0, 1 x0=sym(pinv(A)*B) x0 = 125/131 96/131 -10/131 -39/131,驗證得出的解 a1=randn(1); a2=rand(1); % 取不同分布的隨機數(shù) x=a1*Z(:,1)

24、+a2*Z(:,2)+x0; norm(double(A*x-B) ans = 0 通解 syms a1 a2; x=a1*Z(:,1)+a2*Z(:,2)+x0 x = 2*a1+3*a2+125/131 -5/2*a1-7/2*a2+96/131 a1-10/131 a2-39/131,摩爾彭羅斯廣義逆求解出的方程最小二乘解不滿足原始代數(shù)方程。,9.2.3 線性方程組的直接求解分析,LU分解,格式 l,u,p=lu(A) L是一個單位下三角矩陣，u是一個上三角矩陣， p是代表選主元的置換矩陣。 故：Ax=y = PAx=Py = LUx=Py = PA=LU l,u=lu(A) 其中l(wèi)等于

25、P-1 L，u等于U，所以(P-1 L)U=A,例：對A進行LU分解 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.0000 u = 4.0000 6.0000 7.0000 0 1.0000 -2.5000 0 0 2.5000 p = 0 0 1 0 1 0 1 0 0, l,u=lu(A) lP-1 L l = 0.2500 0.5000 1.0000 0.5000 1.0000 0 1.0000 0 0 u = 4.0000 6.0000 7.0000 0 1.0000

26、 -2.5000 0 0 2.5000,QR分解 將矩陣A分解成一個正交矩陣與一個上三角矩陣的乘積。 求得正交矩陣Q和上三角陣R，Q和R滿足A=QR。,格式： Q,R = qr(A),例： A = 1 2 3;4 5 6; 7 8 9; 10 11 12; Q,R = qr(A) Q = -0.0776 -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

27、 -1.0413 -2.0826 0 0 -0.0000 0 0 0,Cholesky(喬里斯基 )分解 若矩陣A為 n階對稱正定陣,則存在唯一的對角元素為正的三角陣D，使得,格式： 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 3,利用矩陣的LU、QR和cholesky分解求方程組的解,（1）LU分解： A*X=b 變成 L*U*X=b 所以 X=U(Lb) 這樣可以大大提高運算速度。 例：求方程組 的一個特解。 解： A=4 2 -1;3 -1 2;11 3 0; B=

28、2 10 8; D=det(A) D = 0, L,U=lu(A) L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0 U = 11.0000 3.0000 0 0 -1.8182 2.0000 0 0 0.0000, X=U(LB) Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.018587e-017. X = 1.0e+016 * % 結(jié)果中的警告是由于系數(shù)行列式為零產(chǎn)生的。 -0.4053 % 可以通

29、過A*X驗證其正確性。 1.4862 1.3511 A*X ans = 0 8 8,（2）Cholesky分解 若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置的乘積， 方程 A*X=b 變成 R*R*X=b 所以 X=R(Rb) （3）QR分解 對于任何長方矩陣A，都可以進行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A=QR 方程 A*X=b 變形成 QRX=b 所以 X=R(Qb) 這三種分解，在求解大型方程組時很有用。其優(yōu)點是運算速度快、可以節(jié)省磁盤空間、節(jié)省內(nèi)存。,三個變換 在線性方程組的迭代求解中，要用到系數(shù)矩陣A的上三角矩陣、對角陣和

30、下三角矩陣。此三個變換在MATLAB中可由以下函數(shù)實現(xiàn)。 上三角變換： 格式 triu(A,1) 對角變換： 格式 diag(A) 下三角變換： 格式 tril(A,-1) 例：對此矩陣做三種變換。, 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); b ans = 1 1 1,9.3 迭代解法的幾種形式9.3.1 Jacobi迭代法,方程組 Ax=b A可寫成 A=D-L-U 其中:D=diaga11,a22,ann, -L、-U分別為A

31、的嚴格下、上三角部分（不包括對角線元素）. 由 Ax=b x=Bx+f 由此可構(gòu)造迭代法： x(k+1)=Bx(k)+f 其中：B=D-1(L+U)=I-D-1A, f=D-1b.,function 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; end n,例：用Jacobi方法求解， 設(shè)x(0)=0,精度為10-6。 a=10 -1 0; -1 10 -2; 0

32、 -2 10; b=9; 7; 6; jacobi(a,b,0;0;0) n = 11 ans = 0.9958 0.9579 0.7916,9.3.2 Gauss-Seidel迭代法,由原方程構(gòu)造迭代方程 x(k+1)=G x(k)+f 其中：G=(D-L)-1 U, f=(D-L)-1 b D=diaga11,a22,ann, -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

33、; while norm(y-x0)=1.0e-6 x0=y; y=G*x0+f; n=n+1; end n,例：對上例用Gauss-Seidel迭代法求解 a=10 -1 0; -1 10 -2; 0 -2 10; b=9; 7; 6; seidel(a,b,0;0;0) n = 7 ans = 0.9958 0.9579 0.7916 例：分別用Jacobi和G-S 法迭代求解,看是否收斂。, a=1 2 -2; 1 1 1; 2 2 1; b=9; 7; 6; jacobi(a,b,0;0;0) n = 4 ans = -27 26 8 seidel(a,b,0;0;0) n = 101

34、1 ans = 1.0e+305 * -Inf Inf -1.7556,9.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通常有經(jīng)驗給出。,function y=sor(a,b,w,x0) D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1); M=(D-w*

35、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; end n,例：上例中，當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 = 8 ans = 0.9958 0.9579 0.7916,9.3.4 兩步迭代法,當線性方程系數(shù)矩陣為對稱正定時，可用一種特殊的迭代法來解決，其迭代公式為： （D-L）x(k+1/2) =U x(k) +b （D-U

36、）x(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 b,function 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 x0=y; y=G1*x0+f1; y=G2*y+f2; n=n+1; e

37、nd n,例：求解方程組 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.9560,9.4 線性方程組的符號解法,在MATLAB的Symbolic Toolbox中提供了線性方程的符號求解函數(shù)，如 linsolve(A,b) 等同于 X = sym(A)sym(b). solve(eqn1,eqn2,.,eqnN,var1,var2,.,varN ),例： A=sym(10,-1,0;-1,10

38、,-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,x,y) x = 1 3 y = 1 -3/2,9.5 稀疏矩陣技術(shù),稀疏矩陣的建立： 格式 S=sparse(i,j,s,m,n

39、) 生成一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 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0,例：n=5,建立主對角線上元素為4，兩條次對角線為1的三

40、對角陣。 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 0 0 1 4 1 0 0 0 1 4 1 0 0 0 1 4,格式 A=spdiags(B,d,m,n) 生成一mxn階的稀疏矩

41、陣，使得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ī)律的稀疏矩陣，可使用此命令由外部數(shù)據(jù)轉(zhuǎn)化為稀疏矩陣。 調(diào)用形式為：先用load函數(shù)加載以行表示對應位置和元素值的.dat文本文件，再用此命令轉(zhuǎn)化為稀疏矩陣。 例：無規(guī)律稀疏矩陣的建立。 首先編制文本文件sp.dat如下： 5 1 5.00 3 5 8.00 4 4

42、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,稀疏矩陣的計算： 同滿矩陣比較，稀疏矩陣在算法上有很大的不同。具體表現(xiàn)在存儲空間減少，計算時間減少。 例：比較求解下面方程組n1000時兩種方法的差別。, n=1000; 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; b=ones(1000,1); tic; x=ab; t1=toc t1 = 0.4800 a=full