1、控制工程中的程序設(shè)計(jì)授課教師:馮肖亮E-mail: 河南工業(yè)大學(xué) 電氣工程學(xué)院1第1頁，共55頁。數(shù)值問題求解第四章2第2頁，共55頁。第四章 數(shù)值問題求解（上）4.1 特殊矩陣4.2 矩陣分析4.3 矩陣分解4.4 秩與線性相關(guān)性4.6 線性方程組的求解3第3頁，共55頁。4.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) 第4頁，共55頁。對角元素矩陣，

2、使用 diag( ) 命令兩類使用規(guī)則： (1) 構(gòu)造對角矩陣已知向量 V 生成對角矩陣 A ： A=diag(V)設(shè)V為具有m個元素的向量，diag(V)將產(chǎn)生一個mm對角矩陣，其主對角線元素即為向量V的元素，其他元素為0。生成主對角線上第k條對角線為向量V的矩陣 A ： A=diag(V,k)其功能是產(chǎn)生一個nn(n=m+abs(k)對角陣，其第k條對角線的元素即為向量V的元素，其他元素為0 ；k可正，可負(fù)。4.1 特殊矩陣第5頁，共55頁。(2) 提取矩陣的對角線元素 已知矩陣 A 提取對角元素列向量 V ： Vdiag(A)設(shè)A為mn矩陣，diag(A)函數(shù)用于提取矩陣A主對角線元素產(chǎn)

3、生一個具有min(m,n)個元素的列向量。已知矩陣A，提取主對角線上第k條對角線為向量V ： Vdiag(A,k)其功能是提取矩陣A主對角線上第k條對角線的元素。4.1 特殊矩陣第6頁，共55頁。 例1：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) % 主對角線上第 2條對角線為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第

4、7頁，共55頁。 例4.2 生成三對角矩陣: 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第8頁，共55頁。 例4.3 先建立55矩陣A，然后將A的第1行元素乘以1，第2行乘以2，第5行乘以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;D=diag(1,2,3,4,5); D*A %用D左乘A，對A的每行乘以一個指定常數(shù)說明：如果要對A的每行元素乘以一個常數(shù)，可以通過用一個對角陣左

5、乘矩陣A。9第9頁，共55頁。 2. 矩陣的三角陣 (1)下三角矩陣求矩陣A的下三角陣的MATLAB函數(shù)是tril(A)。求矩陣A第-k條對角線以下的下三角陣的MATLAB函數(shù)是tril(A,-k)。 (2)上三角矩陣在MATLAB中，提取矩陣A的上三角矩陣的函數(shù)是triu(A) ，其用法與提取下三角矩陣的函數(shù)tril(A) 完全相同。同樣地，triu(A,k)可求矩陣A第k條對角線以上的上三角陣10第10頁，共55頁。 A=1 2 -2;1 1 1;2 2 1; b=diag(A); bans = 1 1 1 triu(A)ans = 1 2 -2 0 1 1 0 0 1 tril(A)an

6、s = 1 0 0 1 1 0 2 2 1求矩陣A的上三角矩陣、對角陣和下三角矩陣 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第11頁，共55頁。4.2.1 矩陣結(jié)構(gòu)變換1. 矩陣的轉(zhuǎn)置轉(zhuǎn)置運(yùn)算符是單撇號()。2. 矩陣的旋轉(zhuǎn)矩陣的旋轉(zhuǎn)利用函數(shù)rot90(A,k)，功能是將矩陣A旋轉(zhuǎn)90的k倍，當(dāng)k為1時可省略。3. 矩陣的左右翻轉(zhuǎn)對矩陣A實(shí)施左右翻轉(zhuǎn)的函數(shù)是fliplr(A)。4. 矩陣的上下翻轉(zhuǎn)對矩陣A實(shí)施上下翻轉(zhuǎn)的函數(shù)是flipud(A)。4.2 矩陣分析12第12頁，共55頁。 4.2.2 矩陣

7、的逆與偽逆1. 矩陣的逆若A為n*n的非奇異方陣，則C=A-1,求解命令為C=inv(A)例4.4 ：求取4階Hilbert矩陣的逆矩陣 format long; H=hilb(4); H1=inv(H)H1 = 1.0e+003 * 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999 -0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984 0.23999999999998 -2.69999999999976 6.47999

