1、第三章第三章 線性方程組及矩陣特征值線性方程組及矩陣特征值 在這一章中我們會學習到線性方程組的解法，在這一章中我們會學習到線性方程組的解法，有有直接求解直接求解和和迭代求解迭代求解兩種方法，線性方程組和兩種方法，線性方程組和矩陣是緊密聯系的，我們先來學習預備知識，有矩陣是緊密聯系的，我們先來學習預備知識，有關矩陣運算的一些關矩陣運算的一些MATLAB命令。命令。3.0 預備知識一、對角陣與三角陣一、對角陣與三角陣1. 對角陣：求矩陣的對角元素用對角陣：求矩陣的對角元素用diag(A)函數，其還有更函數，其還有更進一步的形式進一步的形式diag(A,k)，其功能是提取第，其功能是提取第k條對角線

2、條對角線的元素的元素(注注：k可以小于可以小于0，主對角線為第，主對角線為第0條對角線，條對角線，主對角線上方依次為第主對角線上方依次為第1，第，第2，第，第k條對角線，條對角線，下方依次為第下方依次為第-1，第，第-2，第第-k條對角線條對角線)。(1)提取矩陣的對角線元素提取矩陣的對角線元素 設設A為為mn矩陣，矩陣，diag(A)函數用于提取矩陣函數用于提取矩陣A的的主主對角線元素對角線元素，產生一個具有，產生一個具有min(m,n)個元素的個元素的列向量列向量。 (2)構造對角矩陣構造對角矩陣 設設V為具有為具有m個元素的向量，個元素的向量，diag(V)將產生一個將產生一個m階對角階

3、對角矩陣，其主對角線元素即為向量矩陣，其主對角線元素即為向量V的元素，其它元素為的元素，其它元素為0。 diag(V)函數也有更進一步的形式函數也有更進一步的形式diag(V,k)，其功能是產生一，其功能是產生一個個n階階(n=m+k)對角陣，其第對角陣，其第k k條對角線的元素即為向量條對角線的元素即為向量V的的元素。元素。例例：A=1 2 3;4 5 6diag(A)diag(A,1)V=1 2 3 4;diag(V)diag(V,2)diag(V,-1) 例例1.1 先建立先建立55矩陣矩陣A，然后將，然后將A的第的第1行元素乘行元素乘以以1，第，第2行乘以行乘以2，第，第5行乘以行乘以

4、5。命令如下：命令如下：A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19;%或者或者A=magic(5);D=diag(1,2,3,4,5);D*A 2. 矩陣的三角陣矩陣的三角陣 (1)下三角矩陣下三角矩陣 求矩陣求矩陣A的下三角陣的的下三角陣的MATLAB函數是函數是tril(A)。 tril(A)函數也有更進一步的一種形式，即函數也有更進一步的一種形式，即tril(A,k)，其功能是求矩陣其功能是求矩陣A的第的第k條對角線以下的元素。條對角線以下的元素。 (2)上三角矩陣上三角矩陣 在在MATLAB中，提

5、取矩陣中，提取矩陣A的上三角矩陣的函數的上三角矩陣的函數是是triu(A)和和triu(A,k)，其用法與提取下三角矩陣，其用法與提取下三角矩陣的函數的函數tril(A)和和tril(A,k)完全相同。完全相同。例：例：A=magic(4),B=tril(A), C=tril(A,1),D=triu(A)二、二、 特殊矩陣的生成特殊矩陣的生成1.魔方矩陣魔方矩陣 函數函數magic(n)，其功能是生成一個，其功能是生成一個n階魔方陣，其階魔方陣，其每一行、每一列以及主對角線元素之和都等于每一行、每一列以及主對角線元素之和都等于sum(1:n2)/n，且矩陣元素在，且矩陣元素在1,2，n2中取值

