《第7章--Matlab与线性代数.doc》由会员分享，可在线阅读，更多相关《第7章--Matlab与线性代数.doc（71页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date第7章-Matlab与线性代数第6章 Matlab与线性代数第7章 Matlab在线性代数中的应用7.1 矩阵的基本函数运算矩阵的函数运算是矩阵运算中最为实用的部分。它主要包括矩阵的逆运算、矩阵的行列式、特征值运算、矩阵秩、迹的运算、求正交矩阵运算、向量组的最大无关组等。7.1.1 矩阵的逆运算矩阵的逆运算是矩阵运算中很重要的一种运算。它在线性代数及计算方法中都有很多
2、的论述。而在Matlab中，众多的复杂理论只变成了一个简单的命令inv或A (-1)。例如：a=2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5;inv(a)ans = -0.0471 0.5882 -0.2706 -0.9412 0.3882 -0.3529 0.4824 0.7647 -0.2235 0.2941 -0.0353 -0.4706 -0.0353 -0.0588 0.0471 0.29417.1.2. 矩阵的行列式运算矩阵的行列式的值可由det函数计算得出。例如：a1=det(a);a2=det(inv(a);a1*a2ans = 1.00007.1.
3、3. 迹函数矩阵所有对角线上元素的和称为矩阵的迹，在Matlab中可由trace函数计算得出。例如：trace(A)ans = 157.1.4 向量的点乘（内积）维数相同的两个向量的点乘可由函数dot得出。例如：X=-1 0 2;Y=-2 -1 1;Z=dot(X,Y)Z = 4还可用另一种算法：sum(X.*Y)ans = 47.1.5 向量叉乘在数学上，两个向量的叉乘是一个过两相交向量的交点且垂直于两向量所在平面的向量。在Matlab中，用函数cross实现。例如：c=cross(X,Y)c = 2 -3 17.1.6 混合积混合积由以上两函数实现。例如：求X(Yc)d=dot(X,cro
4、ss(Y,c)d = 147.2 符号矩阵的运算7.2.1 四则运算符号矩阵的四则运算同数值矩阵的四则运算完全相同。7.2.2 其它一些基本运算符号矩阵的其它一些基本运算也与数值矩阵的运算格式相同。这些运算包括矩阵的转置（）、行列式（det）、逆(inv)、秩（rank）、幂（）和指数（exp和expm等）运算。7.2.3 符号矩阵的简化符号工具箱中还提供了符号矩阵因式分解、展开、合并、简化及通分等符号操作函数。下面将一一进行介绍。7.2.3.1 因式分解 factor (s) 符号表达式因式分解函数，s为符号矩阵或符号表达式，常用于多项式的因式分解。例7-1：将x9-1分解因式。syms x
5、factor(x9-1) ans = (x-1)*(x2+x+1)*(x6+x3+1)例7-2：问为何值时，齐次方程组 有非0解？ syms kA=1-k -2 4;2 3-k 1;1 1 1-k;det(A) ans = -6*k+5*k2-k3 factor(det(A) ans = -k*(k-2)*(-3+k)从而得到：当k = 0、k = 2、k = 3时，原方程组有非0解。7.2.3.2 符号矩阵的展开expand (s) 符号表达式展开函数，s为符号矩阵或符号表达式，常用于多项式的因式分解中，也常用于三角函数，指数函数和对数函数的展开中。例7-3：将(x+1)3，sin (x+y
6、)展开。syms x yp=expand(x+1)3) p = x3+3*x2+3*x+1 q=expand(sin(x+y) q = sin(x)*cos(y)+cos(x)*sin(y)7.2.3.3 同类项合并collect (s, v) 将s中的变量v的同幂项系数合并。collect (s) s为矩阵或表达式，此命令对由命令findsym函数返回的默认变量进行同类项合并。7.2.3.4 符号简化simple (s) s为矩阵或表达式, 寻找符号矩阵或符号表达式的最简型。说明：simple (s)将表达式的长度化到最短。若还想让表达式更加精美，可使用函数pretty。pretty (s)
7、 使表达式更加精美。例7-4：计算行列式 的值。 syms a b c dA=1 1 1 1;a b c d;a2 b2 c2 d2;a4 b4 c4 d4;d1=det(A) d1 = b*c2*d4-b*d2*c4-b2*c*d4+b2*d*c4+b4*c*d2-b4*d*c2-a*c2*d4+a*d2*c4+a*b2*d4-a*b2*c4-a*b4*d2+a*b4*c2+a2*c*d4-a2*d*c4-a2*b*d4+a2*b*c4+a2*b4*d-a2*b4*c-a4*c*d2+a4*d*c2+a4*b*d2-a4*b*c2-a4*b2*d+a4*b2*c d2=simple(d1)
8、d2 = -(d-c)*(b-c)*(b-d)*(-c+a)*(a-d)*(a-b)*(a+d+c+b) pretty(d2) -(d - c) (b - c) (b - d) (-c + a) (a - d) (a - b) (a + d + c + b)7.3 秩与线性相关性7.3.1 矩阵和向量组的秩以及向量组的线性相关性矩阵的秩是矩阵A中最高阶非零子式的阶数，向量组的秩通常由该向量组构成的矩阵来计算。rank (A) A为矩阵，rank为矩阵求秩函数。例7-5：求向量组a1= (1,-2,2,3), a2 = (-2,4,-1,3), a3 = (-1,2,0,3), a4 = (0,
9、6,2,3)的秩，并判断其线性相关性。A=1 -2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;rank(A)ans = 3由于秩为3 向量的个数，因此向量组线性相关。7.3.2 向量组的最大无关组矩阵可以通过初等行变换化成行最简形，从而找出列向量组的最大无关组，Matlab将矩阵化成行最简形的命令是：rref。 格式：rref (A) A为矩阵例7-6：求向量组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)的一个最大无关组。a1=1,-2,2,3;a2 =-
10、2,4,-1,3;a3 =-1,2,0,3;a4 =0,6,2,3;a5 =2,-6,3,4;A=a1 a2 a3 a4 a5A = 1 -2 -1 0 2 -2 4 2 6 -6 2 -1 0 2 3 3 3 3 3 4format rat %以有理格式输出B=rref(A) %求A的行最简形B = 1 0 1/3 0 16/9 0 1 2/3 0 -1/9 0 0 0 1 -1/3 0 0 0 0 0 从B中可以得到：向量组a1, a2, a4是其中一个最大无关组。7.4 线性方程组的求解7.4.1 求线性方程组的唯一解或特解这类问题的求法分为两类：一类主要用于解低阶稠密矩阵直接法；另一类
11、是解大型稀疏矩阵迭代法。7.4.1.1 利用矩阵除法求线性方程组的特解（或一个解）方程：AX=b解法：X=Ab说明：如果A为平方矩阵，则除了算法不同，Ab与inv (A)*b大致相当。如果A是一个nn的矩阵，b是一个具有n个元素的列向量或具有多个此类列向量的矩阵，则X=Ab为用Gauss消元法得到的方程AX=b的解。如果A奇异，则会给出警告信息。如果A为mn的矩阵（mn），且b为具有m个元素的列向量或具有多个此类列向量的矩阵，则X=Ab为等式系统AX=b的最小二乘意义上的解。例7-7：求方程组的解 A=5 6 0 0 0;1 5 6 0 0;0 1 5 6 0;0 0 1 5 6;0 0 0
12、1 5;b=1 0 0 0 1;R_A=rank(A)R_A = 5X=AbX = 2.2662 -1.7218 1.0571 -0.5940 0.3188这就是方程组的解。例7-8：求方程组的一个特解 A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;X=AbWarning: Rank deficient, rank = 2 tol = 8.8373e-015.X = 0 0 -8/15 3/5 format shortX = 0 0 -0.5333 0.60007.4.1.2 利用矩阵的LU、QR和cholesky分解求解方程组的解1. LU分解LU分解又称为
13、Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行变换）和上三角矩阵的乘积。即A = LU，L为下三角矩阵的基本变换，U为上三角阵。则：A*X = b 变成L*U*X = b 所以 X = U (L b) 这样可以大大提高运算速度。命令 L, U = lu(A)例7-9：求解方程组的一个特解 A=4 2 -1;3 -1 2;11 3 0;b=2 10 8;D=det(A)D = 0L,U=lu(A)L = 0.3636 -0.5000 1.0000 0.2727 1.0000 0 1.0000 0 0U = 11.0000 3.0000 0 0 -1.8182 2.0000 0
14、 0 0.0000X=U(Lb)Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.018587e-017.X = 1.0e+016 * -0.4053 1.4862 1.3511说明：结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。2. Cholesky分解若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：A = R*R 其中R为上三角矩阵。方程 A*X = b 变成 R*R*X = b所以 X = R (
15、Rb)命令 R = chol (A)3. QR分解对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A = QR方程 A*X = b 变成 QRX = b所以 X = R (Qb)命令 Q, R = qr (A)例7-10：求方程组的一个特解 A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;Q,R=qr(A)Q = -0.3015 0.1421 -0.9428 -0.9045 -0.3553 0.2357 -0.3015 0.9239 0.2357R = -3.3166 -0.9045 6.3317 -0.9045 0
16、5.1168 -7.6752 -8.9544 0 0 -0.0000 -0.0000X=R(Qb)Warning: Rank deficient, rank = 2 tol = 8.8373e-015.X = 0 0 -0.5333 0.6000说明：这三种分解，在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。7.4.2 求线性齐次方程组的通解在Matlab中，函数null用来求解零空间，即满足A *X = 0的解空间，实际上是求出解空间的一组基（基础解系）。格式：z = null (A) %z的列向量为方程组的正交规范基，满足于ZZ = E z = null(A,
17、r) %z的列向量为方程组A X = 0的有理基。例7-11：求解方程组的通解 A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3;format ratB=null(A,r)B = 2 5/3 -2 -4/3 1 0 0 1 写出通解：syms k1 k2x=k1*B(:,1)+k2*B(:,2) %写出方程组的通解x = 2*k1+5/3*k2 -2*k1-4/3*k2 k1 k2 pretty(x) %让通解表达式更加精美 2 k1 + 5/3 k2 -2 k1 - 4/3 k2 k1 k2 7.4.3 求非齐次线性方程组的通解非齐次线性方程组需要先判断方程组是否有解，若有解，再
18、去求通解。因此步骤为：第一步：判断AX = b是否有解，若有解则进行第二步；第二步：求AX = b的一个特解；第三步：求AX = 0的通解；第四步：AX = b的通解 = AX = 0的通解+AX = b的一个特解。例7-12：求解方程组 在Matlab编辑器中建立M文件：LX0601.mA=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(
19、A,r) %求AX=0的基础解系else X=equation no solve %判断无解end运行后结果显示： R_A = 2 R_B = 3 X =equition no solve说明该方程组无解。例7-13：求解方程组的通解 解一：在Matlab编辑器中建立M文件：LX0602.mA=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;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_A In E:matlab6p1workLX0602.m at li
20、ne 11X = 0 0 -8/15 3/5 C = 3/2 -3/4 3/2 7/4 1 0 0 1 所以原方程组的通解为syms k1 k2X=k1*C(:,1)+k2*C(:,2)+XX = 3/2*k1-3/4*k2 3/2*k1+7/4*k2 k1-8/15 k2+3/5 pretty(X) 3/2 k1 - 3/4 k2 3/2 k1 + 7/4 k2 k1 - 8/15 k2 + 3/5 解二：在Matlab编辑器中建立M文件：LX0603.mA=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8;b=1 4 0;B=A b;C=rref(B) %求增广矩阵的行最简形，可
21、得最简同解方程组。运行结果显示为：C = 1 0 -3/2 3/4 5/4 0 1 -3/2 -7/4 -1/4 0 0 0 0 0 对应齐次方程组的基础解系为： 非齐次方程组的特解为： 所以原方程组的通解为： 7.5 特征值与二次型工程技术中的一些问题，如振动问题和稳定性问题，常归结为求一个方阵的特征值和特征向量。7.5.1 方阵的特征值与特征向量在Matlab中，用如下几种调用格式来求A的特征值和特征向量。1. d=eig(A) % d为矩阵A的特征值排成的向量。2. V, D=eig(A) % D为A的特征值对角阵，V的列向量为对应特征值的特征向量（且为单位向量）3. V, D=eig(
22、A, nobalance) % 当A中有小到和截断误差相当的元素时，用nobalance选项，其作用是减少计算误差。例7-14：A=1 2 3;4 5 6;7 8 9;x,y=eig(A)x = -0.2320 -0.7858 0.4082 -0.5253 -0.0868 -0.8165 -0.8187 0.6123 0.4082y = 16.1168 0 0 0 -1.1168 0 0 0 -0.0000其中x为特征向量矩阵，y为特征值矩阵。7.5.2 正交矩阵及二次型7.5.2.1 向量的长度（范数）命令 norm格式 norm (X) %求X的范数7.5.2.2 求矩阵的正交矩阵命令 o
23、rth格式 orth (A) %将矩阵A正交规范化例7-15：A=4 0 0;0 3 1;0 1 3;P=orth(A)P = 0 1 0 -985/1393 0 -985/1393 -985/1393 0 985/1393 format shortPP = 0 1.0000 0 -0.7071 0 -0.7071 -0.7071 0 0.7071Q=P*PQ = 1.0000 0 0.0000 0 1.0000 0 0.0000 0 1.00007.5.2.3 矩阵的schur分解格式 U, T= schur (A) %U为正交矩阵，使得A = UTU和UU=ET = schur (A) %
24、生成schur矩阵T。当A为实对称矩阵时，T为特征值对角阵。例7-16：A=4 0 0;0 3 1;0 1 3;U,T=schur(A)U = 0 0 1.0000 -0.7071 0.7071 0 0.7071 0.7071 0T = 2 0 0 0 4 0 0 0 4这里U就是所求的正交矩阵P，T就是对角矩阵。V,D=eig(A)V = 0 0 1.0000 -0.7071 0.7071 0 0.7071 0.7071 0D = 2 0 0 0 4 0 0 0 4这里V就是所求的正交矩阵P，D就是对角矩阵。说明：对于实对称矩阵，用eig和schur分解效果一样。例7-17：用一个正交变换X
25、 = PY，把二次型 f = 2x1x2+2x1x3-2x1x4-2x2x3+2x2x4+2x3x4化成标准形。解：先写出二次型的实对称矩阵 在Matlab编辑器中建立M文件：LX0603.mA=0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0;P,D=schur(A)syms y1 y2 y3 y4y=y1;y2;y3;y4; %这里不用行向量的转置表示，以免出现复数。X=vpa(P,2)*y %vpa表示可变精度计算，这里取2位精度。f=y1 y2 y3 y4*D*yP = -1/2 390/1351 780/989 780/3691 1/2 -390/1351 780/3691 780/989 1/2 -390/1351 780/1351 -780/1351 -1/2 -1170/1351 0 0 D = -3 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 X = -.50*y1+.29*y2+.79*y3+.21*y4 .50*y1-.29*y2+.21*y3+.79*y4 .50*y1-.29*y2+.56*y3-.56*y4 -.50*y1-.85*y2 f = -3*y12+y22+y32+y42-