8、999999940 -4.19999999999961 -0.13999999999999 1.67999999999984 -4.19999999999961 2.7999999999997413第13頁，共55頁。 4.2.2 矩陣的逆與偽逆2. 矩陣的偽逆MATLAB中，當(dāng)矩陣的逆不存在時，會定義它的偽逆。求一個矩陣偽逆的函數(shù)是pinv(A)。例4.5 求A的偽逆，并將結(jié)果送B。 命令如下：A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)14第14頁，共55頁。 例4.6 用求逆矩陣的方法解線性方程組AX=b。命令如下：A=1,2,3;1,4,9;1,8,27; b

9、=5, -2,6; x=inv(A)*b一般情況下，用左除比求矩陣的逆的方法更有效，即x=Ab。15第15頁，共55頁。 4.2.3 方陣的行列式求方陣A所對應(yīng)的行列式的值的函數(shù)是det(A)。4.2.4 矩陣的秩和跡設(shè)A是一組向量，A的極大無關(guān)組中向量的個數(shù)為A的秩MATLAB中，求矩陣秩的函數(shù)是rank(A)。方陣A的主對角線元素之和為A的跡MATLAB中，求矩陣的跡的函數(shù)是trace(A)。例如， X=2 2 3;4 5 -6;7 8 9，求X的行列式、秩和跡16第16頁，共55頁。17定理 (克萊姆Cramer法則)的系數(shù)行列式不等于零，即那么線性方程組(1)有解，并且解是唯一的，解可

10、以表示為如果線性方程組第17頁，共55頁。18其中 是把系數(shù)行列式 D 中第 j 列的元素用方程組右端的常數(shù)項(xiàng)代替后所得到的 n 階行列式，即注意：在利用克萊姆法則解方程組時，(1)方程組中方程的個數(shù)與未知數(shù)的個數(shù)必須相等；(2)系數(shù)行列式不能等于零. 第18頁，共55頁。 例4.7 用克萊姆(Cramer)方法求解線性方程組Dx=b。程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; %定義系數(shù)矩陣b=4;6;12;6; %定義常數(shù)項(xiàng)向量D1=b,D(:,2:4); %用方程組的右端向量置換D的第1列D2=D(:,1:1),b,D(:,3:4); %用方程

11、組的右端向量置換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,x419第19頁，共55頁。 4.2.5 向量的范數(shù)向量V的范數(shù)為： MATLAB中，求向量范數(shù)的函數(shù)是norm(V , R), norm(V)。4.2.4 矩陣的范數(shù)矩陣A 的范數(shù)為：MATLAB中，求矩陣范數(shù)的函數(shù)是norm(A, R), norm(A) 。20第20頁，共

