6解线性方程组

6.4线性方程组求解6.4线性方程组求解科学运算中最重要问题之一就是求非齐次解线性方程组Ax=b的解•根据系数A矩阵的秩和增广矩阵B＝[A,b]的秩以及未知数个数n关系确定解情况，•1）系数矩阵的秩＝增广矩阵的秩＝n，方程组有唯一解•2）系数矩阵的秩＝增广矩阵的秩＜n，方程组有无穷多解•3）系数矩阵的秩＜增广矩阵的秩，方程组无解6.4.1直接解法1．利用左除运算符的直接解法对于线性方程组Ax=b，可以利用左除运算符“\”求解：x=A\b123123123348379122579xxxxxxxxx系数矩阵：A=134379257常数项矩阵：b＝8129123123123348379122579xxxxxxxxxA=[134;379;257];b=[8129]';B=[A,b];n=3;RA=rank(A);RB=rank(B);ifRA==RB&RA==nX=A\belseifRA==RB&RAnfprintf('方程组有一特解')X=A\belsefprintf('方程组无解')endendA=[134;379;257];b=[8129]';X=A\b对于不滿秩情况，左除给出方程特解，以下列方程组为例：1234123412345980334431xxxxxxxxxxxxA=[15-9-8;3-1-34;11-3-1];b=[041]’;B=[A,b];n=3;RA=rank(A)RB=rank(B)x=A\bsolve()也可以解线性方程,对于滿秩解与左除相同，对于不滿秩情况给出有参数的解，而不象左除给出特解：solve(‘eq1’,‘eq2’,…..’eqn’):对n个方程的默认变量求解solve(‘eq1’,‘eq2’,…..’eqn’，‘var1,var2,…varn’,):对n个方程的var1,varn变量求解S=solve(‘eq1’,‘eq2’,…..’eqn’，‘var1’,‘var2’,…‘varn’,)[x1,x2,..xn]=solve(‘eq1’,‘eq2’,…..’eqn’，‘var1’,‘var2’,…‘varn’,)例：[x,y,z]=solve('x+2*y+3*z-1=0','2*x+y+2*z-1=0','x+3*y+4*z-1=0','x','y','z')[x,y,z,w]=solve('x+5*y-9*z-8*w=0','3*x-y-3*z+4*w-4=0','x+y-3*z-w-1=0','x,y,z,w')2．利用矩阵的分解求解线性方程组矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。(1)LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：[L,U]=lu(X)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。[L,U,P]=lu(X)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组Ax=b的解x=U\(L\b)或x=U\(L\P*b)，这样可以大大提高运算速度。例6-2用LU分解求解线性方程组。2X1+x2-5x3+x4=13X1-5x2+7x4=13X2+x3-x4=13X1+6x2-x3-4x4=13命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2种格式，命令如下：[L,U,P]=lu(A);x=U\(L\P*b)(2)QR分解对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：[Q,R]=qr(X)：产生一个一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。[Q,R,E]=qr(X)：产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解x=R\(Q\b)或x=E*(R\(Q\b))。例6-3用QR分解求解6-2线性方程组。命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[Q,R]=qr(A);x=R\(Q\b)或采用QR分解的第2种格式，命令如下：[Q,R,E]=qr(A);x=E*(R\(Q\b))(3)Cholesky分解如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R，则下三角矩阵为其转置，即X=R'R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为：R=chol(X)：产生一个上三角阵R，使R'R=X。若X为非对称正定，则输出一个出错信息。[R,p]=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足R'R=X(1:q,1:q)。实现Cholesky分解后，线性方程组Ax=b变成R‘Rx=b，所以x=R\(R’\b)。例6--4用Cholesky分解求解例6-2中的线性方程组。命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';R=chol(A)???Errorusing==cholMatrixmustbepositivedefinite命令执行时，出现错误信息，说明A为非正定矩阵。6.4.2迭代解法迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。1．Jacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，即aii≠0(i=1,2,…,n)，则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵和上三角阵，于是Ax=b化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。Jacobi迭代法的MATLAB函数文件Jacobi.m如下：function[y,n]=jacobi(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin3errorreturnendD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D\(L+U);f=D\b;y=B*x0+f;n=1;%迭代次数whilenorm(y-x0)=epsx0=y;y=B*x0+f;n=n+1;end例6-5用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在命令中调用函数文件Jacobi.m，命令如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)2．Gauss-Serdel迭代法在Jacobi迭代过程中，计算时，已经得到，不必再用，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function[y,n]=gauseidel(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin3errorreturnendD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;%迭代次数whilenorm(y-x0)=epsx0=y;y=G*x0+f;n=n+1;end例6-6用Gauss-Serdel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在命令中调用函数文件gauseidel.m，命令如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)例6-7分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。命令如下：a=[1,2,-2;1,1,1;2,2,1];b=[9;7;6];[x,n]=jacobi(a,b,[0;0;0])[x,n]=gauseidel(a,b,[0;0;0])6.4.3超定系统求解（应用方程个数大于自变量个数。）2X1-x2+3x3=33X1+x2-5x3=04X1-x2+x3=3X1+3x2-13x3=-6A=[2-13;31-5;4-11;13-13]b=[303-6]’;Rank(A)X1=A\b验证：A*x1-b超定方程所得解不一定满足Ax=b，x只是最小二乘意义上的解有一组实验数据：ty00.30.81.11.62.20.820.720.630.600.550.50被认为有如下规律：12()tytxex代入数据可得两个未知数的6个方程，用左除方法可以求得x1,x2t=[00.30.81.11.62.2]';y=[0.820.720.630.600.550.50]';e=[ones(size(t))exp(-t)];x=e\y;t1=[0:0.1:2.5]';y1=[ones(size(t1)),exp(-t1)]*xplot(t1,y1,t,y,'ro')00.511.522.50.50.550.60.650.70.750.80.850.96.5非线性方程数值求解6.5.1单变量非线性方程求解在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。（一元函数零点求根）该函数的调用格式为：x=fzero（’fun’，x0）x=fzero(‘fun’,x0,tol,trace)其中fname是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0。例6-8求f(x)=x-10x+2=0在x0=0.5附近的根。步骤如下：(1)建立函数文件funx.m。functionfx=funx(x)fx=x-10.^x+2;(2)调用fzero函数求根。x=fzero('funx',0.5)x=0.37586.5.2非线性方程组的求解对于非线性方程组F(X)=0，用fsolve函数求其数值解。fsolve函数的调用格式为：fsolve(fun,x0),求fun=0在估计值x0附近的解（所用的方法是牛顿法），注意这里的fun是可用由inline()或者由function建立的方程函数，也可用单引