捕食模型(生物数学)

捕食模型

（食饵捕食模型，生物数学重要模型）

假设及建立模型：

假设一个生态系统，其中含有两种生物 A生物和B生物，其中A生物是捕食者，B生物是被捕食者。建立捕食数学模型

1) 在观测数据（DATA1）无误差的情况下，确定模型中的参数，并分析误差。

2) 在观测资料有误差（时间变量不含有误差）的情况下，请分别利用观测数据DATA2和DATA3，确定参数在某种意义下的最优解，并与仿真结果比较，进而改进你们的数学模型。 3) 假设连观测资料的时间变量也含有误差，试利用数据DATA4，建立数学模型，确定参数在某种意义下的最优解。

通过对此生态系统的观测，可以得到相关的观测数据。观测数据的格式依次为：

观测时刻tj

、A生物数目x(tj)、B生物数目y(tj)

对于生态系统中的两种生物A和B，A生物为捕食者，B生物为被捕食者。在某

一段时期内，A生物的数量与B生物的数量之间存在一定的关系。

根据已知条件，可将（15）式改写为如下形式：

dx

=x(α1+α2y)dt

（1）

dy

=y(α3+α4x)dt

（2）

⎧x(t0)=α5⎨

⎩y(t0)=α6

其中αk(1≤k≤6)为模型的待定参数。

进行变换可得：

dyy(α3+α4x)=

dxx(α1+α2y)

（3）

即

积分得：

dy(α1+α2y)dx(α3+α4x)

=

yx

（4）

α（1lny-lny0)+α2(y-y0)+α3(lnx0-lnx)+α4(x0-x)=0

可将上述表达式改写成n元齐次线性方程组的形式，如下所示：

Am⨯nα=0

（5）

上述n元齐次线性方程组有非零解的充分必要条件是系数矩阵的秩R(A)

我们首先用DATA1中的3组数据确定a1，a2，a3，a4， 程序 clear

A=zeros(3,4);

A(1,1)=log(0.[**************]6 /60); A(1,2)= 60-0.[**************]6; A(1,3)=-log(11.[**************] /10); A(1,4)=10-11.[**************] ; A(2,1)=log(7.[**************]/60);

A(2,2)= 7.[**************]-60; A(2,3)=-log(3.[**************]6 /10); A(2,4)=10-3.[**************]6; A(3,1)=log(0.[**************]4/60); A(3,2)= 0.[**************]4-60; A(3,3)=-log(20.[1**********]798/10); A(3,4)=10-20.[1**********]798 ;

r=rank(A); % rank(A)=r

表1 a'k(1≤k≤4)的值

y=(y0+

ααα1ααα

lny0-4x0-3lnx0)-1lny+4x+3lnx （28）α2α2α2α2α2α2

ααααα1α

lny0-4x0-3lnx0)，β1=-1，β2=4，β3=3，

α2α2α2α2α2α2

