第7章MATLAB解方程与函数极值

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,*,第,7,章,MATLAB,解方程与函数极值,7.1,线性方程组求解,7.2,非线性方程数值求解,7.3,常微分方程初值问题的数值解法,7.4,函数极值,7.1,线性方程组求解,7.1.1,直接解法,1,利用左除运算符的直接解法,对于线性方程组,Ax,=,b,，可以利用左除运算符“,”,求解：,x=A,b,例,7-1,用直接解法求解下列线性方程组。,命令如下：,A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;,b=13,-9,6,0;,x=Ab,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(Lb),或,x=U(L,Pb,),，这样可以大大提高运算速度。,例,7-2,用,LU,分解求解例,7-1,中的线性方程组。,命令如下：,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(Lb),或采用,LU,分解的第,2,种格式，命令如下：,L,U,P=,lu(A,);,x=U(LP*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(Qb),或,x=E(R(Qb),。,例,7-3,用,QR,分解求解例,7-1,中的线性方程组。,命令如下：,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(Qb),或采用,QR,分解的第,2,种格式，命令如下：,Q,R,E=,qr(A,);,x=E*(R(Qb),(3),Cholesky,分解,如果矩阵,X,是对称正定的，则,Cholesky,分解将矩阵,X,分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为,R,，则下三角矩阵为其转置，即,X=RR,。,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,的上三角阵，且满足,RR=X(1:q,1:q),。,实现,Cholesky,分解后，线性方程组,Ax=b,变成,RRx,=b,，所以,x=R(Rb),。,例,7-4,用,Cholesky,分解求解例,7-1,中的线性方程组。,命令如下：,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,),?Error using=,chol,Matrix must be positive definite,命令执行时，出现错误信息，说明,A,为非正定矩阵。,7.1.2,迭代解法,迭代解法非常适合求解大型系数矩阵的方程组。在数值分析中，迭代解法主要包括,Jacobi,迭代法、,Gauss-,Serdel,迭代法、超松弛迭代法和两步迭代法。,1,Jacobi,迭代法,对于线性方程组,Ax=b,，如果,A,为非奇异方阵，即,aii0(i=1,2,n),，则可将,A,分解为,A=D-L-U,，其中,D,为对角阵，其元素为,A,的对角元素，,L,与,U,为,A,的下三角阵和上三角阵，于是,Ax=b,化为：,x=D,-1,(L+U)x+D,-1,b,与之对应的迭代公式为：,x(k+1)=D,-1,(L+U)x(k)+D,-1,b,这就是,Jacobi,迭代公式。如果序列,x(k+1),收敛于,x,，则,x,必是方程,Ax=b,的解。,Jacobi,迭代法的,MATLAB,函数文件,Jacobi.m,如下：,function,y,n,=jacobi(A,b,x0,eps),if,nargin,=3,eps,=1.0e-6;,elseif,nargin,=,eps,x0=y;,y=B*x0+f;,n=n+1;,end,例,7-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),-1,Ux,(k),+(D-L),-1,b,该式即为,Gauss-,Serdel,迭代公式。和,Jacobi,迭代相比，,Gauss-,Serdel,迭代用新分量代替旧分量，精度会高些。,Gauss-,Serdel,迭代法的,MATLAB,函数文件,gauseidel.m,如下：,function,y,n,=gauseidel(A,b,x0,eps),if,nargin,=3,eps,=1.0e-6;,elseif,nargin,=,eps,x0=y;,y=G*x0+f;,n=n+1;,end,例,7-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),例,7-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),7.2,非线性方程数值求解,7.2.1,单变量非线性方程求解,在,MATLAB,中提供了一个,fzero,函数，可以用来求单变量非线性方程的根。该函数的调用格式为：,z=fzero(fname,x0,tol,trace),其中,fname,是待求根的函数文件名，,x0,为搜索的起点。一个函数可能有多个根，但,fzero,函数只给出离,x0,最近的那个根。,tol,控制结果的相对精度，缺省时取,tol,=,eps,，,trace,指定迭代信息是否在运算中显示，为,1,时显示，为,0,时不显示，缺省时取,trace=0,。,例,7-8,求,f(x,)=x-10,x,+2=0,在,x0=0.5,附近的根。,步骤如下：,(1),建立函数文件,funx.m,。,function,fx,=,funx(x,),fx,=x-10.x+2;,(2),调用,fzero,函数求根。,z=fzero(funx,0.5),z=,0.3758,7.2.2,非线性方程组的求解,对于非线性方程组,F(X)=0,，用,fsolve,函数求其数值解。,fsolve,函数的调用格式为：,X=fsolve(fun,X0,option),其中,X,为返回的解，,fun,是用于定义需求解的非线性方程组的函数文件名，,X0,是求根过程的初值，,option,为最优化工具箱的选项设定。最优化工具箱提供了,20,多个选项，用户可以使用,optimset,命令将它们显示出来。如果想改变其中某个选项，则可以调用,optimset,(),函数来完成。例如，,Display,选项决定函数调用时中间结果的显示方式，其中,off,为不显示，,iter,表示每步都显示，,final,只显示最终结果。,optimset(Display,off,),将设定,Display,选项为,off,。,例,7-9,求下列非线性方程组在,(0.5,0.5),附近的数值解。,(1),建立函数文件,myfun.m,。,function q=,myfun(p,),x=p(1);,y=p(2);,q(1)=x-0.6*sin(x)-0.3*,cos(y,);,q(2)=y-0.6*cos(x)+0.3*,sin(y,);,(2),在给定的初值,x0=0.5,y0=0.5,下，调用,fsolve,函数求方程的根。,x=fsolve(myfun,0.5,0.5,optimset(Display,off),x=,0.6354,0.3734,将求得的解代回原方程，可以检验结果是否正确，命令如下：,q=,myfun(x,),q=,1.0e-009*,0.2375 0.2957,可见得到了较高精度的结果。,7.3,常微分方程初值问题的数值解法,7.3.1,龙格库塔法简介,7.3.2,龙格库塔法的实现,基于龙格库塔法，,MATLAB,提供了求常微分方程数值解的函数，一般调用格式为：,t,y,=ode23(fname,tspan,y0),t,y,=ode45(fname,tspan,y0),其中,fname,是定义,f(t,y,),的函数文件名，该函数文件必须返回一个列向量。,tspan,形式为,t0,tf,表示求解区间。,y0,是初始状态列向量。,t,和,y,分别给出时间向量和相应的状态向量。,例,7-10,设有初值问题，试求其数值解，并与精确解相比较,(,精确解为,y(t,)=),。,(1),建立函数文件,funt.m,。,function,yp,=,funt(t,y,),yp,=(y2-t-2)/4/(t+1);,(2),求解微分方程。,t0=0;tf=10;,y0=2;,t,y,=ode23(funt,t0,tf,y0);%,求数值解,y1=sqrt(t+1)+1;%,求精确解,t,y,y1,y,为数值解，,y1,为精确值，显然两者近似。,例,7-11,求解著名的,Van,der,Pol,方程。,例,7-12,有,Lorenz,模型的状态方程，试绘制系统相平面图。,7.4,函数极值,MATLAB,提供了基于单纯形算法求解函数极值的函数,fmin,和,fmins,，它们分别用于单变量函数和多变量函数的最小值，其调用格式为：,x=fmin(fname,x1,x2),x=fm