6、。中取值。 例例1.2 將將101125等等25個數填入一個個數填入一個5行行5列的表格列的表格中，使其每行每列及對角線的和均為中，使其每行每列及對角線的和均為565。 命令如下：命令如下： B=100+magic(5) 2.范得蒙矩陣范得蒙矩陣 函數函數vander(V)生成以向量生成以向量V為基礎向量的范得蒙矩為基礎向量的范得蒙矩陣。陣。例例：V=1 2 3 4； vander(V),rot90(vander(V)3.3.希爾伯特矩陣希爾伯特矩陣 生成希爾伯特矩陣的函數是生成希爾伯特矩陣的函數是hilb(n)。MATLAB中，中，有一個專門求希爾伯特矩陣的逆的函數有一個專門求希爾伯特矩陣的

7、逆的函數invhilb(n)，其功能是求其功能是求n階的希爾伯特矩陣的逆矩陣。階的希爾伯特矩陣的逆矩陣。 例例:求求3階希爾伯特矩陣及其逆矩陣。階希爾伯特矩陣及其逆矩陣。 命令如下：命令如下： format rat % 以有理形式輸出以有理形式輸出 H= hilb(3) invH=invhilb(3)4.4.托普利茲矩陣托普利茲矩陣 生成托普利茲矩陣的函數是生成托普利茲矩陣的函數是toeplitz(x,y)toeplitz(x,y)，它生，它生成一個以成一個以x x為第為第1 1列，列，y y為第為第1 1行的托普利茲矩陣行的托普利茲矩陣, ,除第一除第一行第一列外，其他每個元素都與左上角的元

8、素相同。行第一列外，其他每個元素都與左上角的元素相同。這里這里x, yx, y均為向量，二者不必等長。特別地，當均為向量，二者不必等長。特別地，當x=yx=y時，時，用用toeplitz(x)toeplitz(x)。例例： T=toeplitz(1:3)T=toeplitz(1:3)5.5.帕斯卡矩陣帕斯卡矩陣 函數函數pascal(n)pascal(n)生成一個生成一個n n階的帕斯卡矩陣。階的帕斯卡矩陣。例例1.3 求求(x+y)5的展開式。的展開式。在在MATLAB命令窗口，輸入命令：命令窗口，輸入命令：pascal(6)ans = 1 1 1 1 1 1 1 2 3 4 5 6 1 3

9、 6 10 15 21 1 4 10 20 35 5656 1 5 15 35 70 126 1 6 21 56 126 252此矩陣的次對角線上的元素此矩陣的次對角線上的元素1，5，10，10，5，1即為展開即為展開式的系數。式的系數。三、三、 矩陣結構變換矩陣結構變換1. 矩陣的轉置矩陣的轉置 轉置運算符是單撇號轉置運算符是單撇號()。2. 矩陣的旋轉矩陣的旋轉 矩陣的旋轉利用函數矩陣的旋轉利用函數rot90(A,k)，功能是將矩陣，功能是將矩陣A逆時針旋轉逆時針旋轉90的的k倍，當倍，當k為為1時可省略。時可省略。3. 矩陣的左右翻轉矩陣的左右翻轉 對矩陣對矩陣A實施左右翻轉的函數是實施

10、左右翻轉的函數是fliplr(A)。4. 矩陣的上下翻轉矩陣的上下翻轉 對矩陣對矩陣A實施上下翻轉的函數是實施上下翻轉的函數是flipud(A)。四、四、 矩陣的逆與偽逆矩陣的逆與偽逆1. 矩陣的逆矩陣的逆 求一個矩陣的逆非常容易。求方陣求一個矩陣的逆非常容易。求方陣A的逆可調的逆可調用函數用函數inv(A)。例例1.4 用求逆矩陣的方法解線性方程組。用求逆矩陣的方法解線性方程組。 命令如下：命令如下： A=1,2,3;1,4,9;1,8,27; b=5,2,6; x=inv(A)*b 一般情況下，用左除比求矩陣的逆的方法更一般情況下，用左除比求矩陣的逆的方法更有效，即有效，即x=Ab。2.

11、矩陣的偽逆矩陣的偽逆 MATLAB中，求一個矩陣中，求一個矩陣偽逆偽逆的函數是的函數是pinv(A)。如果矩陣如果矩陣A不是一個方陣，或者不是一個方陣，或者A是一個非滿秩的方陣是一個非滿秩的方陣時，矩陣時，矩陣A沒有逆矩陣，但可以找到一個與沒有逆矩陣，但可以找到一個與A的轉置矩的轉置矩陣陣A同型的矩陣同型的矩陣B，使得：，使得： ABA=A,BAB=B此時稱矩陣此時稱矩陣B為矩陣為矩陣A的偽逆。的偽逆。例例1.5 求求A的偽逆，并將結果賦值給的偽逆，并將結果賦值給B。A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)例例1.6 求矩陣求矩陣A的偽逆。的偽逆。 A=0,0,0

12、;0,1,0;0,0,1;pinv(A)五、五、 方陣的行列式方陣的行列式 求方陣求方陣A所對應的行列式的值的函數是所對應的行列式的值的函數是det(A)。例例1.7 用克萊姆用克萊姆(Cramer)法則求解四元線性方程組。法則求解四元線性方程組。程序如下：程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; %定義系數矩陣定義系數矩陣b=4;6;12;6; %定義常數項向量定義常數項向量D1=b,D(:,2:4); %用方程組的右端向量置換用方程組的右端向量置換D的第的第1列列 D2=D(:,1),b,D(:,3:4); %用方程組的右端向量置用方程組的右端

