MATLAB 应用 第十章 常微分方程.ppt

MATLAB,高等数学实验,实验十常微分方程,实验目的掌握用MATLAB求微分方程及方程组解析解的方法,学习求微分方程近似解的欧拉折线法和龙格-库塔法。,10.1学习MATLAB命令,10.1.1求常微分方程的符号解函数,dsolve(eq1,eq2,cond1,cond2,v)输入符号形式的常微分方程(组)eq1,eq2,及初始条件cond1,cond2,;默认的自变量是t,如果要指定自变量v,则在方程组及初始条件后面加v,并用逗号分开。在常微分方程(组)的表达式eqn中,大写字母D表示对自变量(默认为t)的微分算子：D=d/dt,D2=d2/dt2,。算子D后面的字母则表示因变量,即待求解的未知函数。返回的结果中可能会出现任意常数C1,C2等。若命令找不到解析解,则返回一警告信息。,例如,输入下面的方程都可以得到其解：dsolve(Dx=-a*x)%ans=C2/exp(a*t)x=dsolve(Dx=-a*x,x(0)=1,s)%输出x=1/exp(a*s)w=dsolve(D3w=-w,w(0)=1,Dw(0)=0,D2w(0)=0)%输出略f,g=dsolve(Df=f+g,Dg=-f+g,f(0)=1,g(0)=2)%输出略,10.1.2求常微分方程组初值问题的数值解函数,T,Y=ode23(odefun,tspan,y0)T,Y=ode45(odefun,tspan,y0)这两个命令都是求解导数已经解出的一阶微分方程y=f(t,y)满足初始条件y0的数值解。输入参数odefun为M文件定义的函数或inline格式的函数f(t,y);tspan=t0,tf表示求解区间(从t0到tf)。输出的解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时问点t0,t1,t2,上的解,则令tspan=t0,t1,t2,tf(要求是单调的)。,10.2实验内容,10.2.1求微分方程的解析解,对于可以用积分方法求解的微分方程和微分方程组,可用dsolve命令来求其通解或特解。例如要求方程y+y-2y=0的通解,通过输入：y=dsolve(D2y+D1y-2*y=0,x)%一阶导数符号是D1y或Dy,二阶导数符号是D2y执行以后就能得到含有两个任意常数C12和C13的通解：y=C12*exp(x)+C13/exp(2*x)如果想求解微分方程的初值问题：,则只要输入：y=dsolve(D2y+4*Dy+3*y=0,y(0)=6,Dy(0)=10,x)%用单引号对把方程和初始条件分别括起来输出为：y=14/exp(x)-8/exp(3*x),【例1】求微分方程的通解。输入：clear;dsolve(Dy+2*x*y=x*exp(-x2),x)便得到微分方程的通解：ans=C1/exp(x2)+x2/(2*exp(x2)其中C1是任意常数。,【例2】求微分方程在初始条件下的特解。输入：clear;dsolve(x*Dy+y-exp(-x)=0,y(1)=2*exp(1),x)就可以得到特解：ans=(1/exp(1)-1/exp(x)+2*exp(1)/x,【例3】求微分方程的通解。输入：clear;simplify(dsolve(D2y-2*Dy+5*y=exp(x)*cos(2*x),x)便得到通解：(exp(x)*(x*sin(2*x)+4*C13*cos(2*x)+4*C14*sin(2*x)/4,解微分方程组时的命令格式为x,y=dsolve(Dx=f(x,y),Dy=g(x,y)或s=dsolve(Dx=f(x,y),Dy=g(x,y),【例4】求方程组在条件下的特解。输入：clear;x,y=dsolve(Dx=-x-2*y+exp(t),Dy=x+y,x(0)=1,y(0)=0)或x,y=dsolve(Dx=-x-2*y+exp(t),Dy=x+y,x(0)=1,y(0)=0)则得到特解：x=cos(t)y=exp(t)/2-cos(t)/2+sin(t)/2,10.2.2欧拉折线法,首先扼要介绍欧拉折线法的思想。给定微分方程和初始条件,考虑的线性近似：如果在包含的很短的区间内,则函数是的很好近似。欧拉折线法是通过一系列线性近似得到在较长区间内的近似解。设,如果增量很小,则：是对精确值的很好近似。所以从曲线上的点出发首先得到与非常近似的。再利用和斜率,可以进行下一步。,设,再通过线性近似：得到的近似值。第三步,由和斜率得到：连接,的折线就是微分方程满足初始条件的一个近似解。因此用欧拉折线法求近似解的一般步骤是：MATLAB可以轻易完成这些计算。,【例5】用欧拉折线法解初值问题：,(1)计算过程是：,所以用MATLAB进行计算时,输入：clearclfszy=;y=1;szy=szy;y;forx=0.1:0.1:1y=y+(1+y)*0.1;szy=szy;y;endszy输出为用欧拉折线法算出的初值问题(1)在结点x=0:0.1:1处的数值解(存放在数组szy中)。,szy=1.00001.20001.42001.66201.92822.22102.54312.89743.28723.71594.1875,而输入：y=dsolve(Dy=1+y,y(0)=1,x)可以算得初值问题(1)的精确解是：y=2*exp(x)-1为了比较精确解和近似解的误差,编写程序算出精确解在结点x=0:0.1:1处相应的纵坐标(存放在数组jqy中)。jqy=;forx=0:0.1:1y=-1+2.*exp(x);jqy=jqy;y;endjqy,输出为：jqy=1.00001.21031.44281.69971.98362.29742.64423.02753.45113.91924.4366,接下来,把解曲线和折线在同一坐标系中画出来,输入：plot(szy,b)holdonplot(jqy,r)得到图10.1。,图10.1,从输出容易观察到：欧拉折线法在经过多步以后,其误差会积累起来。输入：err=jqy-szy输出为在0,0.1,0.2,1.0处精确解和近似解的差。err=00.01030.02280.03770.05540.07640.10110.13010.16390.20330.2491,可以看到在x=1处,相对误差大约为5.6%。为了减少误差,一个方法是减少步长dx的值。另外在近似算法上也有很多改进办法。例如用改进的欧拉折线法可以算出初值问题(1)在结点x=0:0.1:1处的数值解(存放在数组gjszy中)。x=0:0.1:1;h=0.1;gjszy=;y=1;gjszy=gjszy;y;fori=1:1:length(x)-1y1=y+(1+y)*h;y=y+(1+y+1+y1)*h/2;gjszy=gjszy;y;endgjszyplot(x,gjszy,k+,x,jqy,r),输出为：gjszy=1.00001.21001.44211.69851.98182.29492.64093.02313.44563.91244.4282结果见图10.2。,图10.2,用四阶龙格-库塔法,输入：x=0:0.1:1;h=0.1;szyode=;y=1;szyode=szyode;y;fori=1:1:length(x)-1k1=1+y;k2=1+y+k1*h/2;k3=1+y+k2*h/2;k4=1+y+k3*h;y=y+(k1+2*k2+2*k3+k4)*h/6;szyode=szyode;y;endszyodeplot(x,szyode,k+,x,jqy,r),输出为：szyode=1.00001.21031.44281.69971.98362.29742.64423.02753.45113.91924.4366结果见图10.3。,图10.3,10.2.3求微分方程的数值解,对于不可以用积分方法求解的微分方程初值问题,可以用二三阶龙格-库塔法ode23或四五阶龙格-库塔法ode45命令来求其数值解。例如要求方程,的近似解(),输入：fun=inline(y2+x3,x,y);ode23(fun,0,1.5,0.5)%绘制初值问题的数值解曲线,命令中的0,1.5表示x相应的区间,0.5表示y的初值输出见图10.4,图中圆圈表示计算过程中选取的结点。,图10.4,要得到结点处的坐标,输入：x,y=ode23(fun,0,1.5,0.5);x,y%显示结点处的坐标输出为：ans=00.50000.15000.54070.30000.59040.45000.65660.60000.75260.74080.88930.86581.07350.97631.31641.07321.62851.15702.02081.22872.50601.28953.09941.34063.82061.38334.69461.41885.75221.44847.03211.47288.58131.493010.45781.500011.3080,因为用ode23或ode45命令得到的输出是解的近似值。首先在区间0,1.5内插入一系列点,计算出在这些点上函数的近似值,再通过插值方法得到在区间上的近似解。因此ode23或ode45命令的输出只能用省略的形式给出。可以通过作图对函数作进一步的了解。,【例6】求初值问题在区间1.2,4上的近似解,并作图。输入：fun=inline(1+x*y)*y/(x*y-1),x,y);x,y=ode45(fun,1.2,4,1);x,yode45(fun,1.2,4,1)%或plot(x,y)其输出为数值近似解(插值函数)的形式：,ans=1.20001.00001.20461.04561.20911.08451.21371.11941.21831.15181.22741.21141.23661.26581.24571.31681.25491.36543.829035.11513.886037.20593.943039.41894.000041.7613并得到近似解的图形(见图10.5)。,图10.5,【例7】求VanderPol方程,在区间0,20上的近似解。令,原方程化为方程组。单击File：New：M-file打开命令窗口,输入以下程序,建立并保存函数文件vdp1.m：functiondydt=vdp1(t,y)dydt=y(2);(1-y(1)2)*y(2)-y(1);回到MTLAB命令窗口,输入求VanderPol方程近似解命令：t,y=ode23(vdp1,020,0-0.5);plot(t,y(:,1);可以观察到近似解的图形(见图10.6)。,图10.6,10.2.4一阶微分方程的斜率场,设一阶微分方程的形式为,即研究导数已经解出的一阶方程。其中是已知函数。由导数的几何意义,未知函数y的斜率是平面上点(x,y)的函数。因此可以在xoy平面上的每一点,作出过该点的以为斜率的一条很短的直线(即是未知函数y的切线)。这样得到的一个图形就是微分方程的斜率场。为了便于观察,实际上只要在xoy平面上取适当多的点,作出在这些点的函数的切线。观察这样的斜率场,就可以得到方程的近似积分曲线。,【例8】绘制微分方程的斜率场和解曲线。因为没有给出初始条件,所以先建立并保存成函数文件zxyf.m。functionDy=zxyf(x,y)Dy=x./y;%计算斜率函数编写命令文件：clf,clear;a=0;b=4;c=0;d=4;n=15;x,y=meshgrid(linspace(a,b,n),linspace(c,d,n);%生成区域中的网格f=x./y;Fx=cos(atan(x./y);Fy=sqrt(1-Fx.2);%计算切线斜率矢量quiver(x,y,Fx,Fy,0.5)%在每一网格点画出相应的切线斜率矢量,0.5是控制矢量大小的参数holdonaxis(a,b,c,d)输出见图10.7。,图10.7,如果增加初始条件y(0)=2,则得到的是一条确定的