1、第三章 線性代數(shù)方程組及矩陣特征值,3.0 預(yù)備知識 3.1 直接法 3.2 迭代法 3.3 不可解問題 3.4 病態(tài)問題, 3.0預(yù)備知識,一、對角陣與三角陣 1、對角陣： diag(A) 提取mn的矩陣A 的主對角線上元素，生 成一個具有min(m,n)個元素的列向量 diag(A,k) 提取第k條對角線的元素 diag(V) 設(shè)V為具有m個元素的向量，將產(chǎn)生一個以向量V的元素為主對角線元素的mm對角矩陣diag(V,k) 產(chǎn)生一個以向量V的元素為第k條對角線的元素的nn(n=m+k)對角陣,2、 矩陣的三角陣 下三角矩陣 tril(A) 提取矩陣A的下三角陣 tril(A,k) 提取矩陣

2、A的第k條對角線以下的元素 上三角矩陣 triu(A) 提取矩陣A的上三角矩陣 triu(A,k) 提取矩陣A的第k條對角線以下的元素,二、 用于專門學(xué)科中的矩陣 (1) 魔方矩陣 魔方矩陣是每行、每列及兩條對角線上的元素和都相等的矩陣。對于n階魔方陣，其元素由1,2,3,n2共n2 個整數(shù)組成. magic(n) 生成一個n階魔方陣,各行各列及兩條 對角線和為(1+2+3+.+n2 )/n,例 magic(5) ans = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 例：將101125等25個數(shù)填入一個5

3、行5列的表格中，使其每行每列及對角線的和均為565。命令如下： M=100+magic(5),(2) 范得蒙矩陣 范得蒙(Vandermonde)矩陣最后一列全為1，倒數(shù)第二列為一個指定的向量，其他各列是其后列與倒數(shù)第二列的點(diǎn)乘積。 vander(V) 生成以向量V為基礎(chǔ)向量的范得蒙矩陣。 例 A=vander(1;2;3;5) A = 1 1 1 1 8 4 2 1 27 9 3 1 125 25 5 1,(3) 希爾伯特矩陣（Hilbert） Hilbert矩陣的每個元素hij=1/(i+j-1) hilb(n) 生成n階希爾伯特矩陣 invhilb(n) 專求n階的希爾伯特矩陣的逆矩陣

4、注意1：高階Hilbert矩陣一般為病態(tài)矩陣，所以直接求逆可能出現(xiàn)錯誤結(jié)論，故應(yīng)該采用invhilb(n) 注意2：由于Hilbert矩陣本身接近奇異，所以建議處理該矩陣時建議盡量采用符號運(yùn)算工具箱，若采用數(shù)值解時應(yīng)該驗(yàn)證結(jié)果的正確性,(4) 托普利茲矩陣 (toeplitz) toeplitz矩陣除第一行第一列外，其他每個元素都與左上角的元素相同。 toeplitz(x,y) 生成一個以x為第一列，y為第一行的托普利茲矩陣。這里x, y均為向量，兩者不必等長。toeplitz(x)用向量x生成一個對稱的托普利茲矩陣。 例 T=toeplitz(1:5,1:7) T = 1 2 3 4 5 6

5、 7 2 1 2 3 4 5 6 3 2 1 2 3 4 5 4 3 2 1 2 3 4 5 4 3 2 1 2 3,(5) 帕斯卡矩陣 二次項(xiàng)(x+y)n展開后的系數(shù)隨n的增大組成一個三角形表，稱為楊輝三角形。由楊輝三角形表組成的矩陣稱為帕斯卡(Pascal)矩陣。 pascal(n) 生成一個n階帕斯卡矩陣。,pascal(6) ans = 1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 56 1 5 15 35 70 126 1 6 21 56 126 252 例 求(x+y)5的展開式。 pascal(6)次對角線上的元素1，5，

6、10，10，5，1即為展開式的系數(shù)。,三、向量和矩陣的范數(shù) norm(V)或norm(V,2) 求向量V（或矩陣V )的2范數(shù) norm(V,1) 求向量V（或矩陣V）的1范數(shù) norm(V,inf) 求向量V（或矩陣V)的范數(shù) 例 已知V，求V的3種范數(shù)。 命令如下： V=-1,1/2,1; v1=norm(V,1) %求V的1范數(shù) v2=norm(V) %求V的2范數(shù) vinf=norm(V,inf) %求V的范數(shù),常用的向量范數(shù):,范數(shù)意義下的單位向量: X=x1, x2T,-1,常用的矩陣范數(shù):,例 求矩陣A的三種范數(shù) 命令如下： A=17,0,1,0,15;23,5,7,14,16;

