矩陣分解與線性方程組求解市公開課獲獎?wù)n件

1、MATLAB矩陣分解與線性方程組求解數(shù)學(xué)與信息科學(xué)系 汪遠(yuǎn)征第1頁第1頁4. MATLAB矩陣分解與線性方程組求解4.1 矩陣分解4.2 秩與線性相關(guān)性4.3 線性方程組求解4.4 特性值與二次型第2頁第2頁4.1 矩陣分解4.1.1 LU分解矩陣三角分解又稱LU分解，它目的是將一個矩陣分解成一個下三角矩陣L和一個上三角矩陣U乘積，即A=LU。Matlab使用函數(shù)lu實(shí)現(xiàn)LU分解，其格式為：L,U = lu(A)其中U為上三角陣，L為下三角陣或其變換形式，滿足LU=A。第3頁第3頁4.1 矩陣分解4.1.1 LU分解L,U,P = lu(A)U為上三角陣，L為下三角陣，P為單位矩陣行變換矩陣，

2、滿足LU=PA。例4-1A=1 2 3;4 5 6;7 8 9;L,U=lu(A)L,U,P=lu(A)第4頁第4頁4.1 矩陣分解4.1.2 Cholesky分解假如A為n階對稱正定矩陣，則存在一個實(shí)非奇異上三角陣R，滿足R*R = A，稱為Cholesky分解Matlab使用函數(shù)chol實(shí)現(xiàn)Cholesky分解，其格式為：R = chol(A)若A非正定，則產(chǎn)生錯誤信息。R,p = chol(A)不產(chǎn)生任何錯誤信息，若A為正定陣，則p=0，R與上相同；若A非正定，則p為正整數(shù)，R是有序上三角陣。第5頁第5頁4.1 矩陣分解4.1.2 Cholesky分解例4-2A=pascal(4) %產(chǎn)

3、生4階pascal矩陣R,p=chol(A)第6頁第6頁4.1 矩陣分解4.1.3 QR分解將矩陣A分解成一個正交矩陣Q與一個上三角矩陣R乘積A=QR，稱為QR分解。Matlab使用函數(shù)qr實(shí)現(xiàn)QR分解，其格式為：Q,R = qr(A)Q,R,E = qr(A)求得正交矩陣Q和上三角陣R，E為單位矩陣變換形式，R對角線元素按大小降序排列，滿足AE=QR。Q,R = qr(A,0)產(chǎn)生矩陣A“經(jīng)濟(jì)大小”分解第7頁第7頁4.1 矩陣分解4.1.3 QR分解Matlab使用函數(shù)qr實(shí)現(xiàn)QR分解，其格式為：Q,R = qr(A)Q,R,E = qr(A)求得正交矩陣Q和上三角陣R，E為單位矩陣變換形式

4、，R對角線元素按大小降序排列，滿足AE=QR。Q,R = qr(A,0)產(chǎn)生矩陣A“經(jīng)濟(jì)大小”分解Q,R,E = qr(A,0)E作用是使得R對角線元素降序, 且Q*R=A(:, E)第8頁第8頁4.1 矩陣分解4.1.3 QR分解例4-3 A = 1 2 3;4 5 6; 7 8 9; 10 11 12;Q,R = qr(A)第9頁第9頁4.1 矩陣分解4.1.3 QR分解返回將矩陣A第j列移去后新矩陣qr分解使用函數(shù) qrdelete，其格式為：Q,R = qrdelete(Q,R,j)例4-4A=-149 -50 -154;537 180 546;-27 -9 -25;Q,R=qr(A)

5、Q,R=qrdelete(Q,R,3) %將A第3列去掉后進(jìn)行qr分解第10頁第10頁4.1 矩陣分解4.1.3 QR分解在矩陣A中第j列插入向量x后新矩陣進(jìn)行qr分解使用函數(shù) qrinsert，其格式為：Q,R = qrinsert(Q,R,j,x)若j不小于A列數(shù)，表示在A最后插入列x。例4-5A=-149 -50 -154;537 180 546;-27 -9 -25;x=35 10 7;Q,R=qr(A)Q,R=qrinsert(Q,R,4,x)第11頁第11頁4.1 矩陣分解4.1.4 Schur分解矩陣Schur分解：設(shè)A Rnn ，則存在正交陣Q使 實(shí)Schur型其中Rii至多2

