线形代数方程组的解法

第五章线形代数方程组的解法 学习目标了解求解线性代数方程组的基本定理掌握高斯消元法及其计算机实现理解用矩阵分解法求解线性方程组的原理掌握大型稀疏方程组的迭代法 5.1线形代数方程组的直接解法n阶线形方程组:(5.1)矩阵形式:其中要求A非奇异.解的判断数值求解方法有以下三条途径（三种框架）直接法：利用Gauss消元或矩阵分解，通过有限次运算可求出精确解。迭代法：构造迭代格式，产生迭代序列，通过无限次迭代过程求解。5.1.1高斯消去法高斯消元法是一个古老的直接法,由它改进得到的选主元法,是目前计算机上常用于求低阶稠密矩阵方程组的有效方法,其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题。高斯消元法的求解过程,可大致分为两个阶段:首先,把原方程组化为上三角形方程组,称之为“消元”过程;然后,用逆次序逐一求出上三角方程组(原方程组的等价方程组)的解,称之为“回代”过程. 5.1.1高斯消去法思想:把矩阵A化为一个上三角矩阵，从而将原方程组约化为容易求解的等价三角方程组，再通过回代过程既可逐一求出各未知数.方法:逐列消元.例:矩阵形式:增广矩阵:对A的第一列消元,进行行变换把A21,A31化为零:对A的第二列消元,进行行变换把A32化为零:进一步:计算x,对于n阶方程组,增广矩阵为:对第一列消元:重新记为:对第二列消元得到:继续消元过程,最后得到:计算x的值:或者写成:注1:在消元过程中所有不能为零.MATLAB程序:functionx=gauss1(A,B)%gauss消去法解线性方程组Ax=Bn=length(B);C=[AB];fori=1:nifC(i,i)~=0C(i,:)=C(i,:)/C(i,i);forj=i+1:nC(j,:)=C(j,:)-C(j,i)*C(i,:);endelseerror('对角线元素为零');endendx(n)=C(n,n+1);fori=n-1:(-1):1ss=C(i,i+1:end-1)*x(i+1:end)';x(i)=C(i,n+1)-ss;end>>A=[231;111;1-2-1];>>B=[9;4;-4];>>x=gauss1(A,B)x=121>>A=[031;111;1-2-1];>>x=gauss1(A,B)???Errorusing==>gauss1对角线元素为零注2：在消元过程中，若相对于该列中对角线以下的元素相比，其绝对值很小时，尽管消去运算可以进行下去，但是用其作除数，即使很小的舍入误差也会引起计算结果的严重扩散和失真。例：其准确到小数点后第9位的解为但是若按照高斯消去法且在编程计算过程中采用4位浮点数求解，用第一个方程消去第二个方程的，得回代以后解得与我们的精确结果相比较，结果严重失真。其主要原因是除数太小，导致舍入误差增大，有效数字失效。若我们在消元前交换两个方程的位置，变为：对上述方程组消元得到三角方程组回代得解 5.1.2高斯列主元消去法优点: 解决直接消去法中的问题. 提高计算精度.方法:消元过程中取该列中(尚未消元的行)绝对值最大的作为主元进行消元.例:增广矩阵:对第一列找主元‘381’进行行变换:消元:对第二列选主元,消元计算x 在消元过程中不会出现待消的某列全为零的情况,否则A奇异.如:方程:得:>>A=[1e-162;11];>>B=[2;3];>>x=gauss1(A,B)x=4.00001.0000>>A\Bans=2.00001.0000高斯列主元消去法程序:functionx=gauss2(A,B)%gauss列主元消去法解线性方程组Ax=Bn=length(B);C=[A,B];fori=1:nc=abs(C(i:end,i));k=find((c-max(c))==0);k=k(1)+i-1;Ci=C(i,:);C(i,:)=C(k,:);C(k,:)=Ci;ifC(i,i)~=0C(i,:)=C(i,:)/C(i,i);forj=i+1:nC(j,:)=C(j,:)-C(j,i)*C(i,:);endelseerror('矩阵A奇异');endendx(n)=C(n,n+1);fori=n-1:(-1):1ss=C(i,i+1:end-1)*x(i+1:end)';x(i)=C(i,n+1)-ss;endfind()函数:找出非零元素.k=find(x)返回向量x的非零元素下标.[i,j]=find(X)返回矩阵X的非零元素下标.[i,j,v]=find(X)返回矩阵X的非零元素下标,同时返回非零元素的值.>>x=[11033055];>>k=find(x)k=135>>y=x(find(x))y=113355>>find(x==0)ans=24>>x==0ans=01010>>find(x>20)ans=35>>M=[123;040]M=123040>>[i,j]=find(M==0)i=22j=13>>[i,j,v]=find(M)i=j=v=111122224133>>A=[031;111;1-2-1];>>B=[9;4;-4];>>x=gauss1(A,B)???Errorusing==>gauss1对角线元素为零>>x=gauss2(A,B)x=1.66673.3333-1.0000>>A=[2-10;111;1-2-1]A=2-101111-2-1>>x=gauss2(A,B)???Errorusing==>gauss2矩阵A奇异>>x=gauss1(A,B)???Errorusing==>gauss1对角线元素为零>>A=[1e-162;11];>>B=[2;3];>>x=gauss2(A,B)x=2.00001.0000 5.2线形方程组的迭代解法线性方程组的迭代解法就是根据所给的方程组AX=b， 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 出一个迭代公式，然后将任意选取的初始向量带入迭代公式，求出，再将同一代入迭代公式，求出，如此反复进行，得到向量序列。当收敛时，其极限即为原方程组的解。本节我们介绍一下雅克比（Jacobi）迭代法和高斯－赛德尔（Gauss-Seidel）迭代法。 5.2.1简单迭代法(雅可比迭代法)线形方程组:矩阵形式:如果矩阵A的对角线元素不为零,即Aii0,则可把方程变为:雅可比迭代法:给定的一组x初始值由以下公式用求如果满足中止条件则停止,否则k=k+1,重复2.即:把此公式写为矩阵形式:令:则有:故雅可比迭代公式的矩阵形式为:令:迭代公式:MATLAB程序:function[x,k]=jacobi(A,B)%Jacobi迭代法解方程组Ax=Bv=diag(A);ifall(v)D=diag(v);Di=diag(1./v);M=-Di*(A-D);d=Di*B;x1=v;x2=x1+1;k=0;whilemax(abs(x2-x1))>0.0001ifk>300disp('可能不收敛')break;endx1=x2;x2=M*x1+d;k=k+1;endelseerror('A对角线有零元素');endall(),全不为零则返回1,否则返回0.>>all([123])ans=1>>all([103])ans=0any(),全为零返回0,否则返回1.>>any([000])ans=0>>any([023])ans=1>>A=[10-1-2;1-102;11-5];>>B=[2;-13;-12];>>[x,k]=jacobi(A,B)x=1.00002.00003.0000k=11>>A=[231;111;1-2-1];>>B=[9;4;-4];>>[x,k]=jacobi(A,B)可能不收敛x=1.0e+077*1.0504-2.31283.1358k=301 赛得尔迭代法方程组:转化为:赛得尔迭代法:给定的一组x初始值由以下公式用求由求.依次由求3.如果满足中止条件则停止,否则k=k+1,重复2.即:把此公式写为矩阵形式:令:方程:即:迭代公式:function[x,k]=seidel(A,B)%Seidel迭代法解方程组Ax=Bv=diag(A);n=length(v);ifall(v)DL=tril(A);U=triu(A,1);C=[DLeye(n)];DLi=C(:,n+1:end);M=-DLi*U;d=DLi*B;x1=v;x2=x1+1;k=0;whilemax(abs(x2-x1))>0.0001ifk>300disp('可能不收敛')break;endx1=x2;x2=M*x1+d;k=k+1;endelseerror('A对角线有零元素');end>>A=[10-1-2;1-102;11-5];>>B=[2;-13;-12];>>[x,k]=seidel(A,B)x=1.00002.00003.0000k=7 5.3几个要注意的问题 5.3.1奇异方程组矩阵除法:X=A\B矩阵方程AX=B的解X=A/B矩阵方程XA=B的解A为mn阶矩阵,m=n,求精确解.m>n,求最小二乘解.m<n,求基本解.>>A=randn(3);>>B=randn(3,1);>>A\Bans=-0.35830.6884-0.1500>>A=randn(3,2);>>A\Bans=-0.46970.7990>>A=randn(2,3);>>B=randn(2,1);>>A\Bans=0-1.2119-0.4689>>A=[12;1-3];>>B=[13;3];>>A\Bans=92>>A=[12;24];>>B=[4;3];>>A\BWarning:Matrixissingulartoworkingprecision.(Type"warningoffMATLAB:singularMatrix"tosuppressthiswarning.)ans=InfInf>>B=[4;8];>>A\BWarning:Matrixissingulartoworkingprecision.(Type"warningoffMATLAB:singularMatrix"tosuppressthiswarning.)ans=InfInf>>[x,k]=seidel(A,B)x=-65k=2 5.3.2病态方程组方程组系数的微小变化导致方程组解的巨大变化的方程组称为病态方程组.>>A1=[110;10101];>>B=[11;111];>>x1=gauss2(A1,B)x1=11>>A2=[110;10.1100];>>x2=gauss2(A2,B)x2=10.00000.1000可用矩阵A的条件数度量.cond()矩阵条件数函数,当矩阵条件数接近1时,说明矩阵非病态.>>cond(A1)ans=1.0402e+004>>cond(A2)ans=1.0203e+004>>A=[12;1-3];>>cond(A)ans=2.6180 5.3.3迭代法的收敛性求线形方程组的迭代法公式:对任意初值x0和任意d都收敛的充分必要条件: 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示M的特征值,称为M的谱半径.eig()求矩阵特征值函数.det()求矩阵行列式函数.inv()求矩阵逆函数.>>A=[10-1-2;1-102;11-5];>>v=diag(A);>>D=diag(v);>>L=tril(A,-1)L=000100110>>U=triu(A,1)U=0-1-2002000>>M1=-inv(D)*(L+U);>>eig(M1)ans=-0.2372-0.10000.3372>>M2=-inv(D+L)*U;>>eig(M2)ans=0-0.03180.1258>>A=[231;111;1-2-1];>>v=diag(A);>>D=diag(v);>>L=tril(A,-1);>>U=triu(A,1);>>M1=-inv(D)*(L+U);>>eig(M1)ans=1.8100-0.1683-1.6418 思考:用迭代法求方程: