微分方程建模（人口预测，捕食者猎物）

陌上花开 ~

已于 2022-07-27 23:11:30 修改

阅读量2.7k

点赞数 1

分类专栏：数学建模文章标签：微分方程人口预测捕食者猎物模型

于 2022-01-25 02:36:37 首次发布

本文链接：https://blog.csdn.net/qq_45550375/article/details/119548703

版权

数学建模专栏收录该内容

14 篇文章 64 订阅

订阅专栏

解简单方程

解析解求法

解方程： $y - y^{'} = 2 x$

clear;clc

dsolve('y-Dy=2*x','x')  % 这里要指定自变量为x
% 2*x + C1*exp(x) + 2  (这里的C1表示任意常数，有时候也会出现C2 C3等)

dsolve('y-Dy=2*x')  % 如果不指定自变量的话，会默认自变量为t，x会看成一个常数
% 2*x + C2*exp(t)

% 注意：最新版本的matlab会逐渐淘汰上面那种写法
% 下面这种写法是新版的matlab推荐的方式（和符号运算中解方程的写法类似）

syms y(x)
eqn = (y - diff(y,x) == 2*x);    % 注意原来方程中的“=”一定要改成“==”
dsolve(eqn)

解方程： $y - y^{'} = a x$

%% 如果微分方程中还有其他的未知参数

% 方法1
dsolve('y-Dy=a*x','x')  % a是一个未知的参数

% 方法2
syms y(x) a
eqn = (y - diff(y,x) == a*x);  
dsolve(eqn)

结果：

a + a*x - C1*exp(x)

解方程： $y-y'=2x,y\left( 0\right) =3$

% 方法1
dsolve('y-Dy=2*x','y(0)=3','x')
% 2*x + exp(x) + 2

% 方法2
syms y(x)
eqn = (y - diff(y,x) == 2*x);  
cond = (y(0) == 3);
dsolve(eqn,cond)
% 2*x + exp(x) + 2

人口预测模型

Malthus模型

求解：

clear;clc
x = dsolve('Dx=r*x','x(0)=x0','t')    % x = dsolve('Dx=r*x','x(t0)=x0','t')
% x = x0*exp(r*t)

% 初始人口为1000，画出50年且增长率分别为0.5%，1%，1.5% 一直到5%的人口变化曲线
x0 = 1000;
for i = 1:10
    r = 0.005*i;
    xx = subs(x); % 把上面这个式子中的x0和r替换成确定的值
    fplot(xx,[0,50])   % 绘制表达式的图形
    hold on
end

结果：
在这里插入图片描述

logistic模型

求解：

clear;clc
x = dsolve('Dx=r*(1-x/xm)*x','x(t0)=x0','t')  

% x = -xm/(exp(log(1 - xm/x0) - r*t + r*t0) - 1) 
% mupad  % 把计算出来的结果粘贴过去可以得到直观的表达式
% 高版本可以使用实时脚本

t0 = 0;x0 = 1000;
xm = 10000;
r = 0.05;
xx = subs(x);    %  10000/(exp(log(9) - t/20) + 1)
fplot(xx,[0,200])  

% 也可以这样画图：
t = 0:200;
x = 10000 ./ (exp(log(9) - t/20) + 1);
plot(t,x,'-')

用实时编辑器的效果(实时脚本)：
在这里插入图片描述

结果：
在这里插入图片描述

捕食者猎物模型

地中海鲨鱼问题

考虑人工捕获，食饵（即猎物）和捕食者在时刻 $t$ 的数量分别记作 $x_{1}(t)$ 和 $x_{2}(t)$ ， $\lambda$ 为捕食者捕食能力， $e$ 为人工捕食强度， $r_{1}$ 为食饵的自然增长率， $r_{2}$ 为捕食者的死亡率，有：
$\begin{cases}\dfrac{dx_{1}}{dt}=x_{1}\left[ \left( r_{1}-e\right) -\lambda _{1}x_{2}\right] \\ \dfrac{dx_{2}}{dt}=x_{2}\left[ -\left( r_{2}+e\right) +\lambda _{2}x_{1}\right] \\ x_{1}\left( 0\right) =c_{1},x_{2}\left( 0\right) =c_{2}\end{cases}$

我们尝试用Matlab求解析解，发现求不出来：

dsolve('Dx1=x1*(0.9-0.1*x2)','Dx2=x2*(-0.6+0.02*x1)','x1(0)=25,x2(0)=2','t')

结果：

[ empty sym ]

用ode45函数求数值解：

% 自变量t的范围为0-15年，食饵和捕食者（鲨鱼）初始值分别为25和2 
[t1,x1]=ode45('pre_war',[0 15],[25 2]);  % 战前的微分方程
[t2,x2]=ode45('post_war',[0 15],[25 2]); % 战后的微分方程

pre_war.m：

% 战前
function dx=pre_war(t,x)
    dx=zeros(2,1); 
    dx(1)=x(1)*(0.7-0.1*x(2));
    dx(2)=x(2)*(-0.8+0.02*x(1));
end

post_war.m：

% 战后
function dx=post_war(t,x)
    dx=zeros(2,1); 
    dx(1)=x(1)*(0.9-0.1*x(2));
    dx(2)=x(2)*(-0.6+0.02*x(1));
end

绘制鲨鱼和食饵数量变化的相轨线图：

figure(1)
plot(pre_prey,pre_shark,'--r',post_prey,post_shark,'-b')
legend('战前','战后')
title('鲨鱼和食饵数量变化的相轨线图')
xlabel('食饵数量');   ylabel('鲨鱼数量')

在这里插入图片描述
绘制捕食者和猎物战前战后的时间-数量图：

figure(2)
plot(t1,pre_prey,'--r',t1,pre_shark,'-r',t2,post_prey,'--b',t2,post_shark,'-b')
legend('战前食饵数量','战前鲨鱼数量','战后食饵数量','战后鲨鱼数量')
xlabel('时间');   ylabel('数量')

在这里插入图片描述
绘制战前战后鲨鱼的鲨鱼比例-时间图（鲨鱼比例 = 鲨鱼数量 /（鲨鱼数+食饵数））：

% 鲨鱼比例  =  鲨鱼数量 /（鲨鱼数+食饵数）
pre_rate=pre_shark./(pre_prey+pre_shark);         % 战前的鲨鱼比例    
post_rate=post_shark./(post_prey+post_shark);   % 战后的鲨鱼比例
figure(3)
plot(t1,pre_rate,'--r',t2,post_rate,'-b')
legend('战前的鲨鱼比例','战后的鲨鱼比例')