微分方程数值解(生物种群的相互竞争模型)

1．求微分方程的解析解, 并画出它们的图形， (1)y ’= y + 2x , y (0) = 1, 0<x <1;
用matlab 编程如下： syms x y k
k=dsolve('Dy=y+2*x','y(0)=1','x') ezplot(k,[0,1]) 其结果和图像如下：
ans=-2*x-2+3*exp(x)
2. y ’’+y cos(x ) = 0, y (0)=1, y ’(0)=0; 用matlab 编程如下： syms x y z
z=dsolve('D2y+y*cos*(x)=0','y(0)=1','Dy(0)=0') ezplot('z')
其图像和结果如下：
z =cos(cos^(1/2)*x^(1/2)*t)
3．Rossler 微分方程组：
⎪⎩
⎪
⎨⎧-+=+=--=)('''c x z b z ay x y z y x
当固定参数b=2, c=4时，试讨论随参数a 由小到大变化（如a ∈(0,0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？ 用matlab 编程
建立rossler.m 文件： function r=rossler(t,x) global a; global b; global c;
r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)]; 建立exp4-3.m 如下： global a; global b; global c; b=2; c=4;
t0=[0,200];
for a=0:0.03:0.65
[t,x]=ode45('rossler',t0,[0,0,0]); subplot(1,2,1);
plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');
title('x(红色),y(绿色),z(蓝色)随t 变化情况');xlabel('t'); subplot(1,2,2);
plot3(x(:,1),x(:,2),x(:,3))
title('相图');xlabel('x');ylabel('y');zlabel('z'); end
(1)当a=0时，图像如下：
所以当a=0时，(x,y,z)收敛于（0,0.5,0.5） (2)当a=0.12时，图像如下：
x(红色),y(绿色),z(篮色)随t 变化情况t
相图
z
当a=0.27时，(x,y,z)仍然收敛，但是收敛速度大大降低。

(3)当a=0.27时，图像如下
(4)当a=0.39时，图像如下
(5)当a=0.52时，图像如下
t
z
x(红色),y(绿色),z(篮色)随t 变化情况t
相图
z
x(红色),y(绿色),z(篮色)随t 变化情况t
相图
z
从这一系列变化的图像中可以看出，随着 a 的增大，(x,y,z)接近极限环的速度加快。

4.Apollo 卫星的运动轨迹的绘制
用matlab 编程如下： x0=[1.2;0;0;-1.04935751]; options=odeset('reltol',1e-8);
[t,y]=ode45('appollo',[0,20],x0,options); plot(y(:,1),y(:,3))
title('Appollo 卫星运动轨迹') xlabel('X') ylabel('Y') 绘图如下：
t
z
1133
12
13
3
12112()()2,2,
1/82.45,1,
(0) 1.2,(0)0,(0)0,(0) 1.04935751
x x x y x r r y y
y x y r r r r x x y y μμμμμμμμμ+-=+--=-+-
-
==-======-
5．盐水的混合问题
一个圆柱形的容器，内装350升的均匀混合的盐水溶液。

如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。

开始时，容器内盐的含量为7千克。

求经过时间t 后容器内盐的含量。

做出如下假设：
1.假设在不同浓度的水中的盐扩散速度都相同。

2.假设任何时刻容器内的盐水都是均匀的。

3.用y （t ）表示容器内t 时刻的盐的含量，用W （t ）表示容器内t 时刻的水的总量，用O 表示盐水流出的速度，用I 代表纯水流入的速度，时间变化t ∆后容器内盐的含量为y （t+t ∆）。

考虑在t ∆内流出的盐水的为O t ∆则其流出的盐为
t O t W t y ∆)
()
(. 通过以上假设可以得如下模型：
t O t W t y t y t t y ∆-
=∆+)
()
()()( t
O I t W Odt t y t y dy t y )()()()()(0-+-
=+ 化简可得
t
O I t W O t y dt dy o )()()(-+-
= 使用MATLAB 求解：
syms x y z k;k=dsolve('Dy=-(y*Y)/(T(t0)+(C-Y)*t)','y(0)=7','t') 结果为：
(7*exp((O*log(W(t0)))/(I - O)))/exp((O*log(W(t0) + t*(I - O)))/(I - O)) 由题目可知W （t 0）=350，O=10.5，I=14，从而y(t)=7000000/(t + 100)^3。

6. 两种生物种群竞争模型
两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。

假设有甲、乙
两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic
规律，即有
⎪⎪⎩
⎪
⎪⎨⎧-=-=)1()1(2222211111N x x x N x x x λλ
其中x 1(t), x 2(t)分别为两种生物种群在时刻t 的数量，λ1，λ2分别为其自然增长率，N 1，N 2是它们各自的最大容量。

当两个种群在同一个自然环境下生存时，乙消耗的同一自然资源对甲的增长产生了阻滞作用，设为
甲对乙的阻滞作用设为
由于生物种群的数量很大，可视为时间t 的连续可微函数。

生物种群的相互竞争模型为
1）m 2(m 1)为种群乙(甲)占据甲(乙)的位置的数量，并且 m 2=αx 2; m 1=α2x 1。

计算x 1(t),
x 2(t), 画出图形及相轨迹图。

解释其解变化过程。

2）改变λ1，λ2, N 1, N 2, 而α1，α2不变，计算并分析结果；若α1=1.5，α2=0.7，再分析结果。

由此能得到什么结论。

用matlab 编程：
建立logistic.m 文件： function l=logistic(t,x) r=[1 1]; n=[100 100]; m=[1.5 0.7];
l=[r(1)*x(1)*(1-x(1)/n(1)+m(1)/n(2));r(2)*x(2)*(1-x(2)/n(2)+m(2)/n(1))]; 建立exp4-8.m 文件: x0=[10 10]; t0=[0 10];
[t x]=ode45('logistic',t0,x0);
plot(t,x(:,1),'r',t,x(:,2),'b');xlabel('时间');ylabel('种群数量'); title('种群数量与时间的关系');
text(t(16),x(16,1),'\leftarrow 种群一的数量');
text(t(18),x(18,2),'种群二的数\rightarrow','HorizontalAlignment','right'); pause
2
12x N α1
21x N α1211111221222221(1)(1)x x x x N N x x x x N N λαλα⎧
=--⎪⎪⎨⎪=--⎪⎩
由上图可知，开始阶段，两个种群的数量都随时间而增长；随着时间的推移，种群数量达到饱和。

两者数量基本持平，种群二略占优势。

当时间足够长时，两者都将接近环境的最大容量。

改变参数，有以下实验结果
1.改变种群的增长率
r1=1,r2=0.5
r1=0.5,r2=1
由图可知，增长率改变，图形也改变，但最增长趋势不变。

2. 改变环境容量
n 1=100,n 2=50
n 1=50,n 2=100
当两个物种对应的环境容量改变时，随着时间的推移，两者依旧将趋向于环境的最大
容量n 1,n 2。

3. 改变α、β的值
α=1.5、β=0.7,β<α时，种群一灭绝。

课堂练习：
7.试用微分方程数值解分析传染病模型SIR 的相轨线趋势。