6、階。若1階，其元素即A特性值，若2階其特性值為A一對共軛復(fù)特性值。Matlab使用函數(shù)schur實(shí)現(xiàn)Schur分解，其格式為：R = schur(A)R為schur矩陣, 即R主對角線元素為特性值三角陣第12頁第12頁4.1 矩陣分解4.1.4 Schur分解Matlab使用函數(shù)schur實(shí)現(xiàn)Schur分解，其格式為：R = schur(A)R為schur矩陣, 即R主對角線元素為特性值三角陣R = schur(A,flag)若A有復(fù)特性根，則flag=complex，不然flag=real。Q,R = schur(A,)返回正交矩陣Q和schur矩陣R，滿足A = Q*R*Q。第13頁第13

7、頁4.1 矩陣分解4.1.4 Schur分解例4-6H = -149 -50 -154; 537 180 546; -27 -9 -25 ;Q, R=schur(H)第14頁第14頁4.1 矩陣分解4.1.5 實(shí)Schur分解轉(zhuǎn)化成復(fù)Schur分解將實(shí)舒爾形式轉(zhuǎn)化成復(fù)舒爾形式函數(shù)是rsf2csf，其格式Q, R = rsf2csf(q, r)例4-7 A=1 1 1 3;1 2 1 1;1 1 3 1;-2 1 1 4;q, r=schur (A)Q, R=rsf2csf(q, r)第15頁第15頁4.2 秩與線性相關(guān)性4.2.1 矩陣和向量組秩與向量組線性相關(guān)性矩陣A秩是指矩陣A中最高階非零

8、子式階數(shù)，或是矩陣線性無關(guān)行數(shù)與列數(shù)；向量組秩通常由該向量組構(gòu)成矩陣來計算。在MATLAB中，求矩陣秩函數(shù)是rank。其格式為：k = rank(A)返回矩陣A行(或列)向量中線性無關(guān)個數(shù)k = rank(A,tol)tol為給定誤差第16頁第16頁4.2 秩與線性相關(guān)性4.2.1 矩陣和向量組秩與向量組線性相關(guān)性例4-8 求向量組1 = (1 2 2 3)，2 = (2 4 1 3)，3 = (1 2 0 3)，4 = (0 6 2 3)，5 = (2 6 3 4)秩，并判斷其線性相關(guān)性。解：A=1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4;k=r

9、ank(A)結(jié)果為由于秩為3 向量個數(shù)，因此向量組線性相關(guān)。第17頁第17頁4.2 秩與線性相關(guān)性4.2.2 求行階梯矩陣及向量組基行階梯使用初等行變換，矩陣初等行變換有三條：(1) 互換兩行ri rj (第i、第j兩行互換)(2) 第i行k倍kri(3) 第i行k倍加到第j行上去rj+kri通過這三條變換能夠?qū)⒕仃嚮尚凶詈喰?，從而找出列向量組一個最大無關(guān)組。第18頁第18頁4.2 秩與線性相關(guān)性4.2.2 求行階梯矩陣及向量組基Matlab將矩陣化成行最簡形命令是rref或rrefmovie。其格式為：R = rref(A)R 是A行最簡行矩陣R,jb = rref(A)jb是一個向量，其

10、含義為：r = length(jb)為A秩；A(:, jb)為A列向量基；jb中元素表示基向量所在列。第19頁第19頁4.2 秩與線性相關(guān)性4.2.2 求行階梯矩陣及向量組基Matlab將矩陣化成行最簡形命令是rref或rrefmovie。其格式為：R,jb = rref(A,tol)tol為指定精度rrefmovie(A)給出每一步化簡過程第20頁第20頁4.2 秩與線性相關(guān)性4.2.2 求行階梯矩陣及向量組基例4-9 求向量組a1=(1,-2,2,3)，a2=(-2,4,-1,3)，a3=(-1,2,0,3)，a4=(0,6,2,3)，a5=(2,-6,3,4)一個最大無關(guān)組a1=1 -2

11、 2 3; a2=-2 4 -1 3; a3=-1 2 0 3;a4=0 6 2 3; a5=2 -6 3 4;A=a1 a2 a3 a4 a5R,jb=rref(A)A(:,jb)即：a1 a2 a4為向量組一個基。第21頁第21頁4.3 線性方程組求解我們將線性方程求解分為兩類：一類是方程組求唯一解或求特解，另一類是方程組求無窮解即通解。能夠通過系數(shù)矩陣秩來判斷：若系數(shù)矩陣秩r=n(n為方程組中未知變量個數(shù))，則有唯一解；若系數(shù)矩陣秩rn，則也許有無窮解；線性方程組無窮解 = 相應(yīng)齊次方程組通解+非齊次方程組一個特解；其特解求法屬于解第一類問題，通解部分屬第二類問題。第22頁第22頁4.3

