实验四
常微分方程的数值解法 指令：
[t,y]=ode23(‘fun’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun’,tspan,yo) 高阶微分方程数值方法
其中fun是定义函数的文件名。该函数fun必须以为dx输出量，以t,y为输入量。tspan=[t0 tfina]表示积分的起始值和终止值。yo是初始状态列向量。
考虑到初始条件有
?dS???SI, S(0)?S0?0,??dt (5.24) ??dI??SI??I, I(0)?I?0.0??dt这就是Kermack与McKendrick的SIR仓室模型. 方程(5.24)无法求出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')
任务：
1 画出i（t），
2分析各参数的影响
例57：求解两点边值问题：xy???3y??x,y(1)?0,y(5)?0 。（注意：相应的数值解法比较复杂）。
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y =
-1/3*x^3+125/468+31/468*x^4
2
t2例：用数值积分的方法求解下列微分方程 y''?y?1?设初始时间t0=0；终止时间
2?tf=3*pi；初始条件y|x?0?0,y'|x?0?0。
解：(1)先将高阶微分方程转化为一阶微分方程组，其左端分别为变量x的两个元素的一阶导数。
t2设x1?y,x2?x1'?y',则方程可化为 x1'?x2 , x2'??x1?1?
2??x1'??01??x1??0???0??t2??01?t2?写成矩阵形式为 x'???????x???1???1?2??????10?x??1???1?2??? x'?10??2?????????2?????变量x的初始条件为x(0)=[0；0]，这就是待积分的微分方程组的标准形式。 (2)MATLAB程序
将导数表达式的右端写成一个exf.m函数程序，内容如下
function xdot=exf(t,x) u=1-(t.^2)/(pi^2);
xdot=[0 1;-1 0]*x+[0 1]'*u;
主程序调用中已有的数值积分函数进行积分，其内容如下
clf,t0=0;tf=3*pi;x0t=[0;0]; %给出初始值 [t,x]=ode23('exf',[t0,tf],x0t) %此处显示结果 y=x(:,1), %y为x的第二列
%本题的解析结果为y2(I)=(1+2/(pi^2))*(1-cos(t(I)))-t(I)^2/(pi^2) %在数值积分输出的时间序列点上计算它的值并画图与数值解作比较 for I=1;length(t);
y2(I)=(1+2/(pi^2))*(1-cos(t(I)))-t(I)^2/(pi^2); %解析解计算 end
u=1-(t.^2)/(pi^2);
clf,plot(t,y,'-',t,u,'+',t,y2,'o')
legend('数值积分解','输入量','解析解') %图例标注语
图6.13 上例中数值积分解与解析解的曲线
(4)结果分析
这个数值积分函数是按精度要求自动选择步长的。它的默认精度为10，因此图中的积分结果和解析解看不出差别。可以用长格式显示y和y2，比较他们的微小误差。若要改变精度要求，可在调用命令中增加可选变元。详情可通过help ode23查找。
练习
1. 传染病模型的各个解与图形
I?31 O 1/? 1 S 图 SIR模型的相轨线
2. . 食饵甲和捕食者乙在时刻t的数量分别记作x(t)，y(t)，当甲独立生存时它的（相对）
增长率为r，即x'?rx，而乙的存在使甲的增长率减小，设减小的程度与乙的数量成正比，于是x(t)满足方程
x'(t)?x(r?ay)?rx?axy， （１） 比例系数a反映捕食者掠取食饵的能力。
设乙独自存在时死亡率为d，即y'??dy，甲为乙提供食物相当于使乙的死亡率降低，并促使其增长。设这个作用与甲的数量成正比，于是y(t)满足方程 y'(t)?y(?d?bx)??dy?bxy， （２）
比例系数b反映食饵对捕食者的供养能力。 设食饵和捕食者的初始数量分别为
x(0)?x0,y(0)?y0。 （３）
设r?1,d?0.5,a?0.1,b?0.02,x0?25,y0?2,
求方程（１）、（２）在条件（３）下的数值解，并画出x(t),y(t)的图形。
002040608010012015103025205d2xdx(0)2dx?x?0在初始条件x(0)?1,3：求微分方程2??(1?x)?0情况下的解，并
dtdtdt图示。
?x1x2???1?4 x1(t)?r1x1?1?NN2?1?
21.81.61.41.210.80.60.40.20?xx?x2(t)?r2x2?1??21?2?N1N2??甲乙种群相轨线乙种群密度x200.20.40.60.811.2甲种群密度x11.41.61.82
5.. 已知某双种群生态系统的数学模型
x11x2?x(t)?2x(1???)1??110215 ?xx?x(t)?4x(1?3?1?2)22?1015?其中以x1(t),x2(t)分别表示t时刻甲乙两种群的数量，请问该模型表示哪类生态（相互竞争、相互依存）系统模型，并画出系统的相轨线图。
6. 拓展练习1 bvp
f'''?Mf'?ff''?f'2?0
f?0???0.，5f'?0??1?0.2*f''?0??0.3*f'''?0?，f(?
7. 拓展练习2
n?1)0
nf??f????2mn?m?1ff???m(d2?f?2)?M(d?f?)?0
n?1nf(0)?0, f?(0)?1??f??(0), f?(?)?0.5
m.n取常数