MATLAB解方程和函数极值.ppt

第三讲解方程与函数极值,方程问题和极值问题始终是数学问题中的核心问题！,学习内容,多项式运算线性方程组求解非线性方程数值求解常微分方程初值问题的数值解法函数极值线性插值,（一）多项式的表示方法,对于多项式的表达式约定如下对于多项式对于上述多项式一般用以下行向量表示,（二）多项式的创建,1.系数向量直接输入法由于MATLAB中多项式是以向量形式存储的,因此最简单的多项式输入即向量输入。例:输入多项式p=1-56-33；poly2sym(p)%此函数将多项式转换换为符号多项式,（二）多项式的创建,2.特征多项式输入法多项式创建的另一途径是从矩阵求其特征多项式获得，由函数poly实现例如：a=123;234;345p1=poly(a)%矩阵a对应的特征多项式poly2sym(p1)%将多项式p1转换为符号多项式,（二）多项式的创建,3.由根创立多项式给定的根也可产生其相应的多项式,此功能还是由函数poly实现例：root=-5-3+4i-3-4i;p=poly(root)poly2sym(p)注：若要生成实系数多项式,则根中的复数比为对称共轭复数,（三）多项式运算,1.求多项式的值求多项式的值可以由两种形式(1)输入变量值代入多项式计算时是以数组为单元的,此时函数为polyval(2)以矩阵为计算单元,进行矩阵式运算以求的多项式的值,此时的函数为polyvalm注：这两种计算在数值上由很大差别。当进行矩阵运算时,变量矩阵须为方阵。,实例演示,p=11155125;b=11;11;c=5poly2sym(p)polyval(p,b)poly(p,c)polyvalm(p,b),（三）多项式运算,2.求多项式的根求多项式的根可以由两种方法(1)直接调用MATLAB的函数roots求解多项式的所有根(2)通过建立多项式的伴随矩阵再求其特征值的方法得到多项式的所有根,实例演示,用两种方法解方程的所有根p=2-56-19;roots(p)a=compan(p)%求多项式的特征矩阵eig(a)%求特征矩阵a的特征根,（三）多项式运算,3.多项式的乘除法(1)多项式的乘法由函数conv来实现,此函数同用于向量的卷积(2)多项式的除法由函数deconv来实现,此函数与向量的解卷函数相同,实例演示,计算两多项式的乘除法p=2-56-19;poly2sym(p);d=3-90-18;ploy2sym(d)pd=conv(p,d)%多项式p与d相乘poly2sym(pd)pl=deconv(pd,d)%多项式pd除以d,（三）多项式运算,4.多项式微分多项式的微分可以用函数polyer进行例：p=2-56-19；poly2sym(p);Dp=polyer(p),（三）多项式运算,5.多项式的拟合多项式拟合是多项式运算的一个重要组成部分，在工程及科研工作中得到了广泛的应用，其一方面可以由矩阵的除法解超定方程来进行;另一方面在MATLAB中还提供了专门的拟合函数polyfit,其调用格式如下：(1)polyfit(X,Y,n)%其中X,Y为拟合数据,n为你和多项式的阶数(2)p,s=polyfit(X,Y,n)%其中p为拟合多项式的系数向量，s为拟合多项式系数向量的结构信息,实例演示,例:用五阶多项式对0,pi/2上的正弦函数进行最小二乘拟合x=0:pi/20:pi/2;y=sin(x);a=ployfit(x,y,5);%用五阶多项式拟合y=sin(x),a为你和多项式的系数x1=0:pi/30:pi/2;y1=sin(x1);y2=a(1)*x1.5+a(2)*x1.4+a(3)*x1.3+a(4)*x1.2+a(5)*x1.+a(6)plot(x1,y1,b-,x2,y2,r*)legend(原曲线，拟合曲线)axis(0,7,-1.2,4)%坐标轴的限制,学习内容,多项式运算线性方程组求解非线性方程数值求解常微分方程初值问题的数值解法函数极值线性插值,线性方程组解的结构,一、齐次线性方程组解的结构齐次线性方程组的矩阵形式为AX=0其中A是mn阶矩阵；X为未知向量,二、线性方程组求解,1.直接解法（1）利用左除运算符的直接解法对于线性方程组Ax=b，可以利用左除运算符“”求解:x=Ab例:用直接解法求解下列线性方程组。命令如下：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分解等下面将对他们一一介绍,a)特征值分解,矩阵的特征值分解也调用函数eig,其调用格式为(1)V,D=eig(X)%得到矩阵X的特征值对角矩阵D和其列为对应特征值的特征向量矩阵V,矩阵的特征值分解为XV=VD(2)V,D=eig(X,nobalance)%此形式为关闭平衡算法的求解方法.平衡算法对于某些问题可以得到更高的准确度。(3)V,D=eig(A,B)%对矩阵A和B做广义特征值分解,即AV=BVD,实例演示-矩阵的特征值分解,单位阵的特征值分解a=-149-50-154;537180546;-27-9-25v,d=eig(a)双矩阵的特征值分解b=2102;105-8;2-811;v,d=eig(a,b),附:复数特征值对角阵与实数特征值对角阵的转换,即使对于实阵,其特征值也可能出现复数.而在实际使用中,常需要把这些共轭复数特征值转换为实数块.MATLAB提供两个函数进行转换(1)V,D=cdf2rdf(V,D)%将复数对角型转换为实数块对角型(2)U,T=rdf2cdf(U,T)%将实数对角块转换为复数对角型,实例演示-复数对角型矩阵与实数对角型矩阵的转换,a=1-3;22/3;v,d=eig(a)%将矩阵a进行特征值分解,其结果复对角型vs,ds=cdf2rdf(v,d)%将复对角型转换为实对角型b=vs*ds/vs%其结果与a相等,b)奇异值分解,设矩阵A为一个mn阶的(实)矩阵(mn),则存在正交矩阵V和U,使得:此式称为A的奇异值分解矩阵奇异值分解可由函数svd实现,其调用格式为(1)U,S,V=svd(X)(2)U,S,V=svd(X,0)作用:生成U,S,和V使得X=USV实例演示例：对矩阵a进行奇异值分解a=1;1;U,S,V=svd(a)注：进行奇异值分解时,矩阵的维数没有限制,c)LU分解,LU分解是高斯消去法的基础,它将一矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式.线性代数中已经证明,只要方阵A是非奇异的,LU分解总是可以进行的。lu函数用于对矩阵进行LU分解，其调用格式为(1)L,U=lu(X)-产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意,这里的矩阵X必须是方阵。(2)L,U,P=lu(X)-产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组Ax=b的解x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。,实例演示,例:用LU分解求解例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),d)QR分解,在求解矩阵的特征值时,引入一种分解方法,即实阵A可以写成A=QR,其中Q为一个正交矩阵和R为一个上三角矩阵.QR分解只能对方阵进行.QR分解由函数qr实现，其调用格式为:(1)Q,R=qr(X):产生一个一个正交矩阵Q和一个上三角矩阵R,使之满足X=QR。(2)Q,R,E=qr(X):产生一个一个正交矩阵Q,一个上三角矩阵R以及一个置换矩阵E,使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解x=R(Qb)或x=E(R(Qb)。,实例演示,例:用QR分解求解前例中的线性方程组。命令如下：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),e)Cholesky分解,如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一个对角元素为正的下三角矩阵R和上三角矩阵R的乘积.使得X=RR.函数chol(X)用于对矩阵X进行Cholesky分解,其调用格式为:(1)R=chol(X):产生一个上三角阵R,使RR=X.若X为非对称正定,则输出一个出错信息。(2)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)。,实例演示,例:用Cholesky分解求解前例中的线性方程组。命令如下：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为非正定矩阵?可以执行命令R,p=chol(X)x=R(Rb),2.迭代解法,迭代解法非常适合求解大型系数矩阵的方程组.在数值分析中,迭代解法主要包括:Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。1Jacobi迭代法:对于线性方程组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-1b与之对应的迭代公式为:X(k+1)=D-1(L+U)x(k)+D-1b这就是Jacobi迭代公式.如果序列x(k+1)收敛于x，则x必是方程Ax=b的解,附：Jacobi迭代法的MATLAB函数文件Jacobi.m如下：,functiony,n=jacobi(A,b,x0,eps)ifnargin=3eps=1.0e-6;elseifnargin=epsx0=y;y=B*x0+f;n=n+1;end,例:用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)！其他迭代法在此一律不在详细说明！,学习内容,线性方程组求解非线性方程数值求解常微分方程初值问题的数值解法函数极值线性插值,二、非线性方程组的求解,对于非线性方程组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.,例:求下列非线性方程组在(0.5,0.5)附近的数值解。(1)建立函数文件myfun.m。functionq=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.63540.3734,学习内容,线性方程组求解非线性方程数值求解常微分方程初值问题的数值解法函数极值线性插值,三、常微分方程初值问题的解法,一、解析解求解微风方程（组）的解析解的MATLAB命令为dsolve(eqn1,eqn2,x)其中eqni表示第i个方程与初始条件等式x表示微风方程（组）中的自变量，默认是自变量为t,实例演示,例1：求解一阶微分方程的通解及x=0,y=1时的特解求通解：dsolve(Dy=1+y2,x)求特解：dsolve(Dy=1+y2,y(0)=1,x)例2：求解二阶微分方程,实例演示,输入命令：dsolve(x2*D2y+x*Dy+(x2-(1/2)2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x)结果为：pi(1/2)*2(1/2)/x(1/2)*sin(x)上机练习：求方程组,二、数值解,微分方程数值解的求法有多种,常用的主要是RungeKutta(龙格-库塔法)和欧拉方法基于龙格库塔法，MATLAB提供了使用2阶（3阶）或4阶（5阶）龙格库塔公式，一般调用格式为：t,y=ode23(fname,tspan,y0，options)t,y=ode45(fname,tspan,y0,options)其中fname是定义f(t,y)的函数文件名，tspan的取法有几种，当tspan=t0,tf时,t0,tf分别表示自变量的初值和终值；当tspan=t0,t1,tf时则输出在指定时刻t0,t1,tf处给出；当tspan=t0:k:tf时则输出在t0,tf的等分点给出y0是函数的初值。t,y为输出矩阵，分别表示自变量t和因变量y的取值,实例演示,例求方程的数值解首先建立M函数文件：jie3.mfunctionf=jie3(x,y)f=x2*D2y+x*Dy+(x2-(1/2)2)*y=0x,y=ode23(jie3,pi/2,pi),2,-2/pi,实例演示,用经典的R-K方法求解解：编制函数文件：fun