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

常微分方程的数值解法 指令：

[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)I0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.5，f'010.2*f''00.3*f'''0，f

7. 拓展练习2

n10

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取常数