1、捕食模型（食饵捕食模型，生物数学重要模型）假设及建立模型：假设一个生态系统，其中含有两种生物 A 生物和 B 生物，其中 A 生物是捕食者，B 生物是被捕食者。建立捕食数学模型1) 在观测数据（DATA1）无误差的情况下，确定模型中的参数，并分析误差。2) 在观测资料有误差（时间变量不含有误差）的情况下，请分别利用观测数据 DATA2 和DATA3，确定参数在某种意义下的最优解，并与仿真结果比较，进而改进你们的数学模型。3) 假设连观测资料的时间变量也含有误差，试利用数据 DATA4，建立数学模型，确定参数在某种意义下的最优解。通过对此生态系统的观测，可以得到相关的观测数据。观测数据的格式依次
2、为：观测时刻 、A 生物数目 、B 生物数目jt)(jtx)(jty对于生态系统中的两种生物 A 和 B，A 生物为捕食者， B 生物为被捕食者。在某一段时期内，A 生物的数量与 B 生物的数量之间存在一定的关系。 根据已知条件，可将（15）式改写为如下形式：（1） 12()dxyt（2） 34()yxdt056()xty其中 为模型的待定参数。进行变换可得：16k（3）3412()yxdx（4）3412()() xy即积分得： 10203040ln)()(ln)()yyxx（可将上述表达式改写成 n 元齐次线性方程组的形式，如下所示：mnA0（5）上述 n 元齐次线性方程组有非零解的充分必要
3、条件是系数矩阵的秩 R(A)n。我们首先用 DATA1 中的 3 组数据确定 ， 4321aa程序clearA=zeros(3,4);A(1,1)=log(0.1374480266382216 /60);A(1,2)= 60-0.1374480266382216;A(1,3)=-log(11.750840650304518 /10);A(1,4)=10-11.750840650304518 ;A(2,1)=log(7.108705996120129/60);A(2,2)= 7.108705996120129-60;A(2,3)=-log(3.4133367257849176 /10);A(2,
4、4)=10-3.4133367257849176;A(3,1)=log(0.4251595082899424/60);A(3,2)= 0.4251595082899424-60;A(3,3)=-log(20.80921881438798/10);A(3,4)=10-20.80921881438798 ;r=rank(A); % rank(A)=rn 时,该方程有无穷多个解,求它的一个基本解y=null(A, r )表 1 的值)4(ak1a2 3a4a-0.0478 -0.0042 -0.9925 0.11253 314140002222(lnln)llnyyxyx（28 ）如设： ， ， ，
5、31400022(ll) 1242， = ， = ， ，则（28）式可以写为如下形式；321xlny2x3ln0123yxx（29）对于（29）式中因变量 是自变量 的线性函数。可以建立起因变y123xx量 的多元线性回归模型，y 利用数据文件 data2.txt,首先绘制出 x(t)、y(t)与时间 t 的关系图象，和y(t)对 x(t)的散点图，如图 3 所示（程序见 gg1） 。利用 MATLAB 工具箱中的 regress(y,x)命令，计算回归系数 的最小二乘估计 及其置信区间，计算结果如表 1 所示，计算程序名为 gg5。由它得到的模型为：-139.4322355563643+19
6、.88421774692257 -y 1x9.988290331415989 + 99.80497072015960 2x 3x（51 ）结果分析：表 1 显示 =0.9971088966600568 是指因变量 的 99.71%2Ry可由模型（51 ）确定，F 值远远超过 F 检验的临界值，p 远小于 ，因而模型（30 ）从整体上看是可用的。表 1 的回归系数给出了模型（30）中 ， ， ，012的3估计值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。图 3 y(t)对 x(t)的散点图和 x(t)、y(t)与时间 t 的关系图 得出 的估计值即可确定 （i=1,2,3,4）之
7、间的相互关系： ，i12， ，再利用残差 达到最小值进行3242 (1,50)iiiey最优化计算来确定 的取值，从而确定其它参数值，模型（30）的残差与 的2x散点图如图 4 所示。表 1 data2.txt 中 的计算结果回归参数 参数估计值 参数置信区间0-139.4322355563643 -142.1292050800128 -136.7352660327159119.88421774692257 19.70905262672553 20.059382867119622-9.988290331415989 -10.20688059603675 -9.76970006679522639
8、9.80497072015960 97.65938535619841 101.9505560841208=0.9971088966600568 F=16784.58380705823 p0.00012R从图中可以看出随着 的增加，残差 也逐渐增大。通过优化计算得出在2xie捕食者 A 和被捕食者 B 共同存在，相互竞争达到动态平衡的状况下参数 的最2优值为 =0.1，从而可以确定参数 ， ， 的最优值。2134= -1.988421774692257， =0.1， = 9.980497072015961，1 2= -0.998829033141599， = 4 512.962285633035
9、274， =72.123015833447776利用上述数据对模型（30）的拟合效果如图 5 所示，相关程序为lorenzep.m 和 gg8.m。从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30 ）式所示的模型时，假设回归变量 对因变量 的影响是相123xxy互独立的。然而根据直觉和经验可以猜想， 、 、 之间的交互作用会对产生影响，y不妨简单地利用 、 、 的乘积代表它们的交互作用，于是将模型（30）1x23增加一项得到：012341232(0,)yxxN:（52）图 4 残差 与 的散点图ie2x图 5 模型（30）的数据拟合 利用数据文件 DATA3.txt,首先绘制
10、出 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 的关系图由它得到的模型为：-120.3100735257612+17.908648417141 -8.803113594738514 +y 1x2x87.8929067081046 3x（53 ）结果分析：表 2 显示 =0. 9228752184692365
11、是指因变量 的 92.29%可由2Ry模型（53）确定， F 值远远超过 F 检验的临界值，p 远小于 ，因而模型（30 ）从整体上看是可用的。表 2 的回归系数给出了模型（30）中 ， ，01， 的估计23表 2 DATA3.txt 中 的计算结果回归参数 参数估计值 参数置信区间0-120.3100735257612 -124.4170908611838-116.2030561903386117.9086484171410 17.64601818521920 18.171278649062792-8.803113594738514 -9.138961160791356 -8.4672660
12、28685673387.89290670810462 84.60654520450265 91.17926821170659=0. 9228752184692365 F= 5967.045867867118 p0.00012R值。检查它们的置信区间发现都不包含零点，表明回归变量都很显著。得出 的估计值即可确定 （i=1,2,3,4）之间的相互关系： ，i12， ，再利用残差 达到最小值进行3242 (1,50)iiiey最优化计算来确定 的取值，从而确定其它参数值，模型（30）的残差与 的2x散点图如图 7 所示。从图中可以看出随着 的增加，残差 也逐渐增大。通过优化计算得出在2xie捕食者
13、A 和被捕食者 B 共同存在，相互竞争达到动态平衡的状况下参数 的最2优值为 =0.115，从而可以确定参数 ， ， 的最优值。2134= -2.059494567971215， =0.115， = 10.10768427143203，1 23=-41.012358063394929， =12.82305785349266， =73.4072806908665561利用上述数据对模型（30）的拟合效果如图 8 所示，相关程序为lorenzep1.m 和 gg10.m。图 7 残差 与 的散点图ie2x从图像可以看出拟合的数据仍然有一些偏差，究其原因是我们在建立（30 ）式所示的模型时，假设回归变量 对因变量 的影响是相123xxy互独立的。然而根据直觉和经验可以猜想， 、 、 之间的交互作用会对产生影响，不妨简单地利用 、 、 的乘积代表它们的交互作用，于是将y1x23模型（30）增加三项得到：01234125136234123(0,)yxxxxN:（54 ）
