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 是初始状态列向量。

考虑到初始条件有

00d , (0)0,d d , (0)0.d S

SI S S t

I SI I I I t

ββμ⎧=-=>⎪⎪⎨

⎪=-=≥⎪⎩ (5.24) 这就是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：求解两点边值问题：0)5(,0)1(,32

==='-''y y x y y x 。（注意：相应的数值解法比较复杂）。

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

例：用数值积分的方法求解下列微分方程 π

21''2

t y y -=+设初始时间t0=0；终止时间

tf=3*pi ；初始条件0|',

0|00====x x y y 。

解：(1)先将高阶微分方程转化为一阶微分方程组，其左端分别为变量x 的两个元素的一阶导数。

设,'',121y x x y x ===则方程可化为 21'x x = , π

21'2

12t x x -+-=

写成矩阵形式为 ⎪⎪⎭

⎫

⎝⎛-⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡-=⎪⎪⎭⎫ ⎝⎛-⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡=ππ2110011021100110'''222121t x t x x x x x 变量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)结果分析

这个数值积分函数是按精度要求自动选择步长的。它的默认精度为3

10-，因此图中的积分结果和解析解看不出差别。可以用长格式显示y 和y2，比较他们的微小误差。若要改变精度要求，可在调用命令中增加可选变元。详情可通过help ode23查找。

练习

1. 传染病模型的各个解与图形

2. . 食饵甲和捕食者乙在时刻t 的数量分别记作)(t x ，)(t y ，当甲独立生存时它的（相对）

增长率为r ，即rx x ='，而乙的存在使甲的增长率减小，设减小的程度与乙的数量成正比，于是)(t x 满足方程

axy rx ay r x t x -=-=)()('， （１） 比例系数a 反映捕食者掠取食饵的能力。

设乙独自存在时死亡率为d ，即dy y -='，甲为乙提供食物相当于使乙的死亡率降低，并促使其增长。设这个作用与甲的数量成正比，于是)(t y 满足方程 bxy dy bx d y t y +-=+-=)()('， （２）

比例系数b 反映食饵对捕食者的供养能力。 设食饵和捕食者的初始数量分别为

图 SIR 模型的相轨线

I

O

1

0)0(,)0(y y x x ==。 （３）

设,2,25,02.0,1.0,5.0,100======y x b a d r

求方程（１）、（２）在条件（３）下的数值解，并画出)(),(t y t x 的图形。

3：求微分方程0)1(222=+--x dt dx x dt x d μ在初始条件0)0(,1)0(==dt

dx x 情况下的解，并

图示。

4

00.20.40.6

0.81 1.2

1.4 1.6 1.82

0.20.40.60.811.21.41.61.82甲种群密度x1

乙种群密度x 2

甲乙种群相轨线

5.. 已知某双种群生态系统的数学模型

121112221()2(1)10215

()4(13)

1015x x x t x x x x t x ⎧

=--⨯⎪⎪⎨

⎪=-⨯-⎪⎩

其中以12(),()x t x t 分别表示t 时刻甲乙两种群的数量，请问该模型表示哪类生态（相互竞争、相互依存）系统模型，并画出系统的相轨线图。

12111112()1x x x t r x N N σ⎛⎫

=-- ⎪⎝

⎭12222212()1x x x t r x N N σ⎛⎫

=-- ⎪

⎝

⎭