12、 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解當(dāng)系數(shù)矩陣A為N*N方陣時，MATLAB會自行用高斯消去法求解線性代數(shù)方程組。若右端項(xiàng)b為N*1列向量，則x=Ab可取得方程組數(shù)值解x（N*1列向量）；若右端項(xiàng)b為N*M矩陣，則x=Ab可同時取得同一系數(shù)矩陣A、M個方程組數(shù)值解x（為N*M矩陣），即x(:,j)=Ab(:,j)，j=1,2,M。第23頁第23頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解解法：X=Ab例4-10 求方程組 解。解：A=5 6 0 0 0 1 5 6 0 0 0 1 5 6

13、 0 0 0 1 5 6 0 0 0 1 5; B=1 0 0 0 1;R_A=rank(A) %求秩X=AB %求解第24頁第24頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解解：A=5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5; B=1 0 0 0 1;R_A=rank(A) %求秩X=AB %求解運(yùn)營后結(jié)果下列這就是方程組解。第25頁第25頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解或用函數(shù)rref求解：C=A,B%由系數(shù)矩陣和常

14、數(shù)列構(gòu)成增廣矩陣CR=rref(C)%將C化成行最簡行則R最后一列元素就是所求之解。第26頁第26頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解例4-11 求方程組 一個特解解：A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;B=1 4 0;X=AB%由于系數(shù)矩陣不滿秩, 該解法也許存在誤差X = 0 0 -0.5333 0.6000 一個特解近似值第27頁第27頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解例4-11 求方程組 一個特解解：若用rref求解，則比較準(zhǔn)確：A=

15、1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;B=1 4 0;C=A,B; %構(gòu)成增廣矩陣R=rref(C)由此得解向量X=1.2500 0.2500 0 0(一個特解)。第28頁第28頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解例4-12 解方程組 (1)Ax=b1；(2)Ay=b2 解法1：分別解方程組 (1)Ax=b1；(2)Ay=b2A=1 -1 1;5 -4 3;2 1 1;b1=2;-3;1;b2=3;4;-5;x=Ab1y=Ab2第29頁第29頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣

16、除法求線性方程組AX=b特解例4-12 解方程組 (1)Ax=b1；(2)Ay=b2 解法1：分別解方程組 (1)Ax=b1；(2)Ay=b2得兩個線性代數(shù)方程組解：(1) x1= -3.8， x2= 1.4， x3= 7.2； (2) y1= -3.6， y2= 2.2， y3= 4.4第30頁第30頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解1利用矩陣除法求線性方程組AX=b特解例4-12 解方程組 (1)Ax=b1；(2)Ay=b2 解法2：將兩個方程組連在一起求解：Az = bb=2 3;-3 4;1 -5z=Ab很明顯，這里解z兩個列向量便是前面分別求得兩組解x和y第

17、31頁第31頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(1)LU分解：LU分解又稱Gauss消去分解，可把任意方陣分解為下三角矩陣和上三角矩陣乘積。即A=LU，L為下三角陣，U為上三角陣。則：A*X=b變成L*U*X=b因此X=U(Lb), 這樣能夠大大提升運(yùn)算速度。第32頁第32頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(1)LU分解：例4-13 求方程組 一個特解。解：A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;D=det(A)L,U=lu

18、(A)X=U(LB)第33頁第33頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(1)LU分解：例4-13 求方程組 一個特解。解：闡明 結(jié)果中警告是由于系數(shù)行列式為零產(chǎn)生。能夠通過A*X驗(yàn)證其正確性。第34頁第34頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(2) Cholesky分解若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置乘積，即：A = R*R其中R為上三角陣。方程A*X=b變成R*R*X = b因此X=R(Rb) 第35頁第3

19、5頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(3)QR分解對于任何長方矩陣A，都能夠進(jìn)行QR分解，其中Q為正交矩陣，R為上三角矩陣初等變換形式，即：A = QR方程A*X=b變形成QRX=b因此X=R(Qb)第36頁第36頁4.3 線性方程組求解4.3.1 求線性方程組唯一解或特解2利用LU、QR和cholesky分解求方程組解(3)QR分解上例中Q, R=qr(A)X=R(QB)闡明：這三種分解，在求解大型方程組時很有用。其長處是運(yùn)算速度快、能夠節(jié)約磁盤空間、節(jié)約內(nèi)存。第37頁第37頁4.3 線性方程組求解4.3.2 求線性齊

20、次方程組通解在Matlab中，函數(shù)null用來求解零空間，即滿足AX=0解空間，事實(shí)上是求出解空間一組基(基礎(chǔ)解系)。其格式為：z = null(A)z列向量為方程組正交規(guī)范基，滿足z z = I。z = null(A,r)z列向量是方程AX=0有理基第38頁第38頁4.3 線性方程組求解4.3.2 求線性齊次方程組通解例4-14 求方程組 通解：解：A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;format rat %指定有理式格式輸出B=null(A,r) %求解空間有理基運(yùn)營后顯示結(jié)果下列：第39頁第39頁4.3 線性方程組求解4.3.2 求線性齊次方程組通解例4-14

