MATLAB线性方程组及矩阵特征值.ppt

第三章 线性代数方程组及矩阵特征值,线性方程组,情形1：m=n 在MATLAB中的求解命令有：,情形2：mn（超定方程），多用于曲线拟合。,解线性方程组的一般函数文件如下：function x,y=line_solution(A,b)m,n=size(A);y=;if norm(b,1)0%非齐次方程组 if rank(A)=rank(A,b)%方程组相容 if rank(A)=m%有唯一解 x=Ab;else%方程组有无穷多个解,基础解系 disp(原方程组有有无穷个解，其齐次方程组的基础 解系为y，特解为x);y=null(A,r);x=Ab;end,else%方程组不相容，给出最小二乘解 disp(方程组的最小二乘法解是：);x=Ab;endelse%齐次方程组 if rank(A)=n%列满秩 x=zero(m,1)%0解 else%非0解 disp(方程组有无穷个解，基础解系为x);x=null(A,r);end endreturn,如在MATLAB命令窗口，输入命令 A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6;x,y=line_solution(A,b)及：A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2;x,y=line_solution(A,b)分别显示其求解结果。,求解线性方程组的主要方法有：,一、Gauss消去法,设 有,线性代数：方法不好时工作量非常大，工作量小的方法是 Gauss 消去法。,3.1 解线性方程组的直接法,消 元：,以此类推，最后方程组化为：,回 代：,二、列主元素消去法-计算结果可靠,到此原方程组化为,到此原方程组化为,(上三角方程组)（3.2),(n-1)原方程组化为,以上为消元过程。,(n)回代求解公式,（3.3)是回代过程。,（3.3),说明：(1)也可采用无回代的列主元消去法(叫Gauss-Jordan消去法)，该法同时消去对角元上 下的元素，且仍旧需要选主元，但比有回 代的列主元消去法的乘除运算次数多。GaussJordan消去法的优点之一是用它来 计算逆矩阵的算法非常容易解释。(2)有回代的列主元消去法所进行的乘除运算 次数为，量很小。,例1：在MATLAB上，用Gauss消去法求解方程组：,程序如下：clear；a=-0.04 0.04 0.12 3;0.56-1.56 0.32 1;-0.24 1.24-0.28 0 x=0,0,0;tempo=a(2,:);a(2,:)=a(1,:);a(1,:)=tempo;aa(2,:)=a(2,:)-a(1,:)*a(2,1)/a(1,1);a(3,:)=a(3,:)-a(1,:)*a(3,1)/a(1,1);a,tempo=a(3,:);a(3,:)=a(2,:);a(2,:)=tempo;aa(3,:)=a(3,:)-a(2,:)*a(3,2)/a(2,2);ax(3)=a(3,4)/a(3,3);x(2)=(a(2,4)-a(2,3)*x(3)/a(2,2);x(1)=(a(1,4)-a(1,2:3)*x(2:3)/a(1,1);x,运行得方程组的解为：,例2：用GaussJordan消去法求解上例中的矩阵 的逆矩阵。,clearA=-0.04 0.04 0.12;0.56-1.56 0.32;-0.24 1.24-0.28a=A,eye(3);tempo=a(2,:);a(2,:)=a(1,:);a(1,:)=tempo;a(1,:)=a(1,:)/a(1,1)for i=2:3;a(i,:)=a(i,:)-a(i,1)*a(1,:);end;atempo=a(3,:);a(3,:)=a(2,:);a(2,:)=tempo;aa(2,:)=a(2,:)/a(2,2);a,for i=1:3 if i=2,a(i,:)=a(i,:)-a(i,2)*a(2,:);endendaa(3,:)=a(3,:)/a(3,3)for i=1:3;if i=3,a(i,:)=a(i,:)-a(i,3)*a(3,:);end;end;aA_inv=a(:,4:6)A*A_inv,三、Gauss 全主元消去法：优点-计算结果更可靠；缺点-挑主元花机时更多，次序有变动，程序复杂。,四、应用（1）求行列式（2）求逆矩阵,(以上过程都应选主元),在MATLAB中用命令det(A),在MATLAB中用命令inv(A)或A(-1)或rref(A,E),记,，则,(下三角上三角）（三角因子分解）,Gauss消元，初等行变换，化原方程组为上三角型。,五矩阵三角分解法,定义3.1,叫,的三角（因子）分解，其中 是,是上三角。,下三角,为单位下三角阵（对角元全为1），,为上三角阵，则称,为Doolittle分解;,若 是下三角，,是单位上三角，则称,定理3.1 n阶阵,有唯一Doolittle分解(Crout）,的前n-1个顺序主子式不为0.,三角分解不唯一,为此引入,定义3.2 若,为Crout分解。,为什么要讨论三角分解？若在消元法进行前能实现三角分解,，则,容易回代求解,回代求解很容易，如,在Gauss消去法中，选主元改变了行的次序，尽管对于Gauss消去法来说，这种次序的变换无法事先知道，但是这个变化的影响却可以用一个算子P表示，其中P是一个置换矩阵。用P左乘原始矩阵A得到：PAx=Py 或 对 做Gauss消去法不需要选主元，所以对 做LU分解，同样不需要选主元。,实际上，将选主元Gauss消去法里的行交换同样作用于单位矩阵，所得矩阵即为P。,在MATLAB中，LU分解的命令是lu，有两种格式：（1）l,u,p=lu(A)其中A是待分解矩阵；l,u,p分别代表L,U,和置换矩阵P 满足：（2）l,u=lu(A)其中，所以,六、解三对角方程组追赶法,给定方程组,按行严格对角占优,（三对角方程组）,其求解算法是Gauss消去法的一种变形，称为三对角法（追赶法）。解法如下：,由第一个方程 得,令，则,将其代入第二个方程，得：,再令，得：,以此类推，一般地，令,得,最后，将 代入第 n个方程 可解得：,以上过程称为“追”；,然后按公式 依顺序,进行回代，求得方程组的解，该过程称为“赶”。,具体算法如下：,三对角方程组的求解程序tri_diag.m如下：,%tri_diag(a,b,c,d,n)solves a tridiagonal equation.function f=tri_diag(a,b,c,d,n)for i=2:n r=a(i)/b(i-1);b(i)=b(i)-r*c(i-1);d(i)=d(i)-r*d(i-1);endd(n)=d(n)/b(n);for i=n-1:-1:1 d(i)=(d(i)-c(i)*d(i+1)/b(i);endf=d;,一.简单迭代法 1.迭代法建立.考虑,(矩阵B不唯一),对应写出,3.2 解线性方程组的迭代法,产生向量序列,若收敛,记,则于(3.4)两端取极限有：,上式说明：,是解向量,从而当k充分大时,注意:迭代阵B不唯一,影响收敛性。,解向量,(3.4)叫简单迭代法,B叫迭代矩阵。,2.收敛性.定义3.3 称,为矩阵B的谱半径。,定理3.4,定理3.3 简单迭代法,设有方程组(其中)Ax=b,即,(3.5),作等价变形,(3.6),-Jacobi迭代法,于是有迭代公式,(k=0,1,2,),(3.7),矩阵形式为：,简记为,(3)设方程组（3.5）的系数矩阵A按行严格对角占优 即：,或按列严格对角占优，即,例3,Jacobi迭代算法,A=9-1-1;-1 10-1;-1-1 15;b=7;8;13;x=0;0;0;er=1;k=0;while er0.00005 er=0;k=k+1;for i=1:3 s=0;t=x(i);for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=t;y(i)=(b(i)-s)/A(i,i);er=max(abs(x(i)-y(i),er);end x=y;xend,0.7778 0.8000 0.86670.9630 0.9644 0.97190.9929 0.9935 0.99520.9987 0.9988 0.99910.9998 0.9998 0.99981.0000 1.0000 1.00001.0000 1.0000 1.0000,二、迭代法 设有简单迭代法 即,(3.8),称如下迭代法,(3.9),为与(3.8)对应的 迭代法，其迭代矩阵 可用“代入法”求得。,与（3.7）对应的迭代公式为：,-Seidel迭代法,（1）迭代法(3.9)对任意 收敛（2）若 则 迭代法(3.9)对 任意 收敛；（3）若简单迭代法(3.8)的迭代矩阵 满 足 或，则相应的Seidel迭代法(3.9)对任意 收敛.,迭代法(3.9)的收敛性,例3,Gauss-Seidel迭代算法,A=9-1-1;-1 10-1;-1-1 15;b=7;8;13;x=0;0;0;er=1;k=0;while er0.00005 er=0;k=k+1;for i=1:3 s=0;t=x(i);for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=(b(i)-s)/A(i,i);er=max(abs(x(i)-t),er);end xend,0.7778 0.8778 0.9770 0.9839 0.9961 0.9987 0.9994 0.9998 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000,松弛因子,=1 即Seidel方法(3.10),(3.11)是一种加权平均。,三.超松弛迭代法(SOR法),SOR方法的收敛性如下（不加证明）：,(1)SOR方法对任意 都收敛的必要条件是：(2)若系数矩阵A对称正定，则 时SOR方法 求解 对任意 收敛；(3)若系数矩阵A按行（或按列）严格对角占优，则 时SOR方法对任意 收敛。,如前例，用SOR法求解,A=9-1-1;-1 10-1;-1-1 15;b=7;8;13;x=0;0;0;er=1;k=0;w=1.1;while er0.00005 er=0;k=k+1;for i=1:3 s=0;t=x(i);for j=1:3 if i=j,s=s+A(i,j)*x(j);end end x(i)=w*(b(i)-s)/A(i,i)+(1-w)*x(i);er=max(abs(x(i)-t),er);end xend,0.8556 0.9741 1.08751.0220 1.0146 0.99390.9988 0.9977 1.00040.9999 1.0003 1.00001.0000 1.0000 1.00001.0000 1.0000 1.0000,3.3 不可解问题,线性方程组并不总是数值可解的，考虑如下三个方程组。,在MATLAB中分别求解如下：,有无穷多个解,其中一个解,方程组不相容,最小二乘解,3.4 病态问题,有许多线性方程组理论上是可解的，但实际计算中由于受到舍入误差的干扰而无法得到精确解，此类问题称为病态问题。,常用范数：,矩阵A的条件数记为Cond(A)，定义为：,条件数总满足：,当矩阵是病态时，其条件数一定很大，但并不能直接说明解的误差。,MATLAB中计算条件数的命令是：cond(A),对于病态矩阵，逆矩阵和行列式的计算都会变得不精确。所以具备下列特征的问题可认为是病态的：,例4：Hilbert矩阵：,计算从55到1414的Hilbert矩阵的条件数、行列式和。,clearfor n=5:14 for i=1:n for j=1:n,a(i,j)=1/(i+j-1);end end c=cond(a);d=det(a)*det(a(-1);fprintf(n=%3.0f cond(a)=%e det*det=%en,n,c,d)end,3.5 矩阵特征值,设A是一个 nn 矩阵，,为A的特征多项式,特征多项式的根称为矩阵A的特征值。,在MATLAB中，特征多项式的系数用命令：c=poly(A)求根：roots(c)也可直接由命令 eig(A)求，同时求出特征向量。,例5：设矩阵,输入矩阵：A=3 4-2;3-1 1;2 0 5计算：eig(A)或 roots(poly(A)若 V,D=eig(A)，则返回特征向量阵V和特征值对角阵D，使AVVD 成立。,作业：3.4，3.8，3.12，3.14，3.21，3.22要求：编写程序，上机运行后保存结果。,