13、向量置換換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六、六、 矩陣的秩矩陣的秩 MATLAB中，求矩陣秩的函數是中，求矩陣秩的函數是rank(A)。例例：求上例中方程組系數矩陣：求上例中方程組系數矩陣D的秩，命令是：的秩，命令是： D=2,2,-1,1;4,3,

14、-1,2;8,5,-3,4;3,3,-2,2; r=rank(D) r = 4 說明說明D是一個滿秩矩陣。是一個滿秩矩陣。七、七、 向量向量和矩陣和矩陣的范數的范數1. 計算向量計算向量3種常用范數的函數種常用范數的函數 (1)norm(V)或或norm(V,2) 計算向量計算向量V的的2-范數范數 (2)norm(V,1) 計算向量計算向量V的的1-范數范數 (3)norm(V,inf) 計算向量計算向量V的的-范數范數 例例1.8 已知已知V，求，求V的的3種范數。種范數。 命令如下：命令如下： V=-1,1/2,1; v1=norm(V,1) %求求V的的1范數范數 v2=norm(V)

15、 %求求V的的2范數范數 vinf=norm(V,inf) %求求范數范數nniixxxxx 2111),(maxmax211ninixxxxx常用的向量范數常用的向量范數:22221122|nniixxxxx 范數意義下的單位向量范數意義下的單位向量: X=x1, x2T1-11|X|1 = 111-1-1|X|2 = 1-111-1-1|X| = 12. 2. 矩陣的范數及其計算函數矩陣的范數及其計算函數 MATLAB中提供了求中提供了求3 3種矩陣范數的函數，其種矩陣范數的函數，其函數調用格式與求向量的范數的函數完全相同函數調用格式與求向量的范數的函數完全相同例例1.9 求矩陣求矩陣A的

16、三種范數。的三種范數。命令如下：命令如下：A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;11,18,25,2,19; a1=norm(A,1) %求求A的的1-范數范數 a2=norm(A) %求求A的的2-范數范數 ainf=norm(A,inf) %求求A的的-范數范數現求解線性方程組現求解線性方程組 Axb11 11221121 1222221 122mnnnnmmmnna xa xa xba xa xa xba xaxaxbL L L L L L L L L L L L L L L L LLLL情形1：m=n 在MATLAB中

17、的求解命令有：123211113,2;1,-1;-1,1;0.2000-0.8000 xxAbxA bx如：求解設： 則： 得： 情形2：mn（超定方程） ( )* (-1)* / *()xA bxinv AbxAbxb Axbinv A解線性方程組的一般函數文件如下：function x,y=line_solution(A,b) m,n=size(A);y=; if norm(b)0 %非齊次方程組 if rank(A)=rank(A,b) %方程組的相容性條件 if rank(A)=n %有唯一解 x=Ab; else %方程組有無窮多個解,基礎解系 disp(原方程組有有無窮個解，其齊次

18、方程組的基礎 解系為y，特解為x); y=null(A,r); x=Ab; end else %方程組不相容，給出最小二乘解 disp(方程組的最小二乘法解是：); x=Ab; endelse %齊次方程組 if rank(A)=n %列滿秩 x=zeros(n,1) %0解 else %非0解 disp(方程組有無窮個解，基礎解系為x); x=null(A,r); end endreturn %表示終止M函數，返回到調用語句之后。null(A,r)-用來求齊次線性方程組的基礎解系的，用來求齊次線性方程組的基礎解系的，加上加上r則求出的是一組最小正整數解，如果不加，則求出的是一組最小正整數解，

19、如果不加，則求出的是解空間的規范正交基。則求出的是解空間的規范正交基。如在如在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)分別顯示其求解結果。分別顯示其求解結果。直接法：理論，無舍入誤差，有限步，精確解迭代法：格式，用某種極限過程去逐步逼近線性方程組的精確解求解線性方程組的主要方法有：一、一、 GaussGauss消去法消

20、去法設有設有線性代數：方法不好時工作量非常大，線性代數：方法不好時工作量非常大， 工作量小的方法是工作量小的方法是GaussGauss消去法。消去法。 3.1 3.1 解線性方程組的直接法解線性方程組的直接法1 111 2211 ,12 112 2222 ,11122,11122,1nnnnnniii ninnnnn nnn naxaxaxaaxaxaxaaxaxaxaaxaxaxaLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL消消 元：元： 222222(3,)iiaainaia然后用將化為零：把第2行，加到第 行。11 11211,1222222,1

