1、MATLAB矩阵分解与线性方程组求解 数学与信息科学系 汪远征,4. MATLAB矩阵分解与线性方程组求解,4.1 矩阵分解 4.2 秩与线性相关性 4.3 线性方程的组的求解 4.4 特征值与二次型,4.1 矩阵分解,4.1.1 LU分解 矩阵的三角分解又称LU分解,它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积,即A=LU。 Matlab使用函数lu实现LU分解,其格式为: L,U = lu(A) 其中U为上三角阵,L为下三角阵或其变换形式,满足LU=A。,4.1 矩阵分解,4.1.1 LU分解 L,U,P = lu(A) U为上三角阵,L为下三角阵,P为单位矩阵的行变
2、换矩阵,满足LU=PA。 例4-1 A=1 2 3;4 5 6;7 8 9; L,U=lu(A) L,U,P=lu(A),4.1 矩阵分解,4.1.2 Cholesky分解 如果A为n阶对称正定矩阵,则存在一个实的非奇异上三角阵R,满足R*R = A,称为Cholesky分解 Matlab使用函数chol实现Cholesky分解,其格式为: R = chol(A) 若A非正定,则产生错误信息。 R,p = chol(A) 不产生任何错误信息,若A为正定阵,则p=0,R与上相同;若A非正定,则p为正整数,R是有序的上三角阵。,4.1 矩阵分解,4.1.2 Cholesky分解 例4-2 A=pa
3、scal(4) %产生4阶pascal矩阵 R,p=chol(A),4.1 矩阵分解,4.1.3 QR分解 将矩阵A分解成一个正交矩阵Q与一个上三角矩阵R的乘积A=QR,称为QR分解。 Matlab使用函数qr实现QR分解,其格式为: Q,R = qr(A) Q,R,E = qr(A) 求得正交矩阵Q和上三角阵R,E为单位矩阵的变换形式,R的对角线元素按大小降序排列,满足AE=QR。 Q,R = qr(A,0) 产生矩阵A的“经济大小”分解,4.1 矩阵分解,4.1.3 QR分解 Matlab使用函数qr实现QR分解,其格式为: Q,R = qr(A) Q,R,E = qr(A) 求得正交矩阵
4、Q和上三角阵R,E为单位矩阵的变换形式,R的对角线元素按大小降序排列,满足AE=QR。 Q,R = qr(A,0) 产生矩阵A的“经济大小”分解 Q,R,E = qr(A,0) E的作用是使得R的对角线元素降序, 且Q*R=A(:, E),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),4.1 矩阵分解,4.1.3 QR分解 返回将矩阵A的第j列移去后的新矩阵的qr分解使用函数 qrdelete,其格式为: Q,R = qrdelete(Q,R,j) 例4-4 A=-149 -50 -154;537
5、180 546;-27 -9 -25; Q,R=qr(A) Q,R=qrdelete(Q,R,3) %将A的第3列去掉后进行qr分解,4.1 矩阵分解,4.1.3 QR分解 在矩阵A中第j列插入向量x后的新矩阵进行qr分解使用函数 qrinsert,其格式为: Q,R = qrinsert(Q,R,j,x) 若j大于A的列数,表示在A的最后插入列x。 例4-5 A=-149 -50 -154;537 180 546;-27 -9 -25; x=35 10 7; Q,R=qr(A) Q,R=qrinsert(Q,R,4,x),4.1 矩阵分解,4.1.4 Schur分解 矩阵的Schur分解:设
6、A Rnn ,则存在正交阵Q使实Schur型其中Rii至多2阶。若1阶,其元素即A的特征值,若2阶其特征值为A的一对共轭复特征值。 Matlab使用函数schur实现Schur分解,其格式为: R = schur(A) R为schur矩阵, 即R的主对角线元素为特征值的三角阵,4.1 矩阵分解,4.1.4 Schur分解 Matlab使用函数schur实现Schur分解,其格式为: R = schur(A) R为schur矩阵, 即R的主对角线元素为特征值的三角阵 R = schur(A,flag) 若A有复特征根,则flag=complex,否则flag=real。 Q,R = schur(
7、A,) 返回正交矩阵Q和schur矩阵R,满足A = Q*R*Q。,4.1 矩阵分解,4.1.4 Schur分解 例4-6 H = -149 -50 -154; 537 180 546; -27 -9 -25 ; Q, R=schur(H),4.1 矩阵分解,4.1.5 实Schur分解转化成复Schur分解 将实舒尔形式转化成复舒尔形式的函数是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),4.2 秩与线性相关性,4.2.1
8、矩阵和向量组的秩与向量组的线性相关性 矩阵A的秩是指矩阵A中最高阶非零子式的阶数,或是矩阵线性无关的行数与列数;向量组的秩通常由该向量组构成的矩阵来计算。 在MATLAB中,求矩阵秩的函数是rank。其格式为: k = rank(A) 返回矩阵A的行(或列)向量中线性无关个数 k = rank(A,tol) tol为给定误差,4.2 秩与线性相关性,4.2.1 矩阵和向量组的秩与向量组的线性相关性 例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)的秩,并判断其线性相关性。 解: A=1
9、-2 2 3;-2 4 -1 3;-1 2 0 3;0 6 2 3;2 -6 3 4; k=rank(A) 结果为由于秩为3 向量个数,因此向量组线性相关。,4.2 秩与线性相关性,4.2.2 求行阶梯矩阵及向量组的基 行阶梯使用初等行变换,矩阵的初等行变换有三条: (1) 交换两行ri rj (第i、第j两行交换) (2) 第i行的k倍kri (3) 第i行的k倍加到第j行上去rj+kri 通过这三条变换可以将矩阵化成行最简形,从而找出列向量组的一个最大无关组。,4.2 秩与线性相关性,4.2.2 求行阶梯矩阵及向量组的基 Matlab将矩阵化成行最简形的命令是rref或rrefmovie。
10、其格式为: R = rref(A) R 是A的行最简行矩阵 R,jb = rref(A) jb是一个向量,其含义为:r = length(jb)为A的秩;A(:, jb)为A的列向量基;jb中元素表示基向量所在的列。,4.2 秩与线性相关性,4.2.2 求行阶梯矩阵及向量组的基 Matlab将矩阵化成行最简形的命令是rref或rrefmovie。其格式为: R,jb = rref(A,tol) tol为指定的精度 rrefmovie(A) 给出每一步化简的过程,4.2 秩与线性相关性,4.2.2 求行阶梯矩阵及向量组的基 例4-9 求向量组a1=(1,-2,2,3),a2=(-2,4,-1,3
11、),a3=(-1,2,0,3),a4=(0,6,2,3),a5=(2,-6,3,4)的一个最大无关组 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 a5 R,jb=rref(A) A(:,jb) 即:a1 a2 a4为向量组的一个基。,4.3 线性方程的组的求解,我们将线性方程的求解分为两类:一类是方程组求唯一解或求特解,另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断: 若系数矩阵的秩r=n(n为方程组中未知变量的个数),则有唯一解; 若系数矩阵的秩rn,则可能有无穷解;
12、线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解;其特解的求法属于解的第一类问题,通解部分属第二类问题。,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 当系数矩阵A为N*N的方阵时,MATLAB会自行用高斯消去法求解线性代数方程组。 若右端项b为N*1的列向量,则x=Ab可获得方程组的数值解x(N*1的列向量); 若右端项b为N*M的矩阵,则x=Ab可同时获得同一系数矩阵A、M个方程组数值解x(为N*M的矩阵),即x(:,j)=Ab(:,j),j=1,2,M。,4.3 线性方程的组的求解,4.3.1 求线性方程
13、组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 解法:X=Ab 例4-10 求方程组 的解。 解: A=5 6 0 0 01 5 6 0 00 1 5 6 0 0 0 1 5 60 0 0 1 5; B=1 0 0 0 1; R_A=rank(A) %求秩 X=AB %求解,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 解: A=5 6 0 0 01 5 6 0 00 1 5 6 0 0 0 1 5 60 0 0 1 5; B=1 0 0 0 1; R_A=rank(A) %求秩 X=AB %求解 运行后结果如下 这就
14、是方程组的解。,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 或用函数rref求解: C=A,B %由系数矩阵和常数列构成增广矩阵C R=rref(C) %将C化成行最简行则R的最后一列元素就是所求之解。,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 %由于系数矩阵不满秩, 该解法可能存在误差 X = 0 0 -0.5333 0.6000
15、 一个特解近似值,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 例4-11 求方程组 的一个特解 解: 若用rref求解,则比较精确: A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; B=1 4 0; C=A,B; %构成增广矩阵 R=rref(C)由此得解向量X=1.2500 0.2500 0 0(一个特解)。,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 例4-12 解方程组 (1)Ax=b1;(2)Ay=b2解法1:分别解方程组 (1)Ax=b
16、1;(2)Ay=b2 A=1 -1 1;5 -4 3;2 1 1; b1=2;-3;1; b2=3;4;-5; x=Ab1 y=Ab2,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线性方程组AX=b的特解 例4-12 解方程组 (1)Ax=b1;(2)Ay=b2解法1:分别解方程组 (1)Ax=b1;(2)Ay=b2 得两个线性代数方程组的解: (1) x1= -3.8, x2= 1.4, x3= 7.2; (2) y1= -3.6, y2= 2.2, y3= 4.4,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 1利用矩阵除法求线
17、性方程组AX=b的特解 例4-12 解方程组 (1)Ax=b1;(2)Ay=b2解法2:将两个方程组连在一起求解:Az = b b=2 3;-3 4;1 -5 z=Ab 很明显,这里的解z的两个列向量便是前面分别求得的两组解x和y,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), 这样可以大大提高运算速度。,4.3 线性方程的组的
18、求解,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(A) X=U(LB),4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 2利用LU、QR和cholesky分解求方程组的解 (1)LU分解: 例4-13 求方程组 的一个特解。 解: 说明 结果中的警告是由于系数行列式为零产生的。可以通过A*X验证其正确性。,4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特
19、解 2利用LU、QR和cholesky分解求方程组的解 (2) Cholesky分解 若A为对称正定矩阵,则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积,即: A = R*R 其中R为上三角阵。 方程A*X=b变成R*R*X = b 所以X=R(Rb),4.3 线性方程的组的求解,4.3.1 求线性方程组的唯一解或特解 2利用LU、QR和cholesky分解求方程组的解 (3)QR分解 对于任何长方矩阵A,都可以进行QR分解,其中Q为正交矩阵,R为上三角矩阵的初等变换形式,即: A = QR 方程A*X=b变形成QRX=b 所以X=R(Qb),4.3 线性方程的组的求解,4.3
20、.1 求线性方程组的唯一解或特解 2利用LU、QR和cholesky分解求方程组的解 (3)QR分解 上例中 Q, R=qr(A) X=R(QB) 说明:这三种分解,在求解大型方程组时很有用。其优点是运算速度快、可以节省磁盘空间、节省内存。,4.3 线性方程的组的求解,4.3.2 求线性齐次方程组的通解 在Matlab中,函数null用来求解零空间,即满足AX=0的解空间,实际上是求出解空间的一组基(基础解系)。其格式为: z = null(A) z的列向量为方程组的正交规范基,满足z z = I。 z = null(A,r) z的列向量是方程AX=0的有理基,4.3 线性方程的组的求解,4.
21、3.2 求线性齐次方程组的通解 例4-14 求方程组 的通解: 解:A=1 2 2 1;2 1 -2 -2;1 -1 -4 -3; format rat %指定有理式格式输出 B=null(A,r) %求解空间的有理基 运行后显示结果如下:,4.3 线性方程的组的求解,4.3.2 求线性齐次方程组的通解 例4-14 求方程组 的通解: 解:或通过行最简型得到基: B=rref(A) 即可写出其基础解系(与上面结果一致)。,4.3 线性方程的组的求解,4.3.2 求线性齐次方程组的通解 例4-14 求方程组 的通解: 解:写出通解: syms k1 k2 X=k1*B(:,1)+k2*B(:,2
22、) %写出方程组的通解 pretty(X) %让通解表达式更加精美 运行后结果如下: 下面是其简化形式,4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 非齐次线性方程组需要先判断方程组是否有解,若有解,再去求通解。 因此,步骤为: 第1步: 判断AX=b是否有解,若有解则进行第二步 第2步: 求AX=b的一个特解 第3步: 求AX=0的通解 第4步: AX=b的通解= AX=0的通解+AX=b的一个特解,4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-15 求解方程组解:在Matlab中建立M文件如下: A=1 -2 3 -1;3 -1 5 -3;2
23、1 2 -2; b=1 2 3; B=A b; n=4; R_A=rank(A), R_B=rank(B), format rat if R_A=R_B&R_A=n %判断有唯一解X=Ab elseif R_A=R_B&R_An %判断有无穷解X=Ab %求特解C=null(A,r) %求AX=0的基础解系 else X=equition no solve %判断无解 end,4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-15 求解方程组解:在Matlab中建立函数文件如下:A=1 -2 3 -1;3 -1 5 -3;2 1 2 -2; b=1 2 3; jie(A,
24、b),4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-15 求解方程组解: 运行后结果显示:说明 该方程组无解,4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-16 求方程组 的通解解法一:调用jie函数 A=1 1 -3 -1;3 -1 -3 4;1 5 -9 -8; b=1 4 0; jie(A,b),4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-16 求方程组 的通解解法一:在Matlab编辑器中建立M文件如下: 运行后结果显示为:所以原方程组的通解为X=k1 +k2 +,4.3 线性方程的组的求解,4.3.3 求
25、非齐次线性方程组的通解 例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) %求增广矩阵的行最简形,得最简同解方程组 运行后结果显示为:,4.3 线性方程的组的求解,4.3.3 求非齐次线性方程组的通解 例4-16 求方程组 的通解解法二:用rref求解 运行后结果显示为: 对应齐次方程组的基础解系为: 非齐次方程组的特解为:所以,原方程组的通解为:X=k11+k22 + *。,4.4 特征值与二次型,工程技术中的一些问题,如振动问题和稳定性问题,常归结为求一个方阵的特征值和特
26、征向量。 4.4.1 特征值与特征向量的求法 设A为n阶方阵,如果数“”和n维列向量x使得关系式Ax = x成立,则称为方阵A的特征值,非零向量x称为A对应于特征值的特征向量。,4.4 特征值与二次型,4.4.1 特征值与特征向量的求法 在MATLAB中,计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有3种: (1) E=eig(A):求矩阵A的全部特征值,构成向量E。 (2) V, D=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的列向量。 (3) V, D=eig(A,nobalance):与第2种格式类似,但第2种格式中先对A作相似变换后求矩
27、阵A的特征值和特征向量,而格式3直接求矩阵A的特征值和特征向量,4.4 特征值与二次型,4.4.1 特征值与特征向量的求法 例4-17 求矩阵 的特征值和特征向量 解:A=-2 1 1;0 2 0;-4 1 3; V,D=eig(A) 结果显示: 即:特征值-1对应特征向量(-0.7071 0 -0.7071)T 特征值2对应特征向量(-0.2425 0 -0.9701)T和(0.3015 0.9045 0.3015)T,4.4 特征值与二次型,4.4.1 特征值与特征向量的求法 例4-18 求矩阵 的特征值和特征向量 解: A=-1 1 0;-4 3 0;1 0 2; V,D=eig(A)
28、结果显示为 说明:当特征值为1 (二重根)时,对应特征向量都是 k (0.4082 0.8165 -0.4082)T,k为任意常数。,4.4 特征值与二次型,4.4.1 特征值与特征向量的求法 例4-19 用求特征值的方法解方程: 3x5 - 7x4 + 5x2 + 2x - 18 = 0 解: p=3,-7,0,5,2,-18; A=compan(p); %A的友阵 x1=eig(A) %求A的特征值 x2=roots(p) %直接求多项式p的零点,4.4 特征值与二次型,4.4.2 复对角矩阵转化为实对角矩阵 函数cdf2rdf将复对角阵d变为实对角阵D,在对角线上,用22实数块代替共轭复
29、数对。其格式 V,D = cdf2rdf (v,d) 例4-15 A=1 2 3;0 4 5;0 -5 4; v,d=eig(A) V,D=cdf2rdf(v,d),4.4 特征值与二次型,4.4.3 正交基 函数orth将矩阵A标准正交化,得到矩阵B。其格式为: B=orth(A) B的列与A的列具有相同的空间,B的列向量是正交向量,且满足: B*B = eye(rank(A),4.4 特征值与二次型,4.4.3 正交基 例4-20 将矩阵 标准正交化。 解: A=4 0 0; 0 3 1; 0 1 3; B=orth(A) Q=B*B 显示结果为,4.4 特征值与二次型,4.4.4 二次型
30、 例4-21 求一个正交变换X=PY,把二次型化成标准形。 解:先写出二次型的实对称矩阵,4.4 特征值与二次型,4.4.4 二次型 例4-21 求一个正交变换X=PY,把二次型化成标准形。 解:在Matlab编辑器中建立M文件如下: A=0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0; P,D=schur(A) syms y1 y2 y3 y4 y=y1;y2;y3;y4; X=vpa(P,2)*y %vpa表示可变精度计算, 这里取2位精度 f=y1 y2 y3 y4*D*y,4.4 特征值与二次型,4.4.4 二次型 例4-21 求一个正交变换X=PY,把二次型化成标准形。 解:运行后结果显示如下:即,