猎食者-猎物模型的求解问题

21-22-2学期|工程计算期末报告

一、问题描述

捕食者-猎物模型(predator-prey models).又称寄生物-寄主模型，是表达捕食者-猎物系统内种群数量变化动态的数学方程。最简单的捕食者-猎物模型如下：

$Missing or unrecognized delimiter for \left \left{ \begin{array} { l l } {x^{\prime}(t)=a x-bxy}&{b>0}\{y^{\prime}(t)=-cy+dxy}&{d>0}\{x(0)=x_{0},y(0)=y_{0}}\\end{array}\right.$

其中，x 和y 分别为猎物和捕食者的数量，a 为猎物的自然生长率，c 为捕食者的死亡率。这里认为捕食者仅以猎物为食，且猎物的食物足够多。
科研人员通过 2 年的观察获得下表所示的捕食者和猎物的数量变化情况，请根据下表数据求得模型中的各项参数值，并对模型进行求解，得到这 24 个月内以天为单位的捕食者和猎物种群数量变化数据。
注：用符号运算求解析解时仅保留解的整数部分。

二、数学建模

由题意知，本题并非一般的微分方程组求解，给出的模型中含未知量参数，因此需要根据表中的数据通过估计a,b,c,d的取值。考虑到求解的捕食者-猎物模型不需要表中数据与模型完全吻合，故选择拟合的方法。
整体思路是根据原模型构造新的合适的模型，带入表格数据后利用最小二乘原理将型值进行分析运算，以求出最接近或最能表示数据趋势的“函数”，即确定参数a,b,c,d的值。四个参数的值确定下来后，之后便和正常的常微分方程组求解问题一致，就容易解决了。
所以总结一下，该题主要考察了拟合和常微分方程组求解两部分的内容。

具体分析过程见下：
将原微分方程的模型做如下变换,方便后续积分。
$$
\begin{equation}
\left{\right.
\end{equation}$$

对上式在相邻两时间点上，即上积分，可得：
$$
\begin{equation}
\left{\right.
\end{equation}$$

根据梯形公式，可将上面方程组最右边两项积分符号消去，即有
$$
\begin{equation}
\left{\right.
\end{equation}
$$

为使表示简洁，记代入（2）式，可得
$$
\begin{equation}
\left{\right.
\end{equation}$$

将上式中含i各项写开，再化为方程组的形式，即下面式：
$$
\begin{equation}
\left{\right.
\end{equation}$$

其中，分别为：

即为要求的解，根据最小二乘解原理：，代入表格数据后可通过matlab求解
常微分方程组确认后，这里再通过matlab使用经典四阶龙格库塔求解方程组的数值解，并绘制数值解图像以及相平面图像，具体见下方实现代码及结果分析。

三、实现代码

function finalWork
%题中数据导入
x0=[25,2]; 
T=1:24;
X=[25 55 97 55 10 3 2 2 3 7 14 30 65 100 37 7 3 2 2 4 8 17 37 76];
Y=[2 3 7 25 26 17 11 7 4 3 2 2 3 10 28 24 16 10 6 4 3 2 2 4]; 
%微分方程组求解
[t,x]=Classical_RK4s(@(t,x)Least_squares_BuShiLieWu(x,T,X,Y),[0,50],x0,0.1); 
%微分方程数值解图像以及相平面曲线图像的绘制
h=plot(t,x,'LineWidth',2);
grid on 
xlabel('t')
ylabel('x/y')
title('微分方程数值解')
set(h(2),'LineStyle','-.')
H=legend('{\itx}({\itt})','{\ity}({\itt})');
set(H,'fontname','consolas','fontsize',12)
figure 
plot(x(:,1),x(:,2))
xlabel('x')
ylabel('y')
title('相平面曲线') 
function xd=Least_squares_BuShiLieWu(x,T,X,Y)
% 最小二乘原理拟合参数值，并求出常微分方程组
% 输入参数：
%      ---x：微分方程组中的自变量
%      ---T：时间范围区间
%      ---X：表格中猎物
%      ---Y：表格中捕食者数量
% 输出参数：
%      ---xd：返回的不含参数的常微分方程组
T=T(:);X=X(:);Y=Y(:);
disp(x)
S2=diff(T).*(X(1:end-1)+X(2:end));
S1=diff(T).*(Y(1:end-1)+Y(2:end));
A1=[diff(T) -S1];A2=[-diff(T) S2];
B1=[log(X(2:end)./X(1:end-1))];B2=[log(Y(2:end)./Y(1:end-1))]; 
T1=A1\B1;
T2=A2\B2;
a=T1(1);b=T1(2);c=T2(1);d=T2(2);
xd=[(a-b*x(2)).*x(1);(-c+d*x(1)).*x(2)]; 
function [x,y]=Classical_RK4s(odefun,xspan,y0,h,varargin)
% 经典Runge-Kutta法求解常微分方程组
% 输入参数：
%      ---odefun：微分方程的函数描述
%      ---xspan：求解区间[x0,xn]
%      ---y0：初始条件
%      ---h：迭代步长
%      ---p1,p2,…：odefun函数的附加参数
% 输出参数：
%      ---x：返回的节点，即x=xspan(1):h:xspan(2)
%      ---y：微分方程的数值解
x=xspan(1):h:xspan(2);
y(:,1)=y0(:);
for k=1:length(x)-1
    K1=feval(odefun,x(k),y(:,k),varargin{:});
    K2=feval(odefun,x(k)+h/2,y(:,k)+h/2*K1,varargin{:});
    K3=feval(odefun,x(k)+h/2,y(:,k)+h/2*K2,varargin{:});
    K4=feval(odefun,x(k)+h,y(:,k)+h*K3,varargin{:});
    y(:,k+1)=y(:,k)+h/6*(K1+2*K2+2*K3+K4);
end
x=x';y=y';

四、结果分析

在函数Least_squares_BuShiLieWu中求得
所以原常微分方程组为：
$Missing or unrecognized delimiter for \left \left{ \begin{array} { l l } {x^{\prime}(t)=1.0340 x-0.0520xy}\{y^{\prime}(t)=-0.5310y+0.0106xy}\{x(0)=25,y(0)=2}\\end{array}\right.$
接下来通过函数Classical_RK4s，求出微分方程组的数值解，并绘制相关图像，如下：

分析图象可知，在猎物数量达到顶峰后，捕食者数量开始迅速增加，伴随猎物数量的下降；在猎物数量先下降到一定程度，捕食者数量也开始下降，待到捕食者数量下降到一定程度，猎物数量又开始增加。由相平面曲线闭合可推断出，捕食者-猎物的数量变化是一个循环往返的过程，这个结果也符合已知的生态学规律。
关于算法的时间复杂度分析，通过重复实验发现其主要与四阶龙格库塔的迭代步长有关，通过将迭代步长依次设为0.1，0.01，0.001等等后发现，数值解的图像基本不发生改变，但程序所需时间随步长的减小显著提升，所以这里选择h=0.1以获得更好的时间效果。