sir matlab做相轨线,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 是初始状态列向量。

考虑到初始条件有

00

d , (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