1、控制工程中的程序设计控制工程中的程序设计授课教师:冯肖亮E-mail: 河南工业大学 电气工程学院1数值问题求解数值问题求解第四章第四章2第四章第四章 数值问题求解数值问题求解（上）（上） 4.1 4.1 特殊矩阵特殊矩阵 4.2 4.2 矩阵分析矩阵分析 4.3 4.3 矩阵分解矩阵分解 4.4 4.4 秩与线性相关性秩与线性相关性 4.6 4.6 线性方程组的求解线性方程组的求解34.1 4.1 特殊矩阵特殊矩阵n 数值矩阵的输入数值矩阵的输入n 零矩阵、幺矩阵及单位矩阵零矩阵、幺矩阵及单位矩阵 生成生成n n n n方阵：方阵： A=zeros(n), B=ones(n), C=eye(
2、n)A=zeros(n), B=ones(n), C=eye(n) 生成生成m m n n矩阵：矩阵： A=zeros(m,n), B=ones(m,n), C=eye(m,n)A=zeros(m,n), B=ones(m,n), C=eye(m,n) 生成和矩阵生成和矩阵B B同样位数的矩阵：同样位数的矩阵： A=zeros(size(B) A=zeros(size(B) n 对角元素对角元素矩阵，矩阵， 使用使用 diagdiag( ) ( ) 命令命令n 两类使用规则：两类使用规则： ( (1 1) ) 构造构造对角矩阵对角矩阵已知向量已知向量 V V 生成对角矩阵生成对角矩阵 A A
3、： A=A=diagdiag(V)(V)设设V V为具有为具有m m个元素的向量，个元素的向量，diagdiag(V)(V)将产生一个将产生一个m mm m对角矩阵，其对角矩阵，其主对角线元素即为向量主对角线元素即为向量V V的的元素，其他元素为元素，其他元素为0 0。生成主对角线上第生成主对角线上第k k条对角线为向量条对角线为向量V V的矩阵的矩阵 A A ： A= A=diagdiag( (V,kV,k) )其其功能是产生一个功能是产生一个n nn n(n=(n=m+absm+abs(k)(k)对角阵，其第对角阵，其第k k条对角线的元素条对角线的元素即为向量即为向量V V的的元素元素，
4、其他元素为，其他元素为0 0 ；k k可正，可负。可正，可负。4.1 4.1 特殊矩阵特殊矩阵(2) (2) 提取矩阵的对角线提取矩阵的对角线元素元素 已知矩阵已知矩阵 A A 提取提取对角元素对角元素列向量列向量 V V ： V Vdiagdiag(A)(A)设设A A为为m mn n矩阵，矩阵，diagdiag(A)(A)函数用于提取矩阵函数用于提取矩阵A A主对角线元素产生一个主对角线元素产生一个具有具有min(min(m,nm,n) )个元素的列向量。个元素的列向量。已知矩阵已知矩阵A A，提取主对角线上第，提取主对角线上第k k条对角线为向量条对角线为向量V V ： V Vdiagd
5、iag( (A,kA,k) )其功能是其功能是提取提取矩阵矩阵A A主主对角线上第对角线上第k k条对角线的元素。条对角线的元素。4.1 4.1 特殊矩阵特殊矩阵 例例1 1：diagdiag( )( )函数的不同调用格式函数的不同调用格式 C=1 2 3; V= C=1 2 3; V=diagdiag(C) % (C) % 生成对角矩阵生成对角矩阵V =V = 1 0 0 1 0 0 0 2 0 0 2 0 0 0 3 0 0 3 V1= V1=diagdiag(V) % (V) % 将列向量通过转置变换成行向量将列向量通过转置变换成行向量V1 =V1 = 1 2 3 1 2 3 C=1 2
6、 3; V= C=1 2 3; V=diagdiag(C,2) % (C,2) % 主对角线上第主对角线上第 2 2条条对角线为对角线为C C的矩阵的矩阵V =V = 0 0 1 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 例例4.2 4.2 生成生成三对角矩阵三对角矩阵: : V= V=diagdiag(1 2 3 4)+(1 2 3 4)+diagdiag(2 3 4,1)+(2 3 4,1)+diagdiag(5 4 3,-1)(5 4 3,-1
7、)V =V = 1 2 0 0 1 2 0 0 5 2 3 0 5 2 3 0 0 4 3 4 0 4 3 4 0 0 3 4 0 0 3 41200523004340034 例例4.3 4.3 先先建立建立5 55 5矩阵矩阵A A，然后将，然后将A A的第的第1 1行元素乘以行元素乘以1 1，第，第2 2行乘以行乘以2 2，第，第5 5行乘以行乘以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,A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19
8、,21,3;11,18,25,2,19;2,19;D=D=diagdiag(1,2,3,4,5); (1,2,3,4,5); D D* *A %A %用用D D左乘左乘A A，对，对A A的每行乘以一个指定常数的每行乘以一个指定常数说明：说明：如果要对如果要对A A的的每每行行元素元素乘以一个常数，可以乘以一个常数，可以通过用一个对角通过用一个对角阵左乘阵左乘矩阵矩阵A A。9 2. 2. 矩阵的三角阵矩阵的三角阵 (1)(1)下三角矩阵下三角矩阵求求矩阵矩阵A A的下三角阵的的下三角阵的MATLABMATLAB函数是函数是triltril(A)(A)。求求矩阵矩阵A A第第-k -k条对角线
9、以下的条对角线以下的下三角阵的下三角阵的MATLABMATLAB函数函数是是triltril(A,(A,-k -k) )。 (2)(2)上三角矩阵上三角矩阵在在MATLABMATLAB中，提取矩阵中，提取矩阵A A的上三角矩阵的函数是的上三角矩阵的函数是triutriu(A(A) ) ，其，其用法与提取下三角矩阵的函数用法与提取下三角矩阵的函数triltril(A(A) ) 完全完全相同相同。同样地，同样地，triutriu( (A,A,k k) )可求可求矩阵矩阵A A第第k k条条对角线以上的上三对角线以上的上三角阵角阵10 A=1 2 -2;1 1 1;2 2 1; b=diag(A);
10、 bans = 1 1 1 triu(A)ans = 1 2 -2 0 1 1 0 0 1 tril(A)ans = 1 0 0 1 1 0 2 2 1122111221A求矩阵求矩阵A A的上三角矩阵、对角阵和下三角矩阵的上三角矩阵、对角阵和下三角矩阵 triutriu(A,1)(A,1)ansans = = 0 2 -2 0 2 -2 0 0 1 0 0 1 0 0 0 0 0 0 triltril(A,-1)(A,-1)ansans = = 0 0 0 0 0 0 1 0 0 1 0 0 2 2 0 2 2 04.2.1 4.2.1 矩阵结构变换矩阵结构变换1. 1. 矩阵的转置矩阵的转
11、置转置运算符是单撇号()。2. 2. 矩阵的旋转矩阵的旋转矩阵的旋转利用函数rot90(A,k)，功能是将矩阵A旋转90的k倍，当k为1时可省略。3. 3. 矩阵的左右翻转矩阵的左右翻转对矩阵A实施左右翻转的函数是fliplr(A)。4. 4. 矩阵的上下翻转矩阵的上下翻转对矩阵A实施上下翻转的函数是flipud(A)。4.2 4.2 矩阵分析矩阵分析12 4.2.2 4.2.2 矩阵的逆与伪逆矩阵的逆与伪逆1. 1. 矩阵的逆矩阵的逆若若A A为为n n* *n n的非奇异方阵，则的非奇异方阵，则C C= =A A-1 -1, ,求解命令为求解命令为C=inv(A)C=inv(A)例例4.4
12、 4.4 ：求取：求取4 4阶阶HilbertHilbert矩阵的逆矩阵矩阵的逆矩阵 format long; H=hilb(4); H1=inv(H) format long; H=hilb(4); H1=inv(H)H1 =H1 = 1.0e+003 1.0e+003 * * 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999 0.01600000000000 -0.11999999999999 0.23999999999998 -0.13999999999999 -0.11999999999999 1
13、.19999999999990 -2.69999999999976 1.67999999999984 -0.11999999999999 1.19999999999990 -2.69999999999976 1.67999999999984 0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961 0.23999999999998 -2.69999999999976 6.47999999999940 -4.19999999999961 -0.13999999999999 1.67999999999984 -4.1
14、9999999999961 2.79999999999974 -0.13999999999999 1.67999999999984 -4.19999999999961 2.7999999999997413 4.2.2 4.2.2 矩阵的逆与伪逆矩阵的逆与伪逆2. 2. 矩阵的伪逆矩阵的伪逆MATLAB中，当矩阵的逆不存在时，会定义它的伪逆。求一个矩阵伪逆的函数是pinv(A)。例例4.5 4.5 求求A A的伪逆，并将结果送的伪逆，并将结果送B B。 命令如下：A=3,1,1,1;1,3,1,1;1,1,3,1;B=pinv(A)14 例例4.6 4.6 用求逆矩阵的方法解用求逆矩阵的方法解线
15、性方程组线性方程组AX=bAX=b。命令如下：A=1,2,3;1,4,9;1,8,27; b=5, -2,6; x=inv(A)*b一般情况下，用左除比求矩阵的逆的方法更有效，即x=Ab。15 4.2.3 4.2.3 方阵的行列式方阵的行列式求方阵求方阵A A所对应的行列式的值的函数是所对应的行列式的值的函数是detdet(A)(A)。4.2.4 4.2.4 矩阵的秩和迹矩阵的秩和迹设A是一组向量，A的极大无关组中向量的个数为A A的秩的秩MATLABMATLAB中，求矩阵秩的函数是中，求矩阵秩的函数是rank(A)rank(A)。方阵A的主对角线元素之和为A A的迹的迹MATLABMATLA
16、B中，求矩阵的迹的函数是中，求矩阵的迹的函数是trace(A)trace(A)。例如，例如， X=2 2 3;4 5 -6;7 8 9X=2 2 3;4 5 -6;7 8 9，求，求X X的行列式、秩和迹的行列式、秩和迹1617定理定理 ( (克莱姆克莱姆CramerCramer法则法则) )1(22112222212111212111 nnnnnnnnnnbxaxaxabxaxaxabxaxaxa的系数行列式不等于零，即nnnnnnaaaaaaaaaD212222111211 ,0 那么线性方程组(1)有解，并且解是唯一的，解可以表示为如果线性方程组18.,332211DDxDDxDDxDD
17、xnn 其中 是把系数行列式 D 中第 j 列的元素用方程组右端的常数项代替后所得到的 n 阶行列式，即jDnnjnnjnnnjjnjjjaabaaaabaaaabaaD1,1,121,221,22111, 111, 111 注意：在利用克莱姆法则解方程组时，(1)方程组中方程的个数与未知数的个数必须相等；(2)系数行列式不能等于零. 例例4.7 4.7 用克莱姆用克莱姆(Cramer)(Cramer)方法求解方法求解线性方程组线性方程组DxDx=b=b。程序如下：程序如下：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; %D=2,2,-1,1;4,3,-1,2;
18、8,5,-3,4;3,3,-2,2; %定义系数矩阵定义系数矩阵b=4;6;12;6; %b=4;6;12;6; %定义常数项向量定义常数项向量D1=D1=b,Db,D(:,2:4); %(:,2:4); %用方程组的右端向量置换用方程组的右端向量置换D D的第的第1 1列列D2=D(:,1:1),D2=D(:,1:1),b,Db,D(:,3:4); %(:,3:4); %用方程组的右端向量置换用方程组的右端向量置换D D的第的第2 2列列D3=D(:,1:2),D3=D(:,1:2),b,Db,D(:,4:4); %(:,4:4); %用方程组的右端向量置换用方程组的右端向量置换D D的第的
19、第3 3列列D4=D(:,1:3),b; %D4=D(:,1:3),b; %用方程组的右端向量置换用方程组的右端向量置换D D的第的第4 4列列DD=DD=detdet(D);(D);x1=x1=detdet(D1)/DD;(D1)/DD;x2=x2=detdet(D2)/DD;(D2)/DD;x3=x3=detdet(D3)/DD;(D3)/DD;x4=x4=detdet(D4)/DD;(D4)/DD;x1,x2,x3,x4x1,x2,x3,x419 4.2.5 4.2.5 向量的范数向量的范数向量向量V V的范数为：的范数为： MATLABMATLAB中，中，求向量范数的求向量范数的函数函
20、数是是norm(V , Rnorm(V , R), norm(V), norm(V)。4.2.4 4.2.4 矩阵的范数矩阵的范数矩阵矩阵A A 的范数为：的范数为：MATLABMATLAB中，求矩阵范数的函数是中，求矩阵范数的函数是norm(A, R), norm(A) norm(A, R), norm(A) 。201/11,1,2,maxpnpipii niVvpVv 111max211max()maxnj nijiTni nijjAaAA AAa 例例4.8 4.8 已知已知V V，求，求V V的的3 3种范数。种范数。命令如下：命令如下：V=-1,1/2,1;V=-1,1/2,1;v1
21、=norm(V,1) v1=norm(V,1) % %求求V V的的11范数范数sum(abs(Vi)v2=norm(V) v2=norm(V) %求求V V的的22范数范数sqrt(sum(Vi.2)vinfvinf=norm(=norm(V,infV,inf) ) %求求范数范数max(abs(Vi)例例4.9 4.9 求矩阵求矩阵A A的三种范数。的三种范数。命令如下：命令如下：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=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;
22、10,12,19,21,3;11,18,25,2,19;a1=norm(A,1) %a1=norm(A,1) %求求A A的的11范数范数a2=norm(A) %a2=norm(A) %求求A A的的22范数范数ainfainf=norm(=norm(A,infA,inf) %) %求求A A的的范数范数21 4.2.7 4.2.7 矩阵的特征值与特征向量矩阵的特征值与特征向量MATLABMATLAB中，计算矩阵中，计算矩阵A A的特征值和特征向量的的特征值和特征向量的函数是函数是eigeig(A)(A)，常用的调用格式，常用的调用格式有有2 2种种：(1)E=(1)E=eigeig(A) (
23、A) 求矩阵求矩阵A A的全部特征值，构成向量的全部特征值，构成向量E E。(2)V,D=(2)V,D=eigeig(A) (A) 求矩阵求矩阵A A的全部特征值，构成对角阵的全部特征值，构成对角阵D D，并求并求A A的特征向量构成的特征向量构成V V的列向量的列向量。例：例： 求求A A的特征值和特征向量。的特征值和特征向量。命令如下：A=1,2,2;1,-1,1;4,-12,1;E=eig(A)V,D=eig(A)22例例4.10 4.10 求求对称矩阵对称矩阵A A的特征值的特征值和特征向量和特征向量。命令如下：命令如下：A A=2,1,4,6;1,2,1,5;4,1,3,4;6,5,
24、4,2;=2,1,4,6;1,2,1,5;4,1,3,4;6,5,4,2;Q,D=Q,D=eigeig(A)(A)%若执行如下命令，我们会有新的发现若执行如下命令，我们会有新的发现Q Q* *D D* *Q Qansans = = 2.0000 1.0000 4.0000 6.0000 2.0000 1.0000 4.0000 6.0000 1.0000 2.0000 1.0000 5.0000 1.0000 2.0000 1.0000 5.0000 4.0000 1.0000 3.0000 4.0000 4.0000 1.0000 3.0000 4.0000 6.0000 5.0000 4.
25、0000 2.0000 6.0000 5.0000 4.0000 2.0000结果与结果与A A相等，说明确实将相等，说明确实将A A分解为了分解为了QDQQDQ的乘积的乘积。234.3 4.3 矩阵分解矩阵分解n 4.3.1 4.3.1 LULU分解分解 矩阵的三角分解又称矩阵的三角分解又称LULU分解，它的目的是将一个矩分解，它的目的是将一个矩阵分解成一个下三角矩阵阵分解成一个下三角矩阵L L和一个上三角矩阵和一个上三角矩阵U U的乘的乘积，即积，即A=LUA=LU。 MatlabMatlab使用函数使用函数lu lu实现实现LULU分解，其分解，其格式有格式有2 2种：种： L,U =
26、L,U = lu lu(A)(A)其中其中U U为上三角阵，为上三角阵，L L为下三角阵或其变换形式为下三角阵或其变换形式，满足，满足LU=ALU=A。 L,U,P = L,U,P = lu lu(A)(A)其中，其中，U U为上三角阵，为上三角阵，L L为下三角阵，为下三角阵，P P为单位矩为单位矩阵的行变换矩阵，满足阵的行变换矩阵，满足LU=PALU=PA。4.3 4.3 矩阵分解矩阵分解例例4-11 4-11 实现矩阵实现矩阵A A的的LULU分解分解命令如下：命令如下：A=1 2 3;4 5 6;7 8 9;A=1 2 3;4 5 6;7 8 9;L,U=lu(A)L,U=lu(A)L
27、,U,P=lu(A)L,U,P=lu(A)4.3 4.3 矩阵分解矩阵分解n 4.3.2 4.3.2 CholeskyCholesky分解分解 如果如果A A为为n n阶阶对称正定矩阵对称正定矩阵，则存在一个，则存在一个实的非奇异实的非奇异上三角阵上三角阵R R，满足，满足RR* *R = AR = A，称为，称为CholeskyCholesky分解分解 MatlabMatlab使用函数使用函数cholchol实现实现CholeskyCholesky分解，其格式为：分解，其格式为： R = R = cholchol(A)(A) 若若A A非正定，则产生错误信息。非正定，则产生错误信息。 R,p
28、R,p = = cholchol(A(A) )不不产生任何错误信息，若产生任何错误信息，若A A为正定阵，则为正定阵，则p=0p=0，R R与上相同；若与上相同；若A A非正定，则非正定，则p p为正整数，为正整数，R R是有序是有序的上三角阵。的上三角阵。4.3 4.3 矩阵分解矩阵分解n例例4-12 4-12 实现实现4 4阶帕斯卡阶帕斯卡(pascal)(pascal)矩阵的矩阵的CholeskyCholesky分解分解命令如下：命令如下：A=pascal(4) A=pascal(4) %产生产生4 4阶阶pascalpascal矩阵矩阵R,p=chol(A)R,p=chol(A)4.3
29、 4.3 矩阵分解矩阵分解n 4.3.3 4.3.3 QRQR分解分解 将矩阵将矩阵A A分解成一个正交矩阵分解成一个正交矩阵Q Q与一个上三角矩阵与一个上三角矩阵R R的的乘积乘积A=QRA=QR，称为，称为QRQR分解。分解。 MatlabMatlab使用函数使用函数qrqr实现实现QRQR分解，其分解，其格式有两种：格式有两种： Q,R = Q,R = qrqr(A(A) )其中，求得正交矩阵其中，求得正交矩阵Q Q和上三角阵和上三角阵R R，R R的对角线元素的对角线元素按大小降序排列，满足按大小降序排列，满足AE=QRAE=QR。 Q,R,E = Q,R,E = qrqr(A)(A)
30、其中，求得其中，求得正交矩阵正交矩阵Q Q和上三角阵和上三角阵R R，E E为单位矩阵的为单位矩阵的变换形式，变换形式，R R的对角线元素按大小降序排列，满足的对角线元素按大小降序排列，满足AE=QRAE=QR。4.3 4.3 矩阵分解矩阵分解n 例例4-134-13实现矩阵实现矩阵A A的的QRQR分解分解 命令如下：命令如下：A = 1 2 3;4 5 6; 7 8 9; 10 11 12;A = 1 2 3;4 5 6; 7 8 9; 10 11 12;Q,R = qr(A)Q,R = qr(A)4.4 4.4 秩与线性相关性秩与线性相关性n矩阵矩阵A A的秩是指矩阵的秩是指矩阵A A中
31、最高阶非零子式的阶数，或中最高阶非零子式的阶数，或是矩阵线性无关的行数与列数；是矩阵线性无关的行数与列数；向量组的秩通常由向量组的秩通常由该向量组构成的矩阵来计算该向量组构成的矩阵来计算。n求求行阶梯矩阵及向量组的基行阶梯矩阵及向量组的基行阶梯使用初等行变换，矩阵的初等行变换有三条：行阶梯使用初等行变换，矩阵的初等行变换有三条： (1) (1) 交换两行交换两行r ri i r rj j ( (第第i i、第、第j j两行交换两行交换) ) (2) (2) 第第i i行的行的k k倍倍krkri i (3) (3) 第第i i行的行的k k倍加到第倍加到第j j行上去行上去r rj j+ +k
32、rkri i通过这三条变换可以将矩阵化成通过这三条变换可以将矩阵化成行最简形行最简形，从而找，从而找出出列向量组的一个最大无关组列向量组的一个最大无关组。4.4 4.4 秩与线性相关性秩与线性相关性n 4.4.1 4.4.1 求行阶梯矩阵及向量组的基求行阶梯矩阵及向量组的基 MatlabMatlab将矩阵化成行最简形的命令是将矩阵化成行最简形的命令是rrefrref，其格式有，其格式有两种形式：两种形式： R = R = rrefrref(A)(A) R R 是是A A的行最的行最简形简形矩阵（矩阵（非零行的第一个非零元为1；且这些非零元所在的列的其他元素都是零） R,jbR,jb = = r
33、refrref(A)(A)R R 是是A A的行最简行矩阵，的行最简行矩阵，jb jb是一是一个向量个向量，其含义，其含义为为： r r = length(= length(jb jb) )为为A A的秩的秩； A A(:, (:, jb jb) )为为A A的列向量基的列向量基； jb jb中元素表示基向量所在的列。中元素表示基向量所在的列。4.4 4.4 秩与线性相关性秩与线性相关性n 4.4.1 4.4.1 求行阶梯矩阵及向量组的基求行阶梯矩阵及向量组的基例例4-14 4-14 求向量组求向量组a a1 1=(1,-2,2,3)=(1,-2,2,3)，a a2 2=(-2,4,-1,3)
34、=(-2,4,-1,3)，a a3 3=(-1,2,0,3)=(-1,2,0,3)，a a4 4=(0,6,2,3)=(0,6,2,3)，a a5 5=(2,-6,3,4)=(2,-6,3,4)的一个最大无关组的一个最大无关组a1=1 -2 2 3; a2=-2 4 -1 3; a3=-1 2 0 3;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;a4=0 6 2 3; a5=2 -6 3 4;A=a1 a2 a3 a4 a5A=a1 a2 a3 a4 a5R,jb=rref(A)R,jb=rref(A)A(:,j
35、b)A(:,jb)即：即：a1 a2 a4a1 a2 a4为向量组的一个基为向量组的一个基。4.5 4.5 线性方程的组的求解线性方程的组的求解我们将线性方程的求解分为两类：我们将线性方程的求解分为两类：一类是方程组求唯一一类是方程组求唯一解解或求特解或求特解，另一类是方程组求无穷解即通解另一类是方程组求无穷解即通解。可以通过。可以通过系系数矩阵的秩数矩阵的秩来判断：来判断：l 若系数矩阵的秩若系数矩阵的秩r=n(nr=n(n为方程组中未知变量的个数为方程组中未知变量的个数) )，则，则有唯一解；有唯一解；l 若系数矩阵的秩若系数矩阵的秩rnrn，则可能有无穷解；，则可能有无穷解；线性方程组的
36、无穷解线性方程组的无穷解 = = 对应齐次方程组的通解对应齐次方程组的通解+ +非齐非齐次方程组的一个特解次方程组的一个特解；其特解的求法属于解的第一类；其特解的求法属于解的第一类问题，通解部分属第二类问题。问题，通解部分属第二类问题。4.5 4.5 线性方程的组的求解线性方程的组的求解4.5.1 4.5.1 求线性方程组的唯一解或求线性方程组的唯一解或特解特解1 1利用矩阵除法求线性方程组利用矩阵除法求线性方程组AX=bAX=b的特解的特解1.1 1.1 当当系数矩阵系数矩阵A A为为N N* *N N的方阵时，的方阵时，MATLABMATLAB会自行用高会自行用高斯消去法求解线性代数方程组
37、。斯消去法求解线性代数方程组。 若若右端项右端项b b为为N N* *1 1的列向量的列向量，则，则x=Abx=Ab可获得方程可获得方程组的数值解组的数值解x x（N N* *1 1的列向量）；的列向量）； 若若右端项右端项b b为为N N* *MM的矩阵的矩阵，则，则x=Abx=Ab可同时获得同可同时获得同一系数矩阵一系数矩阵A A、MM个方程组数值解个方程组数值解x x（为（为N N* *MM的矩的矩阵），即阵），即x(:,j)=Ab(:,j)x(:,j)=Ab(:,j)，j=1,2,Mj=1,2,M。4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-15 4-15 求方程组求方
38、程组 的解。的解。 命令如下：命令如下：A=5 6 0 0 0A=5 6 0 0 0 1 5 6 0 0 1 5 6 0 0 0 1 5 6 0 0 1 5 6 0 0 0 1 5 6 0 0 1 5 6 0 0 0 1 5; 0 0 0 1 5; B=1 0 0 0 1;B=1 0 0 0 1;R_A=rank(A) R_A=rank(A) %求秩求秩X=AB X=AB %求解求解121232343454556156056056051xxxxxxxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解4.5.1 4.5.1 求线性方程组的唯一解或特解求线性方程组的唯一解或特解1 1
39、利用矩阵除法求线性方程组利用矩阵除法求线性方程组AX=bAX=b的特解的特解1.2 1.2 用函数用函数rrefrref求解：求解：C=A,BC=A,B%由系数矩阵和常数列构成增广矩阵由系数矩阵和常数列构成增广矩阵C CR=rref(C)R=rref(C) %将将C C化成行最简行，其中化成行最简行，其中R R的最后一列的最后一列元素就是所求之解。元素就是所求之解。4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-16 4-16 求方程组求方程组 的一个特解的一个特解命令命令1 1：高斯消去法：高斯消去法A A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;=1 1 -
40、3 -1;3 -1 -3 4;1 5 -9 -8;B=1 4 0;B=1 4 0;X=AB %X=AB %由于系数矩阵不满秩由于系数矩阵不满秩, , 该解法可能存在误差该解法可能存在误差X = 0 0 -0.5333 0.6000 X = 0 0 -0.5333 0.6000 一个特解一个特解近似值近似值命令命令2 2：函数：函数rrefrrefC=A,B; %C=A,B; %构成增广矩阵构成增广矩阵R=R=rrefrref(C)(C)由此得解向量由此得解向量X=1.2500 0.2500 0 0(X=1.2500 0.2500 0 0(一个特解一个特解) )。1234123412343133
41、445980 xxxxxxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-17 4-17 解方程组解方程组 (1)(1)AxAx= =b b1 1；(2)(2)AyAy= =b b2 2 解法解法1 1：分别解方程组：分别解方程组 (1)(1)AxAx= =b b1 1；(2)(2)AyAy= =b b2 2A=1 -1 1;5 -4 3;2 1 1;A=1 -1 1;5 -4 3;2 1 1;b1=2;-3;1; b2=3;4;-5;b1=2;-3;1; b2=3;4;-5;x=Ab1x=Ab1y=Ab2y=Ab2 解法解法2 2：将两个方程组连在一起求解：将两个方程
42、组连在一起求解：AzAz = b = bb=2 3;-3 4;1 -5b=2 3;-3 4;1 -5z=Abz=Ab111125432321131xxx 111154322113yyy 345得该线性代数方程组的解：得该线性代数方程组的解：(1) x1= -3.8(1) x1= -3.8， x2= 1.4x2= 1.4， x3= 7.2x3= 7.2； (2) y1= -3.6(2) y1= -3.6， y2= 2.2y2= 2.2， y3= 4.4y3= 4.4得该线性代数方程组的解：得该线性代数方程组的解：(1) x1= -3.8(1) x1= -3.8， x2= 1.4x2= 1.4，
43、x3= 7.2x3= 7.2； (2) y1= -3.6(2) y1= -3.6， y2= 2.2y2= 2.2， y3= 4.4y3= 4.44.5 4.5 线性方程的组的求解线性方程的组的求解4.5.1 4.5.1 求线性方程组的唯一解或特解求线性方程组的唯一解或特解2 2利用利用LULU、QRQR和和choleskycholesky分解求方程组的解分解求方程组的解(1 (1) LU) LU分解：分解：LULU分解又称分解又称GaussGauss消去分解，可把消去分解，可把任意方阵任意方阵分解为下三分解为下三角矩阵和上三角矩阵的乘积。角矩阵和上三角矩阵的乘积。即即A=LUA=LU，L L为
44、下三角阵，为下三角阵，U U为上三角阵。则：为上三角阵。则： A A* *X=bX=b变成变成L L* *U U* *X=bX=b所以所以X=U(Lb), X=U(Lb), 这样可以大大提高运算速度。这样可以大大提高运算速度。4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-18 4-18 求方程组求方程组 的一个特解。的一个特解。解：解：A=4 2 -1;3 -1 2;11 3 0;A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;B=2 10 8;D=det(A)D=det(A)L,U=lu(A)L,U=lu(A)X=U(LB)X=U(LB)123123124223
45、2101138xxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解4.5.1 4.5.1 求线性方程组的唯一解或特解求线性方程组的唯一解或特解2 2利用利用LULU、QRQR和和choleskycholesky分解求方程组的解分解求方程组的解(2) (2) CholeskyCholesky分解：分解：若若A A为为对称正定矩阵对称正定矩阵，则，则CholeskyCholesky分解可将矩阵分解可将矩阵A A分解成分解成上三角矩阵和其转置的乘积，上三角矩阵和其转置的乘积，即：即：A = RA = R* *R R，其中，其中R R为上三角阵为上三角阵则方程则方程A A* *X=bX
46、=b变成变成RR* *R R* *X = bX = b所以所以X=R(Rb)X=R(Rb) 。4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-18 4-18 求方程组求方程组 的一个特解。的一个特解。解：解：A=4 2 -1;3 -1 2;11 3 0;A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;B=2 10 8;R=chol(A)R=chol(A)X=R(RB)X=R(RB)1231231242232101138xxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解4.5.1 4.5.1 求线性方程组的唯一解或特解求线性方程组的唯一解或特解2 2
47、利用利用LULU、QRQR和和choleskycholesky分解求方程组的解分解求方程组的解(3) QR(3) QR分解：分解：对于对于任何长方矩阵任何长方矩阵A A，都可以进行，都可以进行QRQR分解，其中分解，其中Q Q为正为正交矩阵，交矩阵，R R为上三角矩阵的初等变换形式，为上三角矩阵的初等变换形式，即：即：A = QRA = QR则方程则方程A A* *X=bX=b变成变成QRX=bQRX=b所以所以X=R(Qb)X=R(Qb) 。说明：这三种分解，在求解大型方程组时很有用。其说明：这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。优点是运算速度
48、快、可以节省磁盘空间、节省内存。4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-18 4-18 求方程组求方程组 的一个特解。的一个特解。解：解：A=4 2 -1;3 -1 2;11 3 0;A=4 2 -1;3 -1 2;11 3 0;B=2 10 8;B=2 10 8;Q,R=qr(A)Q,R=qr(A)X=R(QB)X=R(QB)1231231242232101138xxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解4.3.2 4.3.2 求线性齐次方程组的通解求线性齐次方程组的通解在在MatlabMatlab中，函数中，函数nullnull用来求解零空间，
49、即满足用来求解零空间，即满足AX=0AX=0的解空间，实际上是求出解空间的一组基的解空间，实际上是求出解空间的一组基( (基础解系基础解系) )。其其格式有两种为格式有两种为：z = null(A)z = null(A)z z的列向量为方程组的正交规范基，满足的列向量为方程组的正交规范基，满足z z z z = = I I。z = null(z = null(A,rA,r) )z z的列向量是方程的列向量是方程AX=0AX=0的有理基的有理基4.5 4.5 线性方程的组的求解线性方程的组的求解例例4-18 4-18 求方程组求方程组 的的通解：通解：解法解法1 1：函数：函数nullnullA
50、 A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;format rat %format rat %指定有理式格式输出指定有理式格式输出B=null(B=null(A,rA,r) %) %求解空间的有理基求解空间的有理基1234123412342202220430 xxxxxxxxxxxx4.5 4.5 线性方程的组的求解线性方程的组的求解4.5.3 4.5.3 求非齐次线性方程组的通解求非齐次线性方程组的通解非齐次线性方程组需要先判断方程组是否有解，若有解，再非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。