21、22,1 22,1 nnnnnniinni nnnnnn na xa xa xaa xa xaa xa xaa xa xaL L L L L L L L L L L L L L L LL L L L L L L L L L L L L L L LLLLL111111(2, )1iiaainaia用將化為零：把第 行，加到第 行。L00 ?iiiiaa問題：或以此類推，最后方程組化為：回回 代：代： ,1,111(1,2,.,1n nnnnnkk nkjjj kkkaxaxaa xaknn 11 112211,122222,1 ,1 nnnnnnnnnn na xa xa xaa xa xaa

22、xaL L L L L L L L L L LLL( (上三角方程組上三角方程組) ) 二、列主元素消去法二、列主元素消去法-計算結果可靠計算結果可靠1111111111(1)max112(j1,2,n1)3jjrii njjrjjrraaracaaca 找行號使，對調行：11111111110 :1,2,3,1,2,1iiiijjijaaaiijaaaaainjna 消元：用消為第 行加到第 行 第 行第 個元素成為（；）到此原方程組化為到此原方程組化為11 11211,1222222,1 22,1 22,1 nnnnnniinni nnnnnn na xa xa xaa xa xaa xa

23、 xaa xa xaL L L L L L L L L L L L L L L LL L L L L L L L L L L L L L L LLLLL到此原方程組化為到此原方程組化為222222(2)m ax2.irinraar 找， 使， 對 調行22222222220 (3,4, ) :2 3,4,2,3,1iiiijjijaainaiaaaaainjna 消元：用把消為第 行第 行，則（，；，）11 11221331,1122223322,1 333,133 3 nnnnnnnnnna xa xa xa xaa xa xa xaa xa xaaL L L L L L L L L L L

24、 L L L LLLL,13nnnn nxa xaL( (上三角方程組上三角方程組) ) （3.2)3.2)( (n-1) ) 原方程組化為原方程組化為以上為消元過程。以上為消元過程。L L L L L11 112 211, 122 222, 1 , 1 n nnn nnnn nnna xa xa xaa xa xaa xaL L L L L L L L L L LLL( (n) ) 回代求解公式回代求解公式 ,1,111(1,2,.,1n nnnnnkk nkjjj kkkaxaxaa xaknn （3.3) 3.3) 是回代過程。是回代過程。 （3.3)3.3)說明說明：(1)也可采用無回

25、代的列主元消去法也可采用無回代的列主元消去法(叫叫Gauss-Jordan消去法消去法)，該法同時消去對角線上方和下方，該法同時消去對角線上方和下方的元素，且仍舊需要選主元，但比有回代的列主的元素，且仍舊需要選主元，但比有回代的列主元消去法的乘除運算次數多。元消去法的乘除運算次數多。Gauss-Jordan消去消去法的優點之一是比較合適用它來求一個矩陣的逆法的優點之一是比較合適用它來求一個矩陣的逆矩陣。矩陣。(2)有回代的列主元消去法所進行的乘除運算次數有回代的列主元消去法所進行的乘除運算次數為為 ，計算量很小。，計算量很小。321133nnn例例1 在在MATLAB上，用列主元素消去法求解方

26、程組：上，用列主元素消去法求解方程組：1230.040.040.1230.561.560.3210.241.240.280 xxx 程序如下：程序如下： cleara = -0.04 0.04 0.12 3; 0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0 x = 0,0,0; tempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo;aa(2,:) = a(2,:) - a(1,:)*a(2,1)/a(1,1);a(3,:) = a(3,:) - a(1,:)*a(3,1)/a(1,1); a %看看a的結果再決定下一步的結果再決

