第8讲方程及方程组解法

第八讲方程及方程组解法2本讲教学目标掌握线性方程解法掌握线性方程组的直接解法掌握线性方程组的迭代解法了解非线性方程和方程组的函数解法了解非线性方程和方程组的迭代解法3线性方程组的求解不仅在工程技术领域涉及，而且在其他许多领域也经常碰到，因此它的应用相当广泛。线性方程组的解法一般分为两类：一类是直接法，另一类是迭代法。非线性方程组的解法主要采用迭代法求解，常用的有不动点迭代法、牛顿迭代法和Broyden迭代法等几种。48.1线性方程及方程组的解法8.1.1线性方程的解法通过调用函数roots求解。例1求解方程的根。p=[1,-6,-72,-27]r=roots(p)p=1-6-72-27r=12.1229-5.7345-0.38843267227xxx58.1.2线性方程组的解法1.直接解法(1)利用左除和右除运算符。通过这两个运算符，程序会自动根据输入的系数矩阵判断选用哪种方法进行求解。对线性方程组Ax=b，利用左除运算符“\”求解：x=A\b右除同样。6例2：求解以下方程组。a=[12;23];b=[8;13];x=a\bx=2.00003.00001212282313xxxx7(2)利用矩阵分解方法①LU分解将非奇异的方阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=L*U。Ax=b可变换为x=U\(L\b)或x=U\(L\Pb)分解格式：[L,U]=lu(A)——产生三角阵U和下三角阵L。[L,U,P]=lu(A)——产生三角阵U和下三角阵L以及一个置换矩阵P。8②QR分解：把矩阵A分解为一个正交矩阵Q和一个上三角矩阵R的乘积，即A=Q*R。Ax=b可变换为x=R\(Q\b)或x=E(R\(Q\b))分解格式：[Q,R]=qr(A)——产生正交阵Q和上三角阵R。[Q,R,E]=qr(A)——产生正交阵Q、上三角阵R及置换阵E。9③Cholesky分解：把对称正定的矩阵A分解成一个上三角矩阵R和其转置矩阵R’的乘积，即A=R‘R。Ax=b可变成x=R\(R’\b)分解格式：R=chol(A)——产生一个上三角阵R。[R,p]=chol(A)——A对称正定时，p=0；否则p为一个正整数。如果A满秩，则R为一个阶为p-1的上三角阵，且满足R‘R=A(1:p-1,1:p-1)。10例3：求解方程组。a=[12-33;-163-1;111];b=[15;-13;6];[L,U]=lu(a);x=U\(L\b)x=1.00002.00003.0000123123123123315163136xxxxxxxxx112.迭代解法迭代解法适合求解大型系数矩阵的方程组。迭代解法主要包括：Jacobi迭代法Gauss-Serdel迭代法超松弛迭代法两步迭代法12(1)Jacobi迭代法对线性方程组Ax=b，若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如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。13function[y,n]=jacobi(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin3errorreturnendD=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)=epsx0=y;y=B*x0+f;n=n+1;endJacobi迭代法的Matlab函数文件Jacobi.m如下：14例4：用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)121232310910272106xxxxxxx15(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迭代用新分量代替旧分量，精度会高些。16function[y,n]=gauseidel(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin3errorreturnendD=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)=epsx0=y;y=G*x0+f;n=n+1;endGauss-Serdel迭代法的Matlab函数文件gauseidel.m如下：17例5：用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)121232310910272106xxxxxxx188.2非线性方程及方程组的解法8.2.1调用函数求解1.使用fzero函数求解，其格式为：z=fzero('fname',x0,tol,trace)其中fname是函数文件名，x0为搜索起点。一个函数可能有多个根，但只给出离x0最近的根。tol控制结果的相对精度，缺省取eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取0。19例6：求f(x)=x-10x+2=0在0.5附近的根。步骤如下：x0=0.5;%初值z=fzero('x-10^x+2',x0)z=0.3758x1=[01];%定义包含根的区间z=fzero('x-10^x+2',x1)z=0.3758函数调用式如下：①建立函数文件funx.m。functionfx=funx(x)fx=x-10^x+2;②调用fzero函数求根。z=fzero('funx',0.5)z=0.3758202.使用fsolve函数求解调用格式：x=fsolve('fun',x0,option)其中fun是函数文件名，x0是初值，option为最优化工具箱的选项设定，默认缺省。最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来。如果想改变其中某个选项，则可以调用optimset()函数来完成。21例7：求解方程在[0,10]内的解。首先作图如下：fplot(‘[5*x^2*sin(x)-exp(-x),0]’,[0,10])发现在[0,10]区间中有4个解，分别在0,3,6,9附近，所以用命令：fun=inline(‘5*x.^2.*sin(x)-exp(-x)’);fsolve(fun,[0369],1e-6)ans=0.50183.14076.28329.424825sin0xxxe22例8：求非线性方程组在(0.5,0.5)附近的解。首先编制函数文件fc.m如下：functiony=fc(x)y=[x(1)-0.7*sin(x(1))-0.2*cos(x(2));x(2)-0.7*cos(x(1))+0.2*sin(x(2))];然后在命令窗口输入：x0=[0.50.5];fsolve(‘fc’,x0)或fsolve(@fc,x0)1122120.7sin0.2cos00.7cos0.2sin0xxxxxx结果如下：ans=0.52650.5079238.2.2其他解法非线性方程(1)对分法(2)迭代法非线性方程组(1)不动点迭代法(2)牛顿迭代法(3)Broyden迭代法具体使用格式和方法参考教材P86-94。24例9编制程序，对任何函数求根，在运行该程序时将显示函数与零点的图形，用户在可能的零点附近用鼠标点击采点，程序通过循环计算每个采集点的根，并输出结果。以函数求零点为例，综合叙述相关命令的用法。tbettfat)(sin)(2-10-50510-5-4-3-2-101ty(t)25y=inline('sin(t)^2*exp(-a*t)-b*abs(t)','t','a','b');a=0.1;b=0.5;t=-10:0.01:10;%对自变量采样y_char=vectorize(y);Y=feval(y_char,t,a,b);clf%作人机对话界面plot(t,Y,'r');holdon,plot(t,zeros(size(t)),'k');%画坐标横轴xlabel('t');ylabel('y(t)'),holdofftitle('用鼠标在函数零点附近连续点五点')26%循环使用命令fzero，求鼠标采集点附近的函数根。fori=1:5[t1,y1]=ginput(1);[ttt(i),yyy(i),exitflag]=fzero(y,t1,[],a,b)%在ttt(i)附近搜索0点endA(:,1)=ttt';A(:,2)=yyy';A计算结果为：A=xf(x)-2.00740.0000-2.00740.0000-0.51980.00001.67380.00001.6738-0.000027本节完，谢谢！！