matlab常微分方程的数值解法试验报告

1、实验四常微分方程的数值解法 指令：t,y=ode23( fun ,tspan,yo)2/3阶龙格库塔方法t,y=ode45( fun ,tspan,yo)t,y=ode113( fUn ,tspan,yo)4/5阶龙格库塔方法高阶微分方程数值方法其中fun是定义函数的文件名.该函数fun必须以为dx输出量,以t,y为输入量.tspan=t0tfina表示积分的起始值和终止值.yo是初始状态列向量.考虑到初始条件有=- ：SI, S(0) = S00,(5.24)dt史=SI -I , I(0) = I0 一0.dt0这就是Kermack与McKendrick的SIR仓室模型.方程(5.24)无

2、法求出S(t)和I (t)的解析解.我们先做数值计算.Matlab代码为：function dy=rigid(t,y)dy=zeros(2,1);a=1;b=0.3;dy(1)=a*y(1).*y(2)-b*y(1);dy(2)=-a*y(1).*y(2);ts=0:.5:50;x0=0.02,0.98;T,Y=ode45( 'rigid' ,ts,x0);%plot(T,Y(:,1),'-',T,Y(:,2),'*')plot(Y(:,2),Y(:,1),'b-')xlabel('s')ylabel('i

3、')任务：1画出i (t),2分析各参数的影响例57：求解两点边值问题：xy"-3y' = x2, y(1) = Q y(5) = 0 .(注意：相应的数值解法比拟复杂).y=dsolve('x*D2y-3*Dy=xA2','y(1)=0,y(5)=0','x')/y =-1/3*xA3+125/468+31/468*xA4t2、例：用数值积分的方法求解以下微分方程t0=0;终止时间y''+y=1 设初始时间2 二tf=3*pi ；初始条件 y |x乂 = 0,y'|x旦=0.解：(1)先将高阶微

4、分方程转化为一阶微分方程组,其左端分别为变量x的两个元素的一阶导数.、,t2设 x1二y, x2 = x1'= y',那么方程可化为x1' = x2 ,x2'=-x1+1 2 二写成矩阵形式为x'=l0 Tl+FkCM0 "+即一口 Ml1-1 0_%Ml2DL-1 0 - Ml 2J变量x的初始条件为x(0)=0 ； 0,这就是待积分的微分方程组的标准形式.(2)MATLAB 程序将导数表达式的右端写成一个exf.m函数程序,内容如下function xdot=exf(t,x)u=1-(t.A2)/(piA2);xdot=0 1;-1 0*x

5、+0 1'*u;主程序调用中已有的数值积分函数进行积分,其内容如下clf,t0=0;tf=3*pi;x0t=0;0;惭出初始值t,x=ode23( 'exf' ,t0,tf,x0t)%t 匕处显示结果y=x(:,1),%丫为* 的第二列泌题的解析结果为 y2(I)=(1+2/(piA2)*(1-cos(t(I)-t(I)A2/(piA2)%在数值积分输出的时间序列点上计算它的值并画图与数值解作比拟for I=1;length(t);y2(I)=(1+2/(piA2)*(1-cos(t(I)-t(I)A2/(piA2);%解析解计算endu=1-(t.A2)/(piA2)

6、;clf,plot(t,y,'-',t,u,'+' ,t,y2,'o')legend('数值积分解,输入量,解析解)图例标注语_gJLJLUD11 力012345673910图6.13上例中数值积分解与解析解的曲线(4)结果分析这个数值积分函数是按精度要求自动选择步长的.它的默认精度为10,因此图中的积分结果和解析解看不出差异.可以用长格式显示y和y2,比拟他们的微小误差.假设要改变精度要求,可在调用命令中增加可选变元.详情可通过help ode23查找.练习1 .传染病模型的各个解与图形I图SIR模型的相轨线1.1. 饵甲和捕食者乙在时

7、刻t的数量分别记作 x(t) , y(t),当甲独立生存时它的(相对)增长率为r ,即x' = rx ,而乙的存在使甲的增长率减小,设减小的程度与乙的数量成正 比,于是 x(t)满足方程x'(t) =x(r -ay) =rx axy ,( i)比例系数a反映捕食者掠取食饵的水平.设乙单独存在时死亡率为 d ,即y'=-dy,甲为乙提供食物相当于使乙的死亡率降低,并促使其增长.设这个作用与甲的数量成正比,于是y(t)满足方程y'(t) = y(-d +bx) = -dy +bxy ,(2)比例系数b反映食饵对捕食者的供养水平.设食饵和捕食者的初始数量分别为x(0)

8、 =X0,y(0) = o(3 )Xt, yt的图形.设 r =1耳 =0.5, aO.lb = 0.02,X.= 25, y0 =2, 求方程1 、 2 在条件3 下的数值解,并画出x(0)=1,改6 =0情况下的解,并 dtd 2x ,2 dX3:求微分方程 -(1 -X )+x=0在初始条件 dt2dt图示.XL 一二NiX21 N2X2 (t) =2X2 1 仃XiX22.NiN2J5.某双种群生态系统的数学模型X1=2X1(1 端-：如t3箕相互竞争、其中以Xit,X2t分别表示t时刻甲乙两种群的数量, 请问该模型表示哪类生态 相互依存系统模型,并画出系统的相轨线图.6 .拓展练习1 bvp-2 一f '''_Mf ' ff ''_ f '2 =0f (0)=-0. 5 f'(0)=1+0.2*f (0)+0.3*f '(0)