数学建模——常偏微分方程模型（下）（模型详解和matlab代码）

5.竞争捕食模型

 

竞争模型

问题：一个小池塘，它的环境足以维持一些野生动物的存活，池塘只有两种鱼，鳟鱼和鲈鱼，问两种鱼能否在池塘中共存？

 

假设：1.环境能支持无限条鱼（单一鱼种能够无限制进行繁殖）

      2.两种鱼一同竞争生存空间和公共食饵，会降低双方的增长率

建模：

令

先考虑环境只有鳟鱼，由于环境能支持无限条鳟鱼，得到方程

考虑鳟鱼和鲈鱼共存，鲈鱼会降低鳟鱼的增长率，且鲈鱼、鳟鱼的数量越多，则竞争越激烈，所以降低增长率的一项与鳟鱼、鲈鱼数量都有关

得到鳟鱼的方程

同理得到鲈鱼的方程

模型分析

上两式为自治微分方程,求解下列方程组可以得到微分方程组的静止点

 

求得静止点为

我们可以通过使用MATLAB画出相图查看静止点是否稳定

MATLAB画相图的程序如下：

clear;clc

syms a b m n;

m=2;

n=2;

a=2;

b=2;

[x,y]=meshgrid(linspace(0,5));

streamslice(x,y,a*x-b*x.*y,m*y-n*x.*y);

xlabel('x');ylabel('y');

 

结果如下：

可以由图观察得出只有初始点在 上的时候，才会去到平衡点

 

上面这个模型是假设资源无限的，

若资源是有限的

则为下面进行了改进的方程组

X、Y为鳟鱼和鲈鱼的数量上限，c,q的意义为鲈鱼会也会以c倍消耗鳟鱼的食物，鳟鱼会以q倍消耗鲈鱼的食物

 

 

捕食模型

问题：鲸鱼捕食磷虾，磷虾以海中的浮游生物为食

建模：

令磷虾的数量为x，鲸鱼的数量为y

当只有磷虾单一物种，假设海洋能够支持无限量的磷虾

当有鲸鱼的时候，捕食程度由鲸鱼和磷虾接触的程度成正比，接触程度与x，y都有关

对鲸鱼有方程

 

 

6.战争模型

正规战模型

(i)双方士兵分别为

每次射击的命中率

(ii)双方的非战斗减员率仅与本方兵力成正比，减员率系数为

(iii)双方的兵力增援为

 

通过上面三个假设我们列出微分方程组

 

若不考虑非战斗减员和增援，且初始兵力为

MATALB程序：

syms x y t c b x0 y0

[x,y]=dsolve('Dx=-a*y','Dy=-b*x','x(0)=x0','y(0)=y0')

结果为：

x =(exp(-t*(a*b 1eee0 )^(1/2))*(a*b)^(1/2)*((b*x0 + y0*(a*b)^(1/2))/(2*(a*b)^(1/2)) + (exp(2*t*(a*b)^(1/2))*(b*x0 - y0*(a*b)^(1/2)))/(2*(a*b)^(1/2))))/b

y =exp(-t*(a*b)^(1/2))*((b*x0 + y0*(a*b)^(1/2))/(2*(a*b)^(1/2)) - (exp(2*t*(a*b)^(1/2))*(b*x0 - y0*(a*b)^(1/2)))/(2*(a*b)^(1/2)))

 

 

游击战模型

 

若没有非战斗减员和增援

使用MATLAB进行求解：

使用dsolve函数提示没有通解，那么只能使用solve函数解数值解

 

混合战模型(x方为游击部队，y方为正规部队)

双方都没有非战斗减员和增援

使用MATLAB进行求解

syms x y t c b x0 y0

[x,y]=dsolve('Dx=-c*x*y','Dy=-b*x','x(0)=x0','y(0)=y0')

结果为无解，那么只能使用solve函数解出数值解

 

 

7.军备竞赛经济模型

若国家1没有其它国家的威胁，那么有理由假设防御支出额按已经支出额的比例减少

若有国家2的威胁，那么假设防御支出按国家2的防御支出 的比例增加，比例为  

但从现实情况出发，b不应该是一个常数，例如，在国家1的实力本来就比国家2强的情况下，b可能就比较小，在国家1的实力本来就比国家2弱的情况下，b可能就比较大

考虑到现实情况应该是，即使双方的防御支出都为0，那么国家1仍然有必要进行防御支出，以防范国家2的潜在威胁，因此应该加一个常数项 ，以反应国家1对国家2感到的所有潜在的不安的因素

根据上面的分析，可以对国家1和国家2列出下列微分方程组

 

8.传染病模型

模型的假设：

1. 人群有三类：易感染者（从未被感染）、已感染者（能传播流感）、移出者(包括死亡者,不能再次被感染)，每个人都是上面这三类之一
2. 群体封闭，没有人进入或外出
3. 该年一旦染上流感，不会再次染上
4. 流感疾病平均持续期为53 周，此期间病人为已感染者

 

模型的建立：

由于平均持续期为 周，已感染人群里面感染时间超过 周的人数占比为( - )/ = =0.6,,因此一周之后已感人群中有60%的人转为移出者

单位时间内易感人群减少的数量与已感人群数量和易感人群数量有关，已感人群数量和易感人群数量越大，则单位时间内易感人群增加的数量越大

其中a的估计值为a=0.001407

 

综合上式会有下列微分方程组，假设初始人群为1000人，已感人群为5人，易感人群为995人，移出人群为0人

对于上面的微分方程组，我们可以使用ode45来进行求解，MATLAB代码如下：

F.m文件代码

function dy=F(t,y);

dy=[0.6*y(2);-0.6*y(2)+0.001407*y(2)*y(3);-0.001407*y(3)*y(2)];

%y1代表Rt,y2代表It,y3代表S(t)

 

运行的主程序：

[T,Y]=ode45('F',[0,200],[0;5;995])

plot(Y)

 

结果如下：红色为易感人群，蓝色为移出人群，绿色为已感人群

 

 

9.圆桶下沉模型

圆桶质量m=239.46Kg，体积V=0.2058 ，海水密度ρ=1035.71Kg/m3,水深s=90m，水的阻力f和速度大小v成正比例，比例常数为k=0.6,如果圆桶速度小于12.2m/s 就说明方法是安全可靠的。

1.判断是否安全

对圆桶进行受力分析，有重力G，有浮力H，有阻力

G=mg

H=

因此会有

其中

求解方程组的MATLAB程序

clc,clear

syms m V rho g k

s=dsolve('m*D2s-m*g+rho*g*V+k*Ds','s(0)=0,Ds(0)=0');

s=subs(s,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});

s=vpa(s,10)   %vpa设置精度，有效数字为10位

v=dsolve('m*Dv-m*g+rho*g*V+k*v','v(0)=0');

v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});

v=vpa(v,7)

y=s-90;

tt=solve(y)

vv=subs(v,tt)

 

求解得到

当 ,所以不安全

 

2.一般情况下 ,重新计算

 

MATLAB程序如下：

clc,clear

syms m V rho g k t real

v=dsolve('m*Dv-m*g+rho*g*V+k*v^2','v(0)=0');

v=subs(v,{m,V,rho,g,k},{239.46,0.2058,1035.71,9.8,0.6});

v=vpa(v,7)

v1=real(v),v2=imag(v)  %求速度函数的实部和虚部

tt=solve(v1-12.2)  %求临界的时间

s=int(v1,0,tt)    %使用积分求临界的位移