21、求方程組 通解：解：或通過行最簡型得到基：B=rref(A)即可寫出其基礎(chǔ)解系(與上面結(jié)果一致)。第40頁第40頁4.3 線性方程組求解4.3.2 求線性齊次方程組通解例4-14 求方程組 通解：解：寫出通解：syms k1 k2X=k1*B(:,1)+k2*B(:,2)%寫出方程組通解pretty(X) %讓通解表示式愈加精美運(yùn)營后結(jié)果下列：下面是其簡化形式第41頁第41頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解非齊次線性方程組需要先判斷方程組是否有解，若有解，再去求通解。因此，環(huán)節(jié)為：第1步: 判斷AX=b是否有解，若有解則進(jìn)行第二步第2步: 求AX=b一個特解第3步: 求

22、AX=0通解第4步: AX=b通解= AX=0通解+AX=b一個特解第42頁第42頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-15 求解方程組解：在Matlab中建立M文獻(xiàn)下列：A=1 -2 3 -1;3 -1 5 -3;2 1 2 -2; b=1 2 3; B=A b;n=4; R_A=rank(A), R_B=rank(B), format ratif R_A=R_B&R_A=n%判斷有唯一解 X=Abelseif R_A=R_B&R_An%判斷有無窮解 X=Ab %求特解 C=null(A,r) %求AX=0基礎(chǔ)解系else X=equition no solve %

23、判斷無解end第43頁第43頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-15 求解方程組解：在Matlab中建立函數(shù)文獻(xiàn)下列：A=1 -2 3 -1;3 -1 5 -3;2 1 2 -2; b=1 2 3;jie(A,b)第44頁第44頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-15 求解方程組解：運(yùn)營后結(jié)果顯示：闡明 該方程組無解第45頁第45頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-16 求方程組 通解解法一：調(diào)用jie函數(shù)A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; b=1 4 0; jie(A,b)第46

24、頁第46頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-16 求方程組 通解解法一：在Matlab編輯器中建立M文獻(xiàn)下列：運(yùn)營后結(jié)果顯示為：因此原方程組通解為X=k1 +k2 +第47頁第47頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-16 求方程組 通解解法二：用rref求解A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;B=A b;C=rref(B) %求增廣矩陣行最簡形，得最簡同解方程組運(yùn)營后結(jié)果顯示為：第48頁第48頁4.3 線性方程組求解4.3.3 求非齊次線性方程組通解例4-16 求方程組 通解解法二：用rref求

25、解運(yùn)營后結(jié)果顯示為：相應(yīng)齊次方程組基礎(chǔ)解系為：非齊次方程組特解為：因此，原方程組通解為：X=k11+k22 + *。第49頁第49頁4.4 特性值與二次型工程技術(shù)中一些問題，如振動問題和穩(wěn)定性問題，常歸結(jié)為求一個方陣特性值和特性向量。4.4.1 特性值與特性向量求法設(shè)A為n階方陣，假如數(shù)“”和n維列向量x使得關(guān)系式Ax = x成立，則稱為方陣A特性值，非零向量x稱為A相應(yīng)于特性值特性向量。第50頁第50頁4.4 特性值與二次型4.4.1 特性值與特性向量求法在MATLAB中，計算矩陣A特性值和特性向量函數(shù)是eig(A)，慣用調(diào)用格式有3種：(1) E=eig(A)：求矩陣A所有特性值，構(gòu)成向量

26、E。(2) V, D=eig(A)：求矩陣A所有特性值，構(gòu)成對角陣D，并求A特性向量構(gòu)成V列向量。(3) V, D=eig(A,nobalance)：與第2種格式類似，但第2種格式中先對A作相同變換后求矩陣A特性值和特性向量，而格式3直接求矩陣A特性值和特性向量第51頁第51頁4.4 特性值與二次型4.4.1 特性值與特性向量求法例4-17 求矩陣 特性值和特性向量解：A=-2 1 1;0 2 0;-4 1 3;V,D=eig(A)結(jié)果顯示：即：特性值-1相應(yīng)特性向量(-0.7071 0 -0.7071)T特性值2相應(yīng)特性向量(-0.2425 0 -0.9701)T和(0.3015 0.9045 0.3015)T第52頁第52頁4.4 特性值與二次型4.4.1 特性值與特性向量求法例4-18 求矩陣 特性值和特性向量解： A=-1 1 0;-4 3 0;1 0 2;V,D=eig(A)結(jié)果顯示為闡明：當(dāng)特性值為1 (二重根)時，相應(yīng)特性向量都是k (0.4082 0