7、4,0,13,0,22;10,12,19,21,3;11,18,25,2,19; a1=norm(A,1) %求A的1范數(shù) a2=norm(A) %求A的2范數(shù) ainf=norm(A,inf) %求A的范數(shù),四、 矩陣的逆與偽逆 1、 矩陣的逆（后面研究完Gauss消去法后還將給出求逆的方法） 求方陣A的逆可用 inv(A) 注意：該函數(shù)也適用于符號變量構(gòu)成的矩陣的求逆 例 用求逆矩陣的方法解線性方程組。 命令如下： A=1,2,3;1,4,9;1,8,27; b=5,2,6; x=inv(A)*b 一般情況下，用左除x=Ab比求矩陣的逆的方法更有效(因?yàn)锳奇異或接近奇異時，用inv(A)可

8、能出錯),例：Hilbert矩陣（非常有名的病態(tài)矩陣）：,驗(yàn)證從55到1414的Hilbert矩陣病態(tài)特征,clear format long; for n=5:14 H=hilb(n);norm1=norm(H*inv(H)-eye(size(H); H1=invhilb(n);norm2=norm(H*H1-eye(size(H); fprintf( n=%3.0f norm1=%e norm2=%en, n,norm1,norm2) end,n= 5 norm1=1.409442e-011 norm2=3.637979e-012 n= 6 norm1=2.534893e-009 norm

9、2=1.455203e-011 n= 7 norm1=9.810538e-008 norm2=5.208793e-010 n= 8 norm1=3.123671e-006 norm2=4.822187e-008 n= 9 norm1=8.116595e-005 norm2=2.479206e-007 n= 10 norm1=2.645008e-003 norm2=1.612897e-005 n= 11 norm1=7.200720e-002 norm2=1.305122e-003 Warning: Matrix is close to singular or badly scaled. Res

10、ults may be inaccurate. RCOND = 2.632766e-017. n= 12 norm1=1.176913e+001 norm2=5.576679e-002 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.339949e-018. n= 13 norm1=5.323696e+001 norm2=1.137063e+001 Warning: Matrix is close to singular or badly scaled. Res

11、ults may be inaccurate. RCOND = 1.708191e-019. n= 14 norm1=1.224232e+004 norm2=2.836298e+002,說明1:對于Hilbert求逆時，不建議用inv,可用 invhilb直接產(chǎn)生逆矩陣 說明2：符號工具箱中也對符號矩陣定義了inv( )函數(shù),即使對更高階的奇異矩陣也可以精確的求解出逆矩陣,例：H=sym(hilb(30); norm(double(H*inv(H)-eye(size(H) 結(jié)果為 ans = 0,說明3：對于奇異矩陣用數(shù)值方法的inv( )函數(shù)，會給出錯誤的結(jié)果，但符號工具箱中的inv( )會

12、給出正確的結(jié)論,例 A=16,2,3,13;5,11,10,8; 9,7,6,12;4,14,15,1; det(A) B=inv(A),結(jié)果： ans = 0 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. B = 1.0e+014 * 0.93824992236885 2.81474976710656 -2.81474976710656 -0.93824992236885 2.81474976710656 8.44424930131

13、968 -8.44424930131968 -2.81474976710656 -2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656 -0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885,但用符號工具箱的inv:,A=16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1; A=sym(A) inv(A),結(jié)果： ? Error using = sym.inv Error, (in inverse

14、) singular matrix,2、 矩陣的偽逆 pinv(A) 若A不是一個方陣，或A是一個非滿秩的方陣時，矩陣A沒有逆矩陣，但可以找到一個與A的轉(zhuǎn)置矩陣A同型的矩陣B，使得： ABA=A, BAB=B 此時稱矩陣B為矩陣A的偽逆。 例 求A的偽逆，并將結(jié)果送B A=3,1,1,1;1,3,1,1;1,1,3,1; B=pinv(A) 例 求矩陣A的偽逆 A=0,0,0;0,1,0;0,0,1; pinv(A),五、 求方陣A的行列式: det(A) 例 用克萊姆(Cramer)方法求解線性方程組(不建議使用） 程序如下： D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3

15、,-2,2; %定義系數(shù)矩陣 b=4;6;12;6; %定義常數(shù)項(xiàng)向量 D1=b,D(:,2:4); %用方程組的右端向量置換D的第1列 D2=D(:,1),b,D(:,3:4); %用方程組的右端向量置換D的第2列 D3=D(:,1:2),b,D(:,4:4);%用方程組的右端向量置換D的第3列 D4=D(:,1:3),b;%用方程組的右端向量置換D的第4列,DD=det(D); x1=det(D1)/DD; x2=det(D2)/DD; x3=det(D3)/DD; x4=det(D4)/DD; x1,x2,x3,x4,六、 求矩陣A的秩 rank(A) 七、 求矩陣A的跡 trace(A

16、) 例：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; r=rank(D) t=trace(A) 結(jié)果：r = 4 t= 6,九、求矩陣特征多項(xiàng)式、特征值、特征向量,設(shè)A是一個 nn 矩陣，,為A的特征多項(xiàng)式。,特征多項(xiàng)式的根稱為矩陣A的特征值。,c=poly(A) 求矩陣A的特征多項(xiàng)式的系數(shù) roots(c) 求多項(xiàng)式c的根,八、 求矩陣A的共軛矩陣 conj(A),eig(A) 求矩陣A的特征值 常用的調(diào)用格式有 ： E=eig(A) 求矩陣A的全部特征值，構(gòu)成向量E。 V,D=eig(A) 求矩陣A的全部特征值，構(gòu)成對角陣D，并求A的特征向量構(gòu)成V的列向量。

17、 V,D=eig(A,nobalance) 與第2種格式類似，但第2種格式中先對A作相似變換后求矩陣A的特征值和特征向量，而格式3直接求矩陣A的特征值和特征向量，即A中某項(xiàng)非常小，這樣求出的特征值及特征向量更精確。 V,D=eig(A,B) 計算廣義特征值和特征向量，使 AV=BVD。,例：設(shè)矩陣,A=3 4 -2;3 -1 1;2 0 5; E=eig(A) V,D=eig(A) V,D=eig(A,nobalance),現(xiàn)求解線性方程組,情形1：m=n（恰定方程組） 在MATLAB中的求解命令有：,情形2：mn（超定方程），多用于曲線擬合。,解線性方程組的一般函數(shù)文件如下： functio

18、n x,y=line_solution(A,b) m,n=size(A);y=; if norm(b)0 % b非零，非齊次方程組 if rank(A)=rank(A,b) %方程組相容 （有解） if rank(A)=n %有唯一解 x=Ab; else %方程組有無窮多個解,基礎(chǔ)解系 disp(原方程組有有無窮個解，其齊次方程組的基礎(chǔ) 解系為y，特解為x); y=null(A,r); x=Ab; end,else %方程組不相容（無解） ，給出最小二乘解 disp(方程組的最小二乘法解是：); x=Ab; end else %齊次方程組 if rank(A)=n %列滿秩 x=zero(n

19、,1) %0解 else %非0解，無窮多個解 disp(方程組有無窮個解，基礎(chǔ)解系為x); x=null(A,r); end end return,如在MATLAB命令窗口，輸入命令 A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6; x,y=line_solution(A,b) 及： A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2; x,y=line_solution(A,b) 分別顯示其求解結(jié)果。,求解線性方程組（m=n),用克萊姆法則，理論上可行，但實(shí)際運(yùn)算時工作量大，耗時。下面研究,、 一、Gauss簡單消元法（m=

20、n),設(shè) 有,用線性代數(shù)中的克萊姆法則求解時，工作量非常大， 工作量小的方法是 Gauss消元法。,3.1 解線性方程組的直接法,消 元：,以此類推，最后方程組化為：,回 代：,失效，故應(yīng)選主元,二、列主元素消去法-計算結(jié)果可靠,到此原方程組化為,到此原方程組化為,(上三角方程組) （3.2),(n-1) 原方程組化為,以上為消元過程。,(n) 回代求解公式,（3.3) 是回代過程。,（3.3),例1：在MATLAB上，用Gauss列主元消去法求解方程組：,clear a = -0.04 0.04 0.12 3; 0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0 %先

21、定義增廣矩陣 x = 0,0,0; %先將解設(shè)為零向量，后面重新賦值 tempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo; a %第一次選主元（第一行與第二行交換）,a(2,:) = a(2,:) - a(1,:)*a(2,1)/a(1,1); %將第一個對角元下面的數(shù)字消為0 a(3,:) = a(3,:) - a(1,:)*a(3,1)/a(1,1); a,a = 0.5600 -1.5600 0.3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286,tempo = a(3,

22、:); a(3,:) = a(2,:); a(2,:)=tempo;a %第二次選主元 a(3,:) = a(3,:) - a(2,:)*a(3,2)/a(2,2); a %第二次將第二個對角元下的數(shù)字消為0,x(3) = a(3,4)/a(3,3); %消去完成，現(xiàn)在開始回代 x(2) = (a(2,4) - a(2,3)*x(3)/a(2,2); x(1) = (a(1,4) - a(1,2:3)*x(2:3)/a(1,1); x,運(yùn)行得方程組的解為：,a = 0.5600 -1.5600 0.3200 1.0000 0.0000 0.5714 -0.1429 0.4286 0.0000

23、0 0.1250 3.1250,x = 7.0000 7.0000 25.0000,也可直接用 x=Ab,說明： (1)也可采用無回代的列主元消去法(叫Gauss-Jordan消去法)，該法同時消去對角元上下的元素，且仍舊需要選主元，但比有回代的列主元消去法的乘除運(yùn)算次數(shù)多。GaussJordan消去法的優(yōu)點(diǎn)之一是用它來算逆矩陣的算法非常容易解釋。 (2)有回代的列主元消去法所進(jìn)行的乘除運(yùn)算次數(shù)為 ，計算量很小。,例2：用GaussJordan消去法求解上例中的矩陣 的逆矩陣。,clear A = -0.04 0.04 0.12 ; 0.56 -1.56 0.32 ; -0.24 1.24 -

24、0.28 a=A,eye(3);a tempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo; a%第一次選主元，并交換第一行和第二行 a(1,:) = a(1,:)/a(1,1) %將主元標(biāo)準(zhǔn)化 for i=2:3; a(i,:)=a(i,:) - a(i,1)*a(1,:); end; %將主元下的元素消為零 a,tempo = a(3,:); a(3,:) = a(2,:); a(2,:)=tempo; a %第二次選主元，交換第二行和第三行 a(2,:)=a(2,:)/a(2,2); a % 第二個主元標(biāo)準(zhǔn)化,a = 0.5600 -1.5600 0.

25、3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286,for i=1:3 if i=2, a(i,:)=a(i,:)-a(i,2)*a(2,:); end end % a(2,2)的上下元素消為0 a,a = 1.0000 0 -0.1250 0 3.8750 4.8750 0 1.0000 -0.2500 0 0.7500 1.7500 0 0 0.1250 1.0000 0.1250 0.1250,a(3,:)=a(3,:)/a(3,3) for i=1:3; if i=3, a(i,:)=a(i,:)-a(i,3

26、)*a(3,:); end; end; a A_inv = a(:,4:6) A*A_inv,結(jié)果為：A_inv = 1.0000 4.0000 5.0000 2.0000 1.0000 2.0000 8.0000 1.0000 1.0000,也可用rref命令來求,三、 Gauss 全主元消去法： 優(yōu)點(diǎn)-計算結(jié)果更可靠； 缺點(diǎn)-挑主元花時間更多， 次序有變動，程序復(fù)雜。,四、Gauss消元法的應(yīng)用,1、求行列式：det(A),2、求逆矩陣： inv（A）或 A(-1) rref(A,eye(size(A) 將(A,E)化為行最簡形，其實(shí)就是Gauss-Jordon消去法,（3）求解方程組Ax

27、=b,求逆矩陣的思想: x=inv(A)*b或x=A(-1)*b,左除法（原理上是運(yùn)用高斯消元法求解,但Matlab在實(shí)際執(zhí)行過程中是通過LU分解法進(jìn)行的）：x=Ab,符號矩陣法（此法最接近精確值，但運(yùn)算速度慢）x=sym(A)sym(b),化為行最簡形：C=A,b rref(C),例1：在MATLAB上，用Gauss列主元消去法求解方程組：,A= -0.04 0.04 0.12; 0.56 -1.56 0.32; -0.24 1.24 -0.28; b=3;1;0; x11=inv(A)*b %法1：求逆思想 x12=A(-1)*b %法1：求逆思想 x3=Ab %法2：左除法 x4=sym

28、(A)sym(b) %法3：符號矩陣法 C=A,b;rref(C) %法4：化為行最簡形,定義3.1,叫,的三角（因子）分解，其中 是,是上三角。,下三角,若L為單位下三角陣（對角元全為1），U為上三角陣，則稱 A=LU 為Doolittle分解;,若 L 是下三角，U是單位上三角，則稱A=LU為Crout分解。,定理3.1 n 階陣,有唯一Doolittle分解(Crout）,的所有順序主子式不為0.,三角分解不唯一,為此引入,定義3.2,五利用矩陣三角分解法求解方程組,為什么要討論三角分解？若在消元法進(jìn)行前能實(shí)現(xiàn)三角分解,， 則,容易回代求解,在Gauss消去法中，選主元改變了行的次序，盡

29、管對于Gauss消去法來說，這種次序的變換無法事先知道，但是這個變化的影響卻可以用一個算子P表示，其中P是一個置換矩陣。用P左乘原始矩陣A得到： PAx=Py 或 對 做Gauss消去法不需要選主元，所以對 做LU分解，同樣不需要選主元。,實(shí)際上，將選主元Gauss消去法里的行交換同樣作用于單位矩陣，所得矩陣即為P。,在MATLAB中，LU分解的命令是lu，有兩種格式：,l,u,p=lu(A) 其中A是待分解矩陣；l,u,p分別代表L(主對角線元素為1的下三角矩陣）、U（上三角矩陣）和由單位陣變換出的置換矩陣P 滿足：PA=LU，即 A=P-1LU,l,u=lu(A) 其中l(wèi)= P-1L，所以

30、A=lU=P-1LU 用此格式求出來的 l 并不一定是真正的下三角矩陣， 需要換行后才能是真正的下三角矩陣,實(shí)現(xiàn)LU分解后， 若采用l,u=lu(A) ，則Ax=b的解為 x=u (l b) 若采用l，u，p=lu(A)，則Ax=b的解為 x=u( l ( p*b),例 利用LU分解法求解方程組,A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先錄入系數(shù)矩陣 b=9 23 22 47 L,U=lu(A) %對A矩陣做LU分解 y=Lb %求解方程組Ly=bx=Uy %求解方程組Ux=y得到方程組的最終解,故方程組的最終解為: x =(0.5000,2.0000,3

31、.0000,-1.0000),或 A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先錄入系數(shù)矩陣 b=9 23 22 47 L,U,P=lu(A) %對A矩陣做LU分解 y=LP*b %求解方程組Ly=Pbx=Uy %求解方程組Ux=y得到方程組的最終解,X=QR：X為方陣，Q為正交矩陣， R為上三角矩陣,六、利用矩陣QR分解法求解方程組,Q，R=qr(X)：產(chǎn)生一個一個正交矩陣Q和一個上三角矩陣R，滿足X=QR,Q，R，E=qr(X)：產(chǎn)生一個一個正交矩陣Q、一個上三角矩陣R、置換矩陣E，滿足XE=QR,實(shí)現(xiàn)QR分解后， 若采用Q，R=qr(X) ，則Ax=b的

32、解為 x=R(Qb) 若采用Q，R，E=qr(X)，則Ax=b的解為 x=E(R(Qb),例：用QR分解求解線性方程組。 A=3，1，-4，1;1，-3，0，2;0，2，1，-1; 1，6，- 1，-3; b=12，-6，4，0; Q，R=qr(A); x=R(Qb) 或采用QR分解的第2種格式，命令如下： Q，R，E=qr(A); x=E*(R(Qb) 兩種得到結(jié)果都是 x = -16.4444 20.6667 -1.1111 36.2222,七、利用矩陣Cholesky分解法求解方程組,若矩陣X是對稱正定的， X=RR, R為上三角矩陣,R=chol(X)：產(chǎn)生一個上三角陣R，使RR=X

33、若X為非對稱正定，則輸出一個出錯信息,R，p=chol(X)：這個命令格式將不輸出出錯信息。當(dāng)X為對稱正定的，則p=0，R與上述格式得到的結(jié)果相同；否則p為一個正整數(shù)。 如果X為滿秩矩陣，則R為一個階數(shù)為q=p-1的上三角陣，且滿足 RR=X(1:q，1:q),實(shí)現(xiàn)Cholesky分解后，線性方程組Ax=b變成RRx=b，所以x=R(Rb)。,【例】用Cholesky分解求解線性方程組。 A=3，1，-4，1;1，-3，0，2;0，2，1，-1;1，6，-1，-3; b=12，-6，4，0; R=chol(A) 結(jié)果： ? Error using = chol Matrix must be p

34、ositive definite 命令執(zhí)行時，出現(xiàn)錯誤信息，說明A為非正定矩陣。,八、解三對角方程組追趕法,給定方程組,按行嚴(yán)格對角占優(yōu),（三對角方程組）,其求解算法是Gauss消去法的一種變形，稱為三對角法（追趕法）。解法如下：,1、由第一個方程,令,2、將其代入第二個方程 ，得：,再令,3、將其代入第三個方程,得,再令,以此類推，,令,將,代入最后一個方程,令,以上過程稱為 “追”；,總結(jié)“追”的算法：,Step1：（初始化變量）,Step2：,下面開始“趕”：,Step3：先求,Step4：再求其他的,三對角方程組的求解程序 tri_diag.m 如下：,% tri_diag(a,b,c

35、,d,n) solves a tridiagonal equation. function x = tri_diag(a,b,c,d,n) for i=2:n r=a(i)/b(i-1); d(i)=d(i)-r*d(i-1); b(i)=b(i)-r*c(i-1); end d(n)=d(n)/b(n); for i=n-1:-1:1 d(i)=(d(i)-c(i)*d(i+1)/b(i); end x=d;,迭代法適用于解大型稀疏方程組(萬階以上的方程組,系數(shù)矩陣中零元素占很大比例,而非零元按某種模式分布) 背景: 電路分析、邊值問題的數(shù)值解和數(shù)學(xué)物理方程 問題: (1)如何構(gòu)造迭代格式？

36、 (2)迭代格式是否收斂？ (3)收斂速度如何？ (4)如何進(jìn)行誤差估計？,3.2 解線性方程組的迭代法,一、稀疏矩陣存儲方式的產(chǎn)生與轉(zhuǎn)化,1、由完全存儲方式 轉(zhuǎn)為 稀疏存儲方式命令：,B=sparse(A) 將矩陣A轉(zhuǎn)化為稀疏存儲方式的矩陣B sparse(m,n) 生成一個mn的所有元素都是0的稀疏 矩陣 sparse(u,v,S) u,v,S是3個等長的向量。 S是要建立的稀疏矩陣的非0元素； u(i)、v(i)分別是S(i)的行標(biāo)和列標(biāo)； 該函數(shù)生成一個max(u)行、max(v)列并以S為稀疏元素的稀疏矩陣, u=1,1,4; v=1,2,1; s=1,3,7; sparse(u,v

37、,s),結(jié)果：ans = (1,1) 1 (4,1) 7 (1,2) 3,U,V,S=find(A) 返回矩陣A中非0元素的下標(biāo)和元素，這里產(chǎn)生的U、V、S可作sparse(u,v,s)的參數(shù) full(A) 返回稀疏存儲矩陣A對應(yīng)的完全存儲方式矩陣,相關(guān)操作的函數(shù)：,2、 產(chǎn)生一個稀疏矩陣 把要建立的稀疏矩陣的非0元素及其所在行和列的位置表示出來后由MATLAB自己產(chǎn)生其稀疏存儲方式，這需要使用spconvert函數(shù)。調(diào)用格式為：,B=spconvert(A) 其中A為一個m3或m4的矩陣，每行表示一個非0元素，m是非0元素的個數(shù)，A每個元素的意義是：(i,1) 第i個非0元素所在的行。(i

38、,2) 第i個非0元素所在的列。(i,3) 第i個非0元素值的實(shí)部。(i,4) 第i個非0元素值的虛部，若矩陣的全部元素都是實(shí)數(shù)，則無須第四列。, A=1,1,4; 1,2,1; 1,3,7; spconvert(A),結(jié)果： ans = (1,1) 4 (1,2) 1 (1,3) 7,3、單位稀疏矩陣的產(chǎn)生 eye 產(chǎn)生一個完全存儲方式的單位矩陣 speye(m,n) 產(chǎn)生一個mn的稀疏存儲單位矩陣, speye(2,3) ans = (1,1) 1 (2,2) 1,二、簡單迭代法 1.迭代法建立. 考慮 Ax=b,(矩陣B不唯一),對應(yīng)寫出,產(chǎn)生向量序列,若收斂,記,則于(1)兩端取極限有

39、：,上式說明：,是解向量 ，從而當(dāng)k充分大時,注意: 迭代陣B不唯一，B的選取影響收斂性。,解向量,(1)叫簡單迭代法,B叫迭代矩陣。,2.收斂性. 定義 稱,為矩陣B的譜半徑。,定理,定理 對于簡單迭代法，若迭代矩陣,設(shè)有方程組( 其中 ) Ax = b，即,(2),作等價變形,(3),-Jacobi迭代法,(k=0,1,2,),(4),法1：,法2：,Ax = b,A=D+(L+U),Dx = - (L+U ) x + b,= x = - D-1(L+U)x + D-1 b,= x = BJ x+ f,= 迭代公式：x(k+1) = BJ x(k)+ f , (k=0,1,2),BJ= -

40、D-1（L+U）, f =D-1b,編寫實(shí)現(xiàn)Jacobo迭代法的函數(shù)jacobi.m如下:,function s=jacobi(a,b,x0,err) % jacobi迭代法求解線性方程組，a為系數(shù)矩陣，b為%ax=b右邊的%矩陣b,x0為迭代初值，err為迭代誤差 if nargin=3 err=1.0e-6; elseif nargin3 error return end D=diag(diag(a);%構(gòu)造對角陣D D1=inv(D);%求對角陣D的逆矩陣 L=tril(a,-1);%構(gòu)造嚴(yán)格下三角陣 U=triu(a,1);%構(gòu)造嚴(yán)格上三角陣 B=-D1*(L+U);%求出迭代矩陣,f

41、=D1*b; s=B*x0+f; n=1;%n為迭代次數(shù) while norm(s-x0)=err n=n+1 x0=s; s=B*x0+f; s end n,A=9 -1 -1;-1 10 -1;-1 -1 15; b=7;8;13; x0=0;0;0; err=0.00005; s=jacobi(A,b,x0,err),結(jié)果：n 1 0.7778 0.8000 0.8667 2 0.9630 0.9644 0.9719 3 0.9929 0.9935 0.9952 4 0.9987 0.9988 0.9991 5 0.9998 0.9998 0.9998 6 1.0000 1.0000 1

42、.0000 7 1.0000 1.0000 1.0000,例,三、 Gauss-Seidel 迭代法,對于,Ax = b,Jacobi迭代公式為,(k=0,1,2,),(4),與（4）對應(yīng)的迭代公式為：, Gauss-Seidel迭代法,(5),法1：,法2：,Ax = b,A=D+L+U,(D+L) x = - U x + b,= 迭代公式：x(k+1) = BG-S x(k)+ f , (k=0,1,2),編寫實(shí)現(xiàn)Seidel迭代法的函數(shù)seidel.m如下:,function s=seidel(a,b,x0,err) % seidel迭代法求解線性方程組，a為系數(shù)矩陣，b為%ax=b右邊

43、的%矩陣b,x0為迭代初值，err為迭代誤差 if nargin=3 err=1.0e-6; elseif nargin3 error return end D=diag(diag(a);%構(gòu)造對角陣D L=tril(a,-1);%構(gòu)造嚴(yán)格下三角陣 U=triu(a,1);%構(gòu)造嚴(yán)格上三角陣 C=inv(D+L);,B=-C*U; f=C*b; n=1%n為迭代次數(shù) s=B*x0+f while norm(s-x0)=err n=n+1 x0=s; s=B*x0+f; s end n,A=9 -1 -1;-1 10 -1;-1 -1 15; b=7;8;13; x0=0;0;0; err=0.00005; s=seidel(A,b,x0,err),結(jié)果：n 1 0.7778 0.8778 0.