如设：β0=(y0+

x1=lny，x2=x，x3=lnx，则（28）式可以写为如下形式；

y=β0+β1x1+β2x2+β3x3

（29）

对于（29）式中因变量y是自变量x={x1

x2

x3}的线性函数。可以建立起因变量

y的多元线性回归模型，

① 利用数据文件data2.txt,首先绘制出x(t)、y(t)与时间t的关系图象，和y(t)对x(t)的散点图，如图3所示（程序见gg1）。

利用MATLAB工具箱中的regress(y,x)命令，计算回归系数β的最小二乘估

ˆ及其置信区间，计算结果如表1所示，计算程序名为gg5。 计β

由它得到的模型为：

ˆ=-139.[1**********]43+19.[**************]1-9.[**************]x2y

+ （51）

99.[1**********]960

x3

结果分析：表1显示R2=0.[**************]8是指因变量y的99.71%可由模型（51）确定，F值远远超过F检验的临界值，p远小于α，因而模型（30） 从整体上看是可用的。表1的回归系数给出了模型（30）中β0，β1，β2，β3的 估计值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。

图3 y(t)对x(t)的散点图和x(t)、y(t)与时间t的关系图

得出β的估计值即可确定αi（i=1,2,3,4）之间的相互关系：α1=-β1α2，

ˆi=yi-yˆi(i=1, ,150)达到最小值进行最α3=β3α2，α4=β2α2，再利用残差ei=ε

优化计算来确定α2的取值，从而确定其它参数值，模型（30）的残差与x2的散

点图如图4所示。

表1 data2.txt中β的计算结果

从图中可以看出随着x2的增加，残差ei也逐渐增大。通过优化计算得出在捕食者A和被捕食者B共同存在，相互竞争达到动态平衡的状况下参数α2的最优值为α2=0.1，从而可以确定参数α1，α3，α4的最优值。

α1= -1.[**************]，α2=0.1，α3= 9.[**************]， α4

=

-0.[**************]

，

α5

= 12.[**************]，

α6=72.[1**********]777

利用上述数据对模型（30）的拟合效果如图5所示，相关程序为lorenzep.m和gg8.m。

从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30）式所示的模型时，假设回归变量x={x1

x2

x3}对因变量y的影响是相互独立

的。然而根据直觉和经验可以猜想，x1、x2、x3之间的交互作用会对y产生影响， 不妨简单地利用x1、x2、x3的乘积代表它们的交互作用，于是将模型（30）

增加一项得到：

⎧y=β0+β1x1+β2x2+β3x3+β4x1x2x3+ε

（52） ⎨2

ε N(0,σ)⎩

图4 残差ei与x2的散点图

图5 模型（30）的数据拟合

② 利用数据文件DATA3.txt,首先绘制出x(t)、y(t)与时间t的关系图象，和

y(t)对x(t)的散点图，如图6所示（程序见gg4）。

利用MATLAB工具箱中的regress(y,x)命令，计算回归系数β的最小二乘估

ˆ及其置信区间，计算结果如表2所示，计算程序名为gg6。

计β

图6 y(t)对x(t)的散点图和x(t)、y(t)与时间t的关系图

由它得到的模型为：

∧

y=-120.[1**********]12+17.[1**********]1x-8.[**************]x+

1

2

87.[1**********]46（53）

x3

结果分析：表2显示R2=0. [**************]5是指因变量y的92.29%可由模型（53）确定，F值远远超过F检验的临界值，p远小于α，因而模型（30）从整体上看是可用的。表2的回归系数给出了模型（30）中β0，β1，β2，β3的估计

表2 DATA3.txt中β的计算结果

值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。

得出β的估计值即可确定αi（i=1,2,3,4）之间的相互关系：α1=-β1α2，

ˆi=yi-yˆi(i=1, ,1500)达到最小值进行α3=β3α2，α4=β2α2，再利用残差ei=ε

最优化计算来确定α2的取值，从而确定其它参数值，模型（30）的残差与x2的散点图如图7所示。

从图中可以看出随着x2的增加，残差ei也逐渐增大。通过优化计算得出在捕食者A和被捕食者B共同存在，相互竞争达到动态平衡的状况下参数α2的最优值为α2=0.115，从而可以确定参数α1，α3，α4的最优值。

α1= -2.[**************]，α2=0.115，α3= 10.[1**********]203， α4

=-1.[**************]

，

α5

=12.[1**********]266，

α6=73.[1**********]651

利用上述数据对模型（30）的拟合效果如图8所示，相关程序为lorenzep1.m和gg10.m。

图7 残差ei与x2的散点图

从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30）式所示的模型时，假设回归变量x={x1

x2

x3}对因变量y的影响是相互独立

的。然而根据直觉和经验可以猜想，x1、x2、x3之间的交互作用会对y产生影响，不妨简单地利用x1、x2、x3的乘积代表它们的交互作用，于是将模型（30）增加三项得到：

⎧y=β0+β1x1+β2x2+β3x3+β4x1x2+β5x1x3+β6x2x3+β4x1x2x3+ε （54） ⎨2ε N(0,σ)⎩