27、定下一步tempo = a(3,:); a(3,:) = a(2,:); a(2,:)=tempo;aa(3,:) = a(3,:) - a(2,:)*a(3,2)/a(2,2); a %化為上三角化為上三角x(3) = a(3,4)/a(3,3); %開始回代開始回代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運行得方程組的解為：7725x,1,111(1,2,.,1n nnnnnkk nkjjj kkkaxaxaa xaknn 例例2 用用Gauss-Jordan消去法求解上例中

28、的矩陣消去法求解上例中的矩陣 的逆矩陣。的逆矩陣。clearA = -0.04 0.04 0.12 ; 0.56 -1.56 0.32 ; -0.24 1.24 -0.28a=A,eye(3);atempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo;aa(1,:) = a(1,:)/a(1,1)for i=2:3; a(i,:)=a(i,:) - a(i,1)*a(1,:); end;atempo = a(3,:); a(3,:) = a(2,:); a(2,:)=tempo;aa(2,:)=a(2,:)/a(2,2);afor i=1:3 if i=2,

29、 a(i,:)=a(i,:)-a(i,2)*a(2,:); endendaa(3,:)=a(3,:)/a(3,3)for i=1:3; if i=3, a(i,:)=a(i,:)-a(i,3)*a(3,:); end;end;aA_inv = a(:,4:6)A*A_inv145_212811A inv結果為：三、三、 GaussGauss全主元消去法：全主元消去法： 優點優點-計算結果更可靠；計算結果更可靠； 缺點缺點-挑主元花時間更多，挑主元花時間更多， 次序有變動，程序復雜。次序有變動，程序復雜。1,nxxL 四、應用四、應用 （1）求行列式）求行列式 （2）求逆矩陣）求逆矩陣 ( (以

30、上過程都應選主元以上過程都應選主元) ) 111111111.( 1).( 1).nnmmnnnnnnnaabbbbaab 1().()A EE A在MATLAB中用命令det(A)在MATLAB中用命令inv(A)或A(-1)或rref(A,E) 記記，則，則 ( (下三角下三角上三角）（三角因子分解）上三角）（三角因子分解） Gauss Gauss消元，初等行變換，化原方程組為上三角型。消元，初等行變換，化原方程組為上三角型。五、矩陣三角分解法五、矩陣三角分解法1211211121|1|kkkkkkA bU gL LL LA bU gL LL L AUAL LL LULA=8 1 6.00

31、; 3 5 7.00;4 9 2.00;9 9 0; B=spconvert(A)B = (8,1) 6 (3,5) 7 (4,9) 2clear A=eye(3) B=sparse(A)B = (1,1) 1 (2,2) 1 (3,3) 1sparse(2,3)ans = All zero sparse: 2-by-3i,j,s = find(A),A = sparse(i,j,s)i= j= s= A=1 1 1 (1,1) 12 2 1 (2,2) 13 3 1 (3,3) 1full(A)3. 單位稀疏矩陣的產生單位稀疏矩陣的產生 單位矩陣只有對角線元素為單位矩陣只有對角線元素為1，其

32、他元素都為，其他元素都為 0，是一，是一種具有稀疏特征的矩陣。我們知道，函數種具有稀疏特征的矩陣。我們知道，函數 eye產生一個完全產生一個完全存儲方式的單位矩陣。存儲方式的單位矩陣。 MATLAB還有一個產生稀疏存儲方式的單位矩陣的函還有一個產生稀疏存儲方式的單位矩陣的函數，這就是數，這就是speye。函數。函數speye(m,n)返回一個返回一個mn的稀疏存儲的稀疏存儲單位矩陣。單位矩陣。例：例：C=eye(2,3),D=speye(2,3)C = 1 0 0 0 1 0D = (1,1) 1 (2,2) 1Axb一一. .簡單迭代法簡單迭代法 1.1.迭代法建立迭代法建立. . 考慮考慮

33、AxbxBxg ( (矩陣矩陣B不唯一不唯一) )其迭代格式其迭代格式 (1)( )(0) (0,1,2,) (3.4)kkxBxgkx取定初始向量產生向量序列產生向量序列(1)(2)()(1),kkxxxx若收斂若收斂, ,記記(1)limkkxx 求解線性方程組的迭代法求解線性方程組的迭代法則于則于(3.4)(3.4)兩端取極限有：兩端取極限有： ,xBxg上式說明：上式說明： 是解向量是解向量 ，從而當，從而當k充分大時充分大時(1)kx注意注意: : 迭代陣迭代陣B不唯一，不唯一，B的選取的選取影響收斂性。影響收斂性。 解向量解向量(3.4)(3.4)叫簡單迭代法叫簡單迭代法, ,B叫

34、迭代矩陣。叫迭代矩陣。 2. 2.收斂性收斂性. . 定義定義3.3 3.3 稱稱1( )max |( ) |ii nBB 為矩陣為矩陣B B的譜半徑。的譜半徑。 xxx(0):1Bx注意（ ）時不能說“對任意都不收斂”。(1)( )11111(0) (0,1,2,),|max|1 |max|1 kknnijijj ni nijxBxgkBbBbx 對迭代法若或則其對任意收斂。 定理定理3.4 定理定理3.3 對于簡單迭代法，若迭代式對于簡單迭代法，若迭代式(1)( )(0) 1 kkxBxgxB對任意初始向量都收斂（ ）設有方程組設有方程組( ( 其中其中 ) ) Ax = b，即，即 0i

35、ia 11112211211222221122nnnnnnnnnna x a xa xba x a xa xba x a xa xbLLL L L L L L L L L L L L LL(3.5)(3.5) 作等價變形作等價變形11112211122211222211221()1()1()00 0nnnnnnnnnnnxbx a xa xaxba xxa xaxba x a xxa (3.6)(3.6)11112211211222221122nnnnnnnnnna xa xa xba xa xa xba xa xa xbLLL L L L L L L L L L L L LL(3.5)(3.

36、5)-Jacobi-Jacobi迭代法迭代法22(1)( )( )( )111122111(1)( )( )( )2221 122(1)( )( )( )1 1221(0.)1(0.)1(.0)kkkknnkkkknnkkkknnnnnnnxbxa xa xaxba xxa xaxba xa xxaL L L L L L L L L L L L L于是有迭代公式于是有迭代公式(k=0,1,2,)(3.7)(3.7)(1)( )1nkkiiijjiijj ixba xa矩陣形式為：矩陣形式為：111211( )1111(1)112( )2(1)2122222222( )(1)21200 0nkk

37、knkkknnnnnnnnnnaabaaaxxbaaxxaaaxxbaaaaaLLMMMLLLLL(1)( )kkJxB xg簡記為簡記為對于矩陣方程Ax=b的雅可比迭代格式還可以通過下面的方法來得到：分裂分裂: Ax = b任取任取 x(0), 迭代計算產生向量序列迭代計算產生向量序列:若若*lim)(xxkk則則x* 是方程組是方程組 Ax = b 的解的解迭代矩陣迭代矩陣迭代計算迭代計算格式格式( k=0,1,2, )x(k+1)= BJx(k) + fx(1),x(2),x(k),splittingD + ( A D )x = BJx + fJ BJ=D-1(D-A) ,fJ =D-1

38、bDx = (D A )x + b =D是是A的主對角線上的元素組成的對角矩陣的主對角線上的元素組成的對角矩陣其中，其中，Jacobi迭代矩陣迭代矩陣 0002122111212211nnnnnnJaaaaaaaaaB 0/0/02122222211111112nnnnnnnnJaaaaaaaaaaaaB nnnJabababf/222111(0)1(0):(1)(3.7)()1(2)|1|1,(3.7).JJJJacobiJacobixBBBJacobix迭代法收斂性迭代法對任意收斂若或則迭代法對任意收斂 (3)設方程組（3.5）的系數矩陣A按行嚴格對角占優 即： 1, (i1,2,n)ni

39、iijjj iaa1, (j1,2,n)njjijiijaa 或按列嚴格對角占優，即(0) (3.7) .Jacobix則迭代法對任意收斂(1)( )1nkkiiijjiijj ixba xa例例123123123971081513xxxxxxxxx(1)( )( )312(13)/15kkkxxx(1)( )11(1)( )22(1)( )3301/91/97/91/1001/108/101/151/15013/15kkkkkkxxxxxx(1)( )( )123(7)/9kkkxxx(1)( )( )213(8)/10kkkxxx123213312(7)/9(8)/10(13)/15xxx

40、xxxxxx(0)1(0)2(0)3000 xxx 等價變形Jacobi迭代算法迭代算法A=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;err=1;k=0;while err0.00005 err=0;k=k+1; for i=1:3 s=0;t=x(i); for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=t; y(i)=(b(i)-s)/A(i,i); err=max(abs(x(i)-y(i),err); end x=y;xend123123123971081513xxxxxxxxx0.7778 0.80

41、00 0.86670.9630 0.9644 0.97190.9929 0.9935 0.99520.9987 0.9988 0.99910.9998 0.9998 0.99981.0000 1.0000 1.00001.0000 1.0000 1.0000(1)( )1nkkiiijjiijj ixba xa二、二、 迭代法迭代法 設有簡單迭代法設有簡單迭代法 即即Seidel(1)( )kkxBxg(1)( )( )( )111111122(1)( )( )( )221122222(1)( )( )( )1122kkkknnkkkknnkkkknnnnnnngggxb xb xb xxb

42、xb xb xxb xb xb xLLL L L L L L L L L L L L L L L LL(3.8)(3.8) 稱如下迭代法稱如下迭代法 (3.9) (3.9)(1)( )( )( )111122111(1)(1)( )( )22222122(1)1(1)(1)( )1212kkkkkkknnkkknnkknnnnnnngggxb xb xb xxbb xb xxxbbxbxxLLL L LL 為與為與(3.8)(3.8)對應的對應的Seidel迭代法，其迭代矩陣迭代法，其迭代矩陣 可用可用 “代入法代入法”求得。求得。sB(1)( )( )( )111122111(1)(1)(

43、)( )2221 12222(1)(1)(1)( )1 1221 01 0 1 0kkkknnkkkknnkkkknnnnnnnxbxa xa xaxba xxa xaxba xa xxa于是有下面的迭代公式：- Seidel迭代法1(1)(1)( )11inkkkiiijjijjiijj ixba xa xa (1)( )1nkkiiijjiijj ixba xa- Jacobi迭代法將方程組將方程組Ax = b 的系數矩陣的系數矩陣 A 分解分解 A = D U L nnaaaD22110001,121nnnaaaL1211,000nnnaaUaGauss-Seidel迭代法的矩陣表示迭代

44、法的矩陣表示 nijkjijijkjijikiiixaxabxa1)(11)1()1( nijkjijiijkjijxabxa1)(1)1(i = 1,2,n)()(2)(1, 111221)1()1(2)1(121222111000knkknnnnknkknnnnxxxaaabbbxxxaaaaaa (D L)x(k+1) = b +Ux(k) x(k+1) = (D L)-1b + (D L)-1Ux(k) 記記 BG-S=(D L)-1U, fG-S=(D L)-1b Gauss-Seidel迭代格式迭代格式: x(k+1)=BG-Sx(k)+fG-S 000, 111212122211

45、1nnnnnnnSGaaaaaaaaaB nnnnnSGbbbaaaaaaf21121222111 （1 1） 迭代法迭代法(3.9)(3.9)對任意對任意 收斂收斂（2 2）若）若 則則 迭代法迭代法(3.9)(3.9)對對 任意任意 收斂；收斂；（3 3）若簡單迭代法）若簡單迭代法(3.8)(3.8)的迭代矩陣的迭代矩陣 滿滿 足足 或或 ，則相應的，則相應的SeidelSeidel迭代法迭代法 (3.9)(3.9)對任意對任意 收斂收斂. . SeidelSeidel(0)x; 1)( Bs Seidel(0)x()n nijBb11B1B(0)x迭代法迭代法(3.9)(3.9)的收斂性

46、的收斂性1|1|1JJBB或()1SB1(1)(1)( )11inkkkiiijjijjiijj ixba xa xa 例例123123123971081513xxxxxxxxx(1)(1)(1)312(13)/15kkkxxx(1)( )11(1)( )22(1)( )331000 1/91/97/91/1010001/108/101/151/15 100013/15kkkkkkxxxxxx(1)( )( )123(7)/9kkkxxx(1)(1)( )213(8)/10kkkxxx123213312(7)/9(8)/10(13)/15xxxxxxxxx(0)1(0)2(0)3000 xxx

47、 Gauss-Seidel迭代算法迭代算法123123123971081513xxxxxxxxxA=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;err=1;k=0;while err0.00005 err=0;k=k+1; for i=1:3 s=0;t=x(i); for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=(b(i)-s)/A(i,i); err=max(abs(x(i)-t),err); end xend 0.7778 0.8778 0.9770 0.9839 0.9961 0.9987 0.999

48、4 0.9998 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001(1)(1)( )11inkkkiiijjijjiijj ixba xa xa Gauss-seidel迭代迭代A=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;err=1;k=0;while err0.00005 err=0;k=k+1; for i=1:3 s=0;t=x(i); for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=(b(i)-s)/A(i,i); err=max(abs(x(i)

49、-t),err); end xendJacobi迭代迭代A=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;err=1;k=0;while err0.00005 err=0;k=k+1; for i=1:3 s=0;t=x(i); for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=t; y(i)=(b(i)-s)/A(i,i); err=max(abs(x(i)-y(i),err); end x=y;xend (3.11)是一種加權平均。是一種加權平均。 - -松弛因子，松弛因子， 即即Seidel方法方法=1三三

50、. . 超松弛迭代法超松弛迭代法( (SOR法法) )(1, 2,) in1(1)(1)( )11inkkkiiijjijjiijj ixba xa xa Seidel迭代法迭代法1(1)( )(1)( )11(1)inkkkkiiiijjijjiijj ixxba xa xa SOR法法(3.11)特點：特點： SOR法是Gauss-Seidel方法的一種加速方法，是解大型稀疏矩陣方程組的有效方法之一，具有計算公式簡單，程序設計容易，占用計算機內存少等優點，但需要選擇好的加速因子即最佳松弛因子。SOR方法的收斂性如下（不加證明）：方法的收斂性如下（不加證明）：(1)(1)SOR方法對任意方法

51、對任意 都收斂的必要條件是：都收斂的必要條件是： (2)(2)若系數矩陣若系數矩陣A對稱正定，則對稱正定，則 時時SOR方法方法 求解求解 對任意對任意 收斂；收斂；(3)(3)若系數矩陣若系數矩陣A按行（或按列）嚴格對角占優，按行（或按列）嚴格對角占優， 則則 時時SOR方法對任意方法對任意 收斂。收斂。02Axb(0)x0102(0)x(0)x如前例，用如前例，用SOR法求解法求解123123123971081513xxxxxxxxxA=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;err=1;k=0;w=1.1;while err0.00005 e

52、rr=0;k=k+1; for i=1:3 s=0;t=x(i); for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=w*(b(i)-s)/A(i,i)+(1-w)*x(i); err=max(abs(x(i)-t),err); end xend0.8556 0.9741 1.08751.0220 1.0146 0.99390.9988 0.9977 1.00040.9999 1.0003 1.00001.0000 1.0000 1.00001.0000 1.0000 1.00001( 1)( )( 1)( )111(1)inkkkkiiiij jij

53、 jjj iiixxbaxaxa 3.3 3.3 不可解問題不可解問題 線性方程組并不總是數值可解的，考慮如下三線性方程組并不總是數值可解的，考慮如下三個方程組。個方程組。121212121212121( )2221( )022( )120 xxAxxxxBxxxxCxxxx 在MATLAB中分別求解如下：( ) 1,1; 2,21;2warning:Matrix is singular to working precision. ans= NaN NaN 1,11ans10( ) 1,1; 1,11;0warning:Matrix is singular to working precisi

54、on. ans=-Inf -IAB 若僅求解第一個方程：nf( )1,2; 1,1;2, 1 2;1;0ansC -0.6-0.6有無窮多個解其中一個解方程組不相容最小二乘解 3.4 3.4 病態問題病態問題 有許多線性方程組理論上是可解的，但實際計算有許多線性方程組理論上是可解的，但實際計算中由于受到舍入誤差的干擾而無法得到精確解，此類中由于受到舍入誤差的干擾而無法得到精確解，此類問題稱為問題稱為病態問題病態問題。如何來判斷它是否是病態問題呢？。如何來判斷它是否是病態問題呢？為此引入矩陣條件數的概念。先給出矩陣范數的定義為此引入矩陣條件數的概念。先給出矩陣范數的定義和性質：和性質：常用矩陣范

55、數：1111/2maxmax211maxmaxnijj niTniji njAaAA AAa (B)表示矩陣B的最大特征值范數性質：范數性質：矩陣條件數的概念：矩陣條件數的概念： 方程組方程組 Ax = b, 右端項右端項 b 有一擾動有一擾動 ， 引起方程組解引起方程組解 x 的擾動的擾動 。bx設設 x 是方程組是方程組 Ax = b 的解的解,則有則有bbxxA )(化簡化簡,得得bxAbAx1|1bAx 由由 Ax = b 得得 |xAb |1bAx 所以所以|)|(|1bbAAxx 定義條件數定義條件數: Cond(A) = |A1 | |A| 或或 C(A) = |A1 | |A|

56、當條件數很大時當條件數很大時,方程組方程組 Ax = b是病態問題是病態問題;當條件當條件數較小時數較小時,方程組方程組 Ax = b是良態問題。是良態問題。|)(|1|)(|11111 AAIAAAAI IAAIAAI 111)( 注注:|11|)(|111AAAAI 11111)()( AAIAAIAAI （下面會用到）類似類似,設方程組設方程組 Ax = b,矩陣矩陣A 有一擾動有一擾動 時時,將引起方程組解將引起方程組解x的擾動的擾動 。Ax設設 x 是方程組是方程組 Ax = b 的解，則有的解，則有bxxAA )( 化簡化簡,得得AxxAA )(AxAAAIx 111)( 取范數取

57、范數|)(|111xAAAAIx |1|11AAAAAAxx 矩陣矩陣A的條件數記為的條件數記為Cond(A)，定義為：，定義為：1( )Cond AAA條件數總滿足：條件數總滿足：( )1Cond A 注：當矩陣是病態時，其條件數一定很大。注：當矩陣是病態時，其條件數一定很大。MATLAB中計算條件數的命令是：中計算條件數的命令是： cond(A) 對于對于病態矩陣病態矩陣，逆矩陣和行列式的計算都會變，逆矩陣和行列式的計算都會變得不精確。所以具備下列特征的問題可認為是病態得不精確。所以具備下列特征的問題可認為是病態的：的：11det( )*det()12( );(3)*( )4( *( )1*( )AAinv inv AAA inv Acond A inv AA inv A（）的計算值與 有較大偏離；（ ）的計算值明顯不等于的計算值與單位矩陣有較大偏離；（ ）與 的偏離可作為與單位陣的偏差的精確度量。條件數的性質：條件數的性質：a) 1)(Acond； b) 對于R )0(，)()(AcondAcond； c) 對于正交陣nnRQ， )()()(AcondAQcondQAcond； 例：例：Hilbert矩陣（非常有名的病態矩陣）：