常微分方程和人口模型实验.docVIP

常微分方程和人口模型实验 实验8 常微分方程和人口模型 实验目的: (1)了解常微分方程的基本概念。 (2)了解常微分方程的解析解。 (3)了解常微分方程的数值解。 (4)学习、掌握MATLAB软件有关的命令。 实验内容: 常微分方程模型的建立及求解。 一、常微分方程的解析解 调用格式:dsolve(S,s1,s2,...,x) 其中，S为方程或方程组，方程S中用D表示求导数，D2，D3，…表示二阶、三阶等高阶导数，方程间用逗号分隔;s1，s2,…为初始条件，x为自变量，初始条件和自变量都可以省略;初始条件缺省时给出带任意常数C1，C2，…的通解;自变量缺省时设定为t(表8.1给出了这个命令的基本用法，请读者运行表中的命令，并观察所得结果。 表8.1 dsolve命令用法 常微分方程问题 求解命令 2yy1,， y=dsolve(?Dy=1+y^2?) 2,yy1,， y=dsolve(?Dy=1+y^2?,?y(0)=1?,?x?) ,y(0)1,, txxxe22，，, x=dsolve(?D2x+2*D1x+2*x=exp(t)?,?x(0)=1?,?Dx(0)=0?) 2,xyyx3,, y=dsolve(?x*D2y-3*Dy=x^2?,?y(1)=0,y(0)=0?,?x?) ,yy(1)(0)0,,, xxy34,，, [x,y]=dsolve(?Dx=3*x+4*y?,?Dy=-4*x+3*y?) ,yxy43,,，, xxy,，, , yxy,,[x,y]=dsolve(?Dx=x+y,Dy=y-x?,?Dx(0)=1,Dy(0)=1?) , ,xy(0)(0)1,,, 二、常微分方程的数值解 在自然科学的众多领域都会遇到常微分方程的求解问题，然而，我们知道，只是有少数十分简单的常微分方程可用简单的方法给出解析解，大部分微分方程只能用数值解法求近似解。下面用一个例子简要介绍常微分方程数值解法的基本思想。 考虑一阶微分方程初值问题 dy,(,()),,,,ftytttt,0f (8.1) dt, ,()yty,00, t,t,？,t,t所谓数值解法，就是寻求上述初值问题的解在一系列离散节点处的01nf yt()yyy,,,h,t,thttn,,()值的的近似值。称为步长，通常取为常量。kk，1kk01nf0最简单的数值解法是欧拉法。欧拉法的思路很简单:在节点处用差商近似代替导数 ,1ythytyt()[()()],, kkk，1 将上式代入(8.1)就得到近似计算公式(称为欧拉公式) yyhftyk,，,(,),0,1,2, kkkk，1 根据初值问题对步长的要求，把常微分方程进一步划分为刚性方程和非刚性方程两类。所谓刚性方程，就是只有步长充分小时才能得到满意近似解的方程。 MATLAB提供了若干直接求解一阶常微分方程(组)初值问题的命令，调用格式如下: [x,y]=odeN(?funname?,[t0,tf],y0,tol,trace) 其中，odeN可以是ode23,ode45,ode113,ode15s,ode23s中的任意一个命令;?funname?是定义 设置解的精函数的文件名;[t0,tf]是常微分方程的求解区间;y0是初始状态列向量;tol用于度，tol可以缺省，缺省时，在ode23中tol=1E-3，在ode45中为tol=1E-6;trace设定是否显示求解的中间结果，缺省值为trace=0，表示不显示中间结果。 注意:命令ode23,ode45,ode113用于求解非刚性方程，命令ode15s,ode23s用于求解刚性方程;这五个命令采用不同的算法，是用于求解不同的方程，表8.2简要给出了这些命令介绍。 表8.2 odeN命令的简要介绍 ode23 基于Runge-Kutta(2,3)法，单步法，低阶，精度较低，速度较快。 Ode45 基于Runge-Dutta(4,5)，单步法，中阶，精度较高，速度较慢。 Ode113 采用Adams-Bashforth-moullton法，各(可变)阶方法，多步法。 Ode15s 采用多步法，精度较低。 Ode23s 采用二阶改进的Rosenbrock法，单步法，速度较快。 对于高阶微分方程，必须先用换元法将其化为化为一阶微分方程组，即“状态方程”，再加以调用Matlab命令求解。 练习1 求下列微分方程的数值解，并作解函数的曲线图。 2x,y,y,,0,x,20,y , ,y(0),1, 1) 建立M函数文件dyex801.m function dy=dyex801(x,y) dy=[y-2*x/y]; 2) 调用程序文件ex801.m [x,y]=ode45(dyex801,[0