12、55頁。 例4.8 已知V，求V的3種范數(shù)。命令如下：V=-1,1/2,1;v1=norm(V,1) %求V的1范數(shù)sum(abs(Vi)v2=norm(V) %求V的2范數(shù)sqrt(sum(Vi.2)vinf=norm(V,inf) %求范數(shù)max(abs(Vi)例4.9 求矩陣A的三種范數(shù)。命令如下：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范數(shù)a2=norm(A) %求A的2范數(shù)ainf=norm(A,inf) %求A的范數(shù)21第21頁，共55頁。 4.2.7

13、矩陣的特征值與特征向量MATLAB中，計(jì)算矩陣A的特征值和特征向量的函數(shù)是eig(A)，常用的調(diào)用格式有2種：(1)E=eig(A) 求矩陣A的全部特征值，構(gòu)成向量E。(2)V,D=eig(A) 求矩陣A的全部特征值，構(gòu)成對角陣D，并求A的特征向量構(gòu)成V的列向量。例： 求A的特征值和特征向量。命令如下：A=1,2,2;1,-1,1;4,-12,1;E=eig(A)V,D=eig(A)22第22頁，共55頁。例4.10 求對稱矩陣A的特征值和特征向量。命令如下：A=2,1,4,6;1,2,1,5;4,1,3,4;6,5,4,2;Q,D=eig(A)%若執(zhí)行如下命令，我們會有新的發(fā)現(xiàn)Q*D*Qan

14、s = 2.0000 1.0000 4.0000 6.0000 1.0000 2.0000 1.0000 5.0000 4.0000 1.0000 3.0000 4.0000 6.0000 5.0000 4.0000 2.0000結(jié)果與A相等，說明確實(shí)將A分解為了QDQ的乘積。23第23頁，共55頁。4.3 矩陣分解4.3.1 LU分解矩陣的三角分解又稱LU分解，它的目的是將一個矩陣分解成一個下三角矩陣L和一個上三角矩陣U的乘積，即A=LU。Matlab使用函數(shù)lu實(shí)現(xiàn)LU分解，其格式有2種：L,U = lu(A)其中U為上三角陣，L為下三角陣或其變換形式，滿足LU=A。L,U,P = lu(

15、A)其中，U為上三角陣，L為下三角陣，P為單位矩陣的行變換矩陣，滿足LU=PA。第24頁，共55頁。4.3 矩陣分解例4-11 實(shí)現(xiàn)矩陣A的LU分解命令如下：A=1 2 3;4 5 6;7 8 9;L,U=lu(A)L,U,P=lu(A)第25頁，共55頁。4.3 矩陣分解4.3.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與上

16、相同；若A非正定，則p為正整數(shù)，R是有序的上三角陣。第26頁，共55頁。4.3 矩陣分解例4-12 實(shí)現(xiàn)4階帕斯卡(pascal)矩陣的Cholesky分解命令如下：A=pascal(4) %產(chǎn)生4階pascal矩陣R,p=chol(A)第27頁，共55頁。4.3 矩陣分解4.3.3 QR分解將矩陣A分解成一個正交矩陣Q與一個上三角矩陣R的乘積A=QR，稱為QR分解。Matlab使用函數(shù)qr實(shí)現(xiàn)QR分解，其格式有兩種： Q,R = qr(A)其中，求得正交矩陣Q和上三角陣R，R的對角線元素按大小降序排列，滿足AE=QR。 Q,R,E = qr(A)其中，求得正交矩陣Q和上三角陣R，E為單位矩陣

17、的變換形式，R的對角線元素按大小降序排列，滿足AE=QR。第28頁，共55頁。4.3 矩陣分解例4-13實(shí)現(xiàn)矩陣A的QR分解 命令如下：A = 1 2 3;4 5 6; 7 8 9; 10 11 12;Q,R = qr(A)第29頁，共55頁。4.4 秩與線性相關(guān)性矩陣A的秩是指矩陣A中最高階非零子式的階數(shù)，或是矩陣線性無關(guān)的行數(shù)與列數(shù)；向量組的秩通常由該向量組構(gòu)成的矩陣來計(jì)算。求行階梯矩陣及向量組的基行階梯使用初等行變換，矩陣的初等行變換有三條：(1) 交換兩行ri rj (第i、第j兩行交換)(2) 第i行的k倍kri(3) 第i行的k倍加到第j行上去rj+kri通過這三條變換可以將矩陣化

18、成行最簡形，從而找出列向量組的一個最大無關(guān)組。第30頁，共55頁。4.4 秩與線性相關(guān)性4.4.1 求行階梯矩陣及向量組的基Matlab將矩陣化成行最簡形的命令是rref，其格式有兩種形式：R = rref(A) R 是A的行最簡形矩陣（非零行的第一個非零元為1；且這些非零元所在的列的其他元素都是零）R,jb = rref(A)R 是A的行最簡行矩陣，jb是一個向量，其含義為： r = length(jb)為A的秩； A(:, jb)為A的列向量基； jb中元素表示基向量所在的列。第31頁，共55頁。4.4 秩與線性相關(guān)性4.4.1 求行階梯矩陣及向量組的基例4-14 求向量組a1=(1,-2

19、,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 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為向量組的一個基。第32頁，共55頁。4.5 線性方程的組的求解我們將線性方程的求解分為兩類：一類是方程組求唯一解或求特解，另一類是方程組求無窮解即通解?？梢酝ㄟ^系數(shù)矩陣的秩來判斷：若系數(shù)矩陣的秩r=n(n為方程組中未知變量的個數(shù))，則有唯一解；

20、若系數(shù)矩陣的秩rn，則可能有無窮解；線性方程組的無窮解 = 對應(yīng)齊次方程組的通解+非齊次方程組的一個特解；其特解的求法屬于解的第一類問題，通解部分屬第二類問題。第33頁，共55頁。4.5 線性方程的組的求解4.5.1 求線性方程組的唯一解或特解1利用矩陣除法求線性方程組AX=b的特解1.1 當(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。

21、第34頁，共55頁。4.5 線性方程的組的求解例4-15 求方程組 的解。命令如下：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 %求解第35頁，共55頁。4.5 線性方程的組的求解4.5.1 求線性方程組的唯一解或特解1利用矩陣除法求線性方程組AX=b的特解1.2 用函數(shù)rref求解：C=A,B%由系數(shù)矩陣和常數(shù)列構(gòu)成增廣矩陣CR=rref(C)%將C化成行最簡行，其中R的最后一列元素就是所求之解。第36頁，共55頁。4.5 線性方程的組的求解例4-16 求方程組 的一

22、個特解命令1：高斯消去法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 一個特解近似值命令2：函數(shù)rrefC=A,B; %構(gòu)成增廣矩陣R=rref(C)由此得解向量X=1.2500 0.2500 0 0(一個特解)。第37頁，共55頁。4.5 線性方程的組的求解例4-17 解方程組 (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=Ab

23、1y=Ab2解法2：將兩個方程組連在一起求解：Az = bb=2 3;-3 4;1 -5z=Ab得該線性代數(shù)方程組的解：(1) x1= -3.8， x2= 1.4， x3= 7.2； (2) y1= -3.6， y2= 2.2， y3= 4.4得該線性代數(shù)方程組的解：(1) x1= -3.8， x2= 1.4， x3= 7.2； (2) y1= -3.6， y2= 2.2， y3= 4.4第38頁，共55頁。4.5 線性方程的組的求解4.5.1 求線性方程組的唯一解或特解2利用LU、QR和cholesky分解求方程組的解(1) LU分解：LU分解又稱Gauss消去分解，可把任意方陣分解為下三角

24、矩陣和上三角矩陣的乘積。即A=LU，L為下三角陣，U為上三角陣。則： A*X=b變成L*U*X=b所以X=U(Lb), 這樣可以大大提高運(yùn)算速度。第39頁，共55頁。4.5 線性方程的組的求解例4-18 求方程組 的一個特解。解：A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;D=det(A)L,U=lu(A)X=U(LB)第40頁，共55頁。4.5 線性方程的組的求解4.5.1 求線性方程組的唯一解或特解2利用LU、QR和cholesky分解求方程組的解(2) Cholesky分解：若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置的乘積，即：A =

25、 R*R，其中R為上三角陣則方程A*X=b變成R*R*X = b所以X=R(Rb) 。第41頁，共55頁。4.5 線性方程的組的求解例4-18 求方程組 的一個特解。解：A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;R=chol(A)X=R(RB)第42頁，共55頁。4.5 線性方程的組的求解4.5.1 求線性方程組的唯一解或特解2利用LU、QR和cholesky分解求方程組的解(3) QR分解：對于任何長方矩陣A，都可以進(jìn)行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A = QR則方程A*X=b變成QRX=b所以X=R(Qb) 。說明：這三種分解，在求解大

26、型方程組時很有用。其優(yōu)點(diǎn)是運(yùn)算速度快、可以節(jié)省磁盤空間、節(jié)省內(nèi)存。第43頁，共55頁。4.5 線性方程的組的求解例4-18 求方程組 的一個特解。解：A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;Q,R=qr(A)X=R(QB)第44頁，共55頁。4.5 線性方程的組的求解4.3.2 求線性齊次方程組的通解在Matlab中，函數(shù)null用來求解零空間，即滿足AX=0的解空間，實(shí)際上是求出解空間的一組基(基礎(chǔ)解系)。其格式有兩種為：z = null(A)z的列向量為方程組的正交規(guī)范基，滿足z z = I。z = null(A,r)z的列向量是方程AX=0的有理基第45頁，共55

27、頁。4.5 線性方程的組的求解例4-18 求方程組 的通解：解法1：函數(shù)nullA=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;format rat %指定有理式格式輸出B=null(A,r) %求解空間的有理基第46頁，共55頁。4.5 線性方程的組的求解4.5.3 求非齊次線性方程組的通解非齊次線性方程組需要先判斷方程組是否有解，若有解，再去求通解。因此，步驟為：第1步: 判斷AX=b是否有解，若有解則進(jìn)行第二步第2步: 求AX=b的一個特解第3步: 求AX=0的通解第4步: AX=b的通解= AX=0的通解+AX=b的一個特解第47頁，共55頁。4.3 線性方程的組的求解例

28、4-19 求解方程組解：在Matlab中建立M文件如下：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 %判斷無解end第48頁，共55頁。4.5 線性方程的組的求解例4-19 求解方程組解：在Matlab中建立函數(shù)M文件如下：function X=jie(A,b)m,n