鱼类物种的数量模型

1、鱼类物种的数量模型摘要本文针对种鱼以及其捕食者种鱼各自随时间的数量变化， 估计出种鱼和种鱼之间的相互影响， 在此基础上， 依据问题中给出的微分方程组， 进行数据拟合以及修改，从而最终得出两物种的数量模型。针对问题一， 根据已给数据， 运用绘制图形， 发现种鱼的数量符合阻滞增长模型（模型），由此建立与的数学关系为，于是运用拟合，计算出种鱼的自然增长率以及环境资源容许的种群最大数量，总体拟合程度较高。针对问题二，将 a,b, c, d 数值以及初始值带入，运用计算出微分方程的数值解 x(t) 和 y(t) 并作图。针对问题三，根据已知数据，首先分别绘制x(t) 及 y(t ) 的图像，再将问题二中

2、 a,b,c, d 数值带入，按照问题二的方法计算出两组数据 x(t) , y(t ) , 并与原数据比较，发现图形走势基本一致， 再根据数据初步拟合的结果显示， 拟合后的图像与原图像横坐标相差倍，故通过遍历 () 邻域内的值优化求解。然后通过控制变量，采用最小二乘拟合法循环逐个拟合 a, b,c, d 四个参数，用平方和残差值的最小值来评价参数的优劣程度， 从而筛选出最优解， 并代入微分方程组，将再次计算出的数据与原数据比较， 发现拟合程度比较理想。 拟合过程中还采用微元的思想，将整个时间区间分成无限个小区间进行积分求和运算，得到 a,b,c, d 的数值解为 ，最后根据 a, b, c,

3、d 求出两个初值分别为。针对问题四，通过观察问题三绘制出的拟合图像， 发现计算出的数据比原数据略大，说明原微分方程组是理想模型， 没有考虑起物种自身阻滞作用的项， 于是在原有的微分方程基础上增加了项，提高了数据的拟合程度。针对问题五，分析了微分方程组中的各参数代表的含义，得出a, b, c, d 不依赖于数据的结论， 即无论 y 取何种初值， 图像的形状不变， 只是在原有基础上进行平移，基于这个结论， a, b, c, d 的数值可解，但对于 y 的初值无法给定。本文亮点首先在于问题三中的数据拟合过程， 通过对拟定的未优化数据与原数据各自形成的图像特征的比较，确定 a,b,c,d 准确值存在的

4、大致范围，运用最小二乘拟合法逐个拟合优化， 最终取得最优解。 然后在问题五中， 清楚地分析了 a, b, c, d 与 x(t ) 、 y(t ) 的关系中最重要的一点就是和各时刻对应的数值无关，而是和图像的形状、周期有关。关键词：阻滞作用、最小二乘拟合、一、问题重述在某水域生长的一种鱼类（A 种鱼），该水域出现了一种以A 种鱼为食的另一种鱼类（ B 种鱼），记 t 时刻 A 种鱼的数量为 x(t ) ，B种鱼的数量为 y(t) 。A 种鱼与 B 种鱼的数量相互影响，本题要求通过建立数学模型解决以下五个问题：1 根据表格中记录的一段时间内所观测到的A种鱼在若干观测观测的数量，描述 A 种鱼的数

5、量随时间变化的情况。2 如果这两种鱼的数量变化符合以下的微分方程组：dxbxyaxdtdydxycydt根据题中给出的及初始值，计算t=10,20, ,100 时两种鱼的数量。3 根据在一段时间内观测到的两种鱼的数量，如果两种鱼的数量变化符合问题二中的微分方程组，求其中的常数a,b,c, d 以及初始时刻两种鱼的数量x(0), y(0) 。4 问题二中的微分方程组能否很好的反映问题三数据？请你修改微分方程组的形式，使得改进后的微分方程组更好的反映问题三中的数据。5 如果在问题三中，无法观测到 B 种鱼的数据（ y(t ) ），只能观测到 A 种鱼的数据（ x(t ) ），是否可以求解与问题三与

6、问题四相应的问题。二、模型假设1. A 种鱼的数量变化符合阻滞增长模型；2. 除了食饵、捕食者关系以及自身阻滞作用， 其他因素对两种鱼的数量无影响；三、符号说明N：环境资源所能容纳的最大种群数量；r ：种群数量的自然增长率；： A 种鱼在初始时刻的种群数量；： B 种鱼在初始时刻的种群数量；：平方和残差值的最小值；：表示的是单位数量 A 种鱼提供的供养 B 种鱼的食物量为单位数量的 B 种鱼消耗的供养 B 的食物量的倍；：表示单位数量 B 种鱼消耗 A 种鱼作为食物量为单位数量的 A 种鱼自身减少数量的倍。四、模型的建立与求解4.1 问题一4.1.1 问题一的分析当水域中只有 A 种鱼时便构成

7、了单种群模型，为了清楚的认识问题一中 A 种鱼的数量变化规律，根据 A 种鱼在从 0 到 19 之间不同时刻的数量，首先运用绘制出 A 种鱼数量随时间变化的图像如下（如图 4.1.1-1 ），它大体上描述了该种群数量随时间变化的规律， 从而判断出 A 种鱼符合的模型， 并根据数据拟合出模型中的未知参数。图4.1.1-14.1.2 模型的建立根据图 4.1.1-1，不难看出A 种鱼的增长规律符合阻滞增长模型（logistic模型）。将 t 时刻的种群数量记为x(t) ，并将 x(t ) 视为连续、可微函数，记初始时刻的数量为 ,假设 A 种鱼的自然增长率为常数r，可以注意到 r 随着数量 x(t

8、) 的增加而下降，即表现为阻滞越来越大。若将r 表示为x 的函数，则它是减函数，有，对作一个简单的假定，为此，引进环境资源所能容纳的最大种群数量为N，设带入上述公式，得上述式右端的因子体现了种群数量自身的增长趋势， 因子则体现了种群数量增长的阻滞作用，式构成阻滞增长模型。 14.1.3 模型的求解问题一要求 A 种鱼的数量变化，应用分析出的规律可得等式其中 =25.81 ，运用拟合所给的数据得，拟合图像如图 4.1.3-1 ：图 4.1.3-1从图 4.1.3-1 中可以看到， 它不仅大致描述了 A 种鱼随时间的变化情况，而且更好地与真实图像相拟合。4.1.4 模型的结果分析由图像的拟合程度可

9、以直观地判断出，当水域中只有 A 种鱼时，该种群的数量变化满足模型，即自身增长阻滞因素的存在是肯定的。4.2 问题二4.2.1 问题二的分析问题二中出现了以A 种鱼为食的B 种鱼，要求通过问题中给出的微分方程组，相关系数的值以及两种鱼数量的初始值计算出两种鱼在不同时刻的数量， 这属于中微分方程的初值解问题，故直接根据其中数据编写程序求解即可。4.2.2 模型的求解dx已知微分方程组的形式为dtaxbxydycydxydt，根据问题二中的数据， a=2，b=0.01，c=1，d=0.01，另在初始时刻 t=0 时，有。求解如下 (精确到小数点后 13 位)：106.695610622932092

10、42.8269687425382043.4135062005784555.56049274398930402.368367290486124.6846053281384096.307582411949135.17813481510405019.302277240457459.7529087274689606.50808159252891150.1411739936237012.2164610147701396.89782586901080322.124619684216432.70427018826290182.57493739181940.505424275731910035.1046386

11、31399344.5736463355181并利用作出如下图像：图 4.2.2-14.2.3 模型的结果分析图 4.2.2-1 只选取了 10 组数据进行绘制，由图像的大概走势可以估计出 A 种鱼和 B 种鱼数量变化关系满足所给定的微分方程模型的具体表现。4.3 问题三4.3.1 问题三的分析首先，根据问题三的数据作出数量变化的直观图像，如图4.3.1-1 ：图 4.3.1-1可以观察到， 该图像呈平滑的周期性曲线，且不存在明显的坏点， 可以认为各点均较好地坐落在函数曲线上。另外，容易发现，图 4.3.1-1 与图 4.2.2-1 走势相近，但是数据初步拟合的结果显示两者横坐标密度相差 10

12、倍，因此，记问题二中参数为 (), 考虑将问题二中缩小 10 倍后得() ，以时的数值为初值，利用微元法，将分成无数个极小的区间，并在每个小区间内对微分方程进行积分处理， 通过编程还原拟合后的数值， 即时的数值， 所得数据结果如下 ( 此处只列出前 20 组数据，精确到小数点后 13 位) ：tx(t)y(t)145.790000000000041.4000000000000253.716356967720439.3666985966953363.128614232103137.7610995249215474.294928452374836.5926355735163587.52080460

13、8288735.89010074330116103.14545976202135.71215140919517121.54182305817336.14730632021508143.06106264111237.35879432925819168.16799934750839.495438265894510197.21330086710842.843281425403811230.25693548283147.950682837306512267.06817207076655.628317265363613306.90550068617166.978236135377514347.80197

14、809090384.121052912449215385.749575134144110.17318125515616414.364741995805148.96946198322017424.887801816950205.06275076590118407.479644798270282.04770131633919358.590063197010375.25628833583720286.887513576624469.433817327986另外，利用问题二中的参数得到拟合曲线，如图4.3.1-2：图 4.3.1-2作出关于 x(t ) 比较图像，如图4.3.1-3 ：图 4.3.1-

15、3作出关于 y(t) 的比较图像，如图4.3.1-4 ：图 4.3.1-4观察图像 , 可以判断以 () 为参数拟合的方程所得到的数据与问题三的数据有着较高的拟合度，因此可以在 () 各参数的相应邻域内考察更优解。4.3.2 模型的求解根据以上分析，可以通过遍历 () 邻域内的值求得最优解，由于分析中认为问题三数据较好地坐落在函数曲线上， 因此，任意数据作为初值计算而得的曲线具有相同的参数与相似的表达式，仅相位有差异，不妨选取处的值 , 计算相应参数值对应的时的值， 并规定用平方和残差值的最小值来评价参数的优劣程度， 其中 n=90 为拟合数据长度。利用以下式子：以及 MATLAB软件进行最小

16、二乘拟合 2 。若对于各遍历 n 个数据，则时间复杂度为，超出了可接受范围，考虑到已经拥有拟合度比较高的 () 为，不妨在考虑 1 个参数时令其余个不变， 得到较优参数后立即替代原有参数， 并逐个多次进行拟合修正， 令拟合结果不断逼近真实最优结果，则时间复杂度将为。利用上述算法，求得较优参数，并得到优化参数拟合序列中的90 组数据如下（精确到小数点后13 位）：tx(t )y(t)tx(t)y(t )1145.790000000000041.40000000000005615.138861457279969.00309568748901254.360795919177839.378646675

17、04955717.422343520928163.45918163106331364.672210450977637.81424460432735820.169980662240758.50728347922591477.066277838805836.72247234591745923.482194752630154.10524540527241591.933719213754536.14393592777826027.481263677354550.217510783100216109.70421118358836.15935520908276132.311321335517846.815

18、122187207917130.84638583274336.88956388021856238.138357640069943.875721392695618155.77628917853738.55466514941876345.136713722011241.402188592139819185.06346638619441.36381849196196453.531307510573039.402859694033220219.01922459118845.73290922348496563.661063058341337.831650334543821257.399150557887

19、52.45299846131746675.890794672368136.687433354563922299.41248660768962.67737517886776790.607178035897636.018749193867223343.45032550993277.868865866179368108.21875020836735.923805891108124384.385849491488101.69953278647569129.15590962540436.550479083822025414.728642761561138.09323862032570153.870916

20、09883238.096312008425426425.268162794879190.41725398555471182.83789081666440.808515500215627405.249485644033261.34831600190172216.33357034984144.929902318401728351.024016487800346.05565198239273254.66464046030551.204413980456329274.641783901213428.36968861623874297.08027467818861.165672455961830196.

21、751522338243490.10805128961375340.79986744442876.907637653628231132.941595890738522.40912493943776381.009198336474101.0833417589683288.0038873030939527.20529120827677410.860432068282136.9048892342973358.6034216805480512.85359084402478421.472118490313188.1434568187323439.8453802685919487.252409565602

22、79402.615608861006257.7546085486483528.1860215083149455.95261725743480350.596735684413340.7529514802883620.6983113177825422.67457147855181276.336736263480421.9499857999543715.8307206710633389.44584580155882199.352233884202483.9504060728033812.6053131914126357.39485067349183133.700452851551518.753781

23、9472613910.4113961448450327.1058997676898488.2828551636930524.357614303775408.90954296930337298.8399933572798558.9715920647978510.263630920899417.88108718175481272.6731142710608640.2082156612421484.976897532318427.18653888312947248.5816280925678728.5160333628630453.981513400348436.73431499277001226.

24、4932146991138820.9848902487768420.979015245146446.47277107575796206.2949568138268916.0778258606198387.988162216631456.36771240761249187.8629300749419012.8230217278798356.140163089553466.39894786128350171.0714301236519110.6064174083346326.024374711171476.55515937956694155.796320241832929.086609598069

25、96297.907333751208486.83299283564912141.915797341273938.04579547838872271.868764055616497.23594100965142129.313878280482947.34319024571603247.888308198598507.77381022211353117.885882250924956.88671058645735225.896141750932518.46222428117785107.534499072907966.62359151767271205.782313839817529.322472

26、1691988198.1691962642994976.51954866447124187.4243467651255310.381508042742989.7062190405298986.55476218018557170.6979425876145411.671951232588582.0685903145883996.71798829256552155.4805915669265513.240968613729775.1877806758082007.00555086926924141.652416308371同时经过程序计算， t 从 0 时刻到 10 时刻对应的 x（t ）与 y（

27、t ）的数值如下：tx(t)y(t)tx(t)y(t)08.68581100.015618.61259.241919.6501591.3295721.660754.6921210.828383.4872825.338950.6613312.259576.4178929.776647.1178413.991970.05731035.130744.0366516.085264.3489从中也可以得到x(t) 与 y(t) 的初始值分别为 =8.68581，=100.0154.3.3 模型的可信度检验根据上述得到的数据， 可以得到以下通过优化的拟合参数图，如图 4.3.2-1:图 4.3.2-1优化

28、拟合后的 x(t) 与原数据的比较图，如图4.3.2-2 ：图 4.3.2-2优化拟合后的 y(t ) 与原数据的比较图，如图4.3.2-3 ：图 4.3.2-3同时两组数据的相轨线图像也显示了较为理想的拟合程度，如图4.3.2-4 ：图 4.3.2-4借助这些图像可以直观地表达出， 拟合后的参数使得计算结果与所给数据更加吻合。这也进一步证明了模型 的可信度比较高， 所得的参数值及变量初值与真实值有着较高的拟合度。同时由于相轨线是封闭曲线，可以确定的是和 y(t)均满足周期函数的特性。4.4 问题四观察问题三中的图像，发现拟合图像总是比真实图像高出一些，结合问题一中采用的人口环境阻滞增长模型，

29、 不难发现，任务问题二中的模型缺少了环境阻滞参数。因此提出以下的优化模型 IV ：（4.4.1 ）类比于问题二给出的微分方程组，可以得出：(4.4.2)而方程组中添加的项分别为与.由于 A 种鱼为 B 种鱼提供食物，此处表示的则是单位数量A 种鱼（相对于）提供的供养 B 种鱼的食物量为单位数量的B 种鱼（相对于）消耗的供养 B 的食物量的倍。新的微分方程中出现了新的参数、和、，这两组参数拟合并不困难，且以问题三中的拟合程度来看，和对于最终结果的影响有限，究其结果，可以理解为B种鱼已构成了造成A 种鱼无法按 J 型曲线繁殖的大部分环境因素，此处为除了 B种鱼以外其余影响A 种鱼生存的环境制约因素

30、， 因此制约效果有限， 同理也是如此。对于和可以直接通过图像求得，而与则因为a,b, c,d可求 , 故根据（ 4.4.2）式可求。4.5 问题五作为种群数量阻滞增长的离散模型， 虽然种群数量一直在随时间变化， 但是由前几问的结论可以判断出较为明显的周期性。通过上述几个问题的讨论与求解，我们发现，参数组决定了 A 种鱼与 B 种鱼数量的走向趋势，而作为这两个种群的基本特征，不会随着时间的变化而变化，例如 a 是 A 种鱼的固有增长率， 而 c 则是 B 种鱼的死亡率， 这两个参数是不可能随时间产生较为显著的改变的。 因此这种趋势无关于初值 () 的选取，例如问题三中的原始数据图，如图 4.5.

31、1:图 4.5.1假设序列已知，序列未知，的曲线在图中的形状、周期都是确定的， t=11 时， B 种鱼依赖于 A 种鱼生存，因此，只需要知道任意一个 t 对应的数值，即可通过问题三中提到的在多个小区间内求积分的方法求得整个序列， 且序列参数一定。不妨假设一个 t=11 时的值，只要的取值在 y 的值域之内，则通过平移的图像，使得 t=11 时，图像与轴交点为，而不影响图像的形状、周期。五、模型的评价及推广优点：对于问题一的讨论， 就指出了种群数量增长的阻滞作用， 提出的单种群模型也较为新颖，并较为理想地刻画了 A 种鱼的数量随时间的变化情况， 这为此后问题四的解决提前做好了铺垫； 第三题对于

32、模型的参数进行拟合的过程中， 多种软件的辅助使得优势互补，很大程度上提高了拟合结果的精确程度。缺点：问题一、三的拟合过程都不可避免地存在误差。推广：典型的具有捕食关系的两个鱼类种群的数量变化模型可以被推广到生物领域的构成食饵捕食者关系的生态系统中， 用于观察并检测各类物种数量的阻滞增长模型，从而在更多领域得到应用。六、参考文献1 姜启源，谢金星，叶俊，数学模型（第四版），北京，高等教育出版社， 2011.1 ；2 龙辉平 , 习胜丰 , 侯新华 . 实验数据的最小二乘拟合算法与分析 J. 计算技术与自动化 . 2008(03)七、附录软件名称： matlab R2011b1stOpt 1.0V

33、C程序源码：1. 任务问题一 1stOpt 画图Parameters r,n; Variable t,x ;Function x=25.81*(1+r*(1-x/n)t ; Data;/t:x:0 25.811 35.472 51.083 60.454 65.875 92.826 120.647 149.208 184.549 222.6210 257.6811 292.5312 325.2613 352.3814 385.6315 406.9216 427.3017 433.8718 439.4219 451.381. 4.1.3.mt=0:1:19;x=25.81,35.47,51.08,

34、60.45,65.87,92.82,120.64,149.20,184.54,222.62,25 7.68,292.53,325.26,352.38,385.63,406.92,427.30,433.87,439.42,451.38; f=(r,t) 25.81*(1+r(1)*(1-t(:,2)/r(2).t(:,1)-t(:,2);tt=t'xx=x'x=tt,xx;r0=0.5,500;warning offr=nlinfit(x,zeros(size(x,1),1),f,r0);disp('拟合得到的 r 为')disp(num2str(r);plot(

35、x(:,1),x(:,2),'*');hold on;g=(x,y) 25.81*(1+r(1)*(1-y/r(2)x-y;ezplot(g);2. f2.mfunction f=f2(t,x)a=2;b=0.01;c=1;d=0.01;f=a*x(1)-b*x(1)*x(2);-c*x(2)+d*x(1)*x(2);end2.f3.mfunction f=f3(t,x)a=0.2;b=0.001;c=0.1;d=0.001;f=a*x(1)-b*x(1)*x(2);-c*x(2)+d*x(1)*x(2);end3.fun.mx0=10,100;t,x=ode45('f

36、2',0:10:100,x0);plot(t,x);legend('x','y');xlabel('t');figure;plot(x(:,1),x(:,2);xlabel('x');ylabel('y');3.fun3.mz0=45.79,41.40;t,x=ode45('f3',11:1:100,z0);plot(t,x);legend('x','y');xlabel('t');figure;plot(x(:,1),x(:,2);xlabe

37、l('x');ylabel('y');3.count.cpp#include <iostream>#include <fstream>using namespace std;#define a 2.2#define b 0.012#define c 1#define d 0.01double fx(double x,double y)double p;p=a-b*y;p=x-x*p*0.00001;return p;double fy(double x,double y)double q;q=-c+d*x;q=y-y*q*0.00001;

38、return q;int main()double x090;double y090;x00=41.59;y00=41.40;double xt,yt;double m=x00;double n=y00;for(int g=0;g<11;g+)for(int h=0;h<10000;h+)xt=fx(m,n);yt=fy(m,n);m=xt;n=yt;x0g+1=xt;y0g+1=yt;ofstream f1("d:4.txt");for(int i=0;i<90;i+)f1<<"t="<<22-(i+11)&l

39、t;<",x="<<x0i<<",y="<<y0i<<endl;for(int i=0;i<90;i+)cout<<"t="<<22-(i+11)<<",x="<<x0i<<",y="<<y0i<<endl;system("pause");return 0;3.main.cpp#include <iostream>usin

40、g namespace std;double fx(double x,double y,double a,double b);double fy(double x,double y,double c,double d);double s(double x,double y);int main()doublex90=45.79,53.03,64.05,75.40,90.36,107.14,127.79,150.77,179.65,211.82,249.91,291.31,334.95,380.67,420.28,445.56,447.63,414.04,347.04,265.33,187.57,

41、128.00,85.25,57.17,39.96,29.22,22.30,16.52,14.41,11.58,10.41,10.17,7.86,9.23,8.22,8.76,7.90,8.38,9.53,9.33,9.72,10.55,13.05,13.58,16.31,17.75,20.11,23.98,28.51,31.61,37.13,45.06,53.40,62.39,72.89,86.92,103.32,121.70,144.86,171.92,202.51,237.69,276.77,319.76,362.05,400.11,427.79,434.56,410.31,354.18,

42、278.49,203.72,141.06,95.08,66.76,45.41,33.13,25.89,20.51,17.11,12.69,11.76,11.22,10.29,8.82,9.51,8.69,9.53,8.68,10.82; doubley90=41.40,38.90,36.78,36.04,33.78,35.40,34.68,36.61,37.71,41.98,45.72,53.10,65.44,83.00,108.74,150.01,205.61,281.60,364.56,440.30,489.68,512.95,510.01,491.06,462.22,430.15,396

43、.95,364.87,333.16,304.97, 277.73,253.16,229.66,209.53,190.07,173.58,156.40,143.05,130.75,117.49,108.16,98.08,88.91,82.28,75.42,69.58,62.58,59.22,54.91,49.79,45.94,43.41,41.30,40.28,37.71,36.58,36.98,36.65,37.87,39.63,42.97,46.95,54.93,64.61,81.28,105.50,143.03,192.45,260.84,339.39,413.79,466.94,494.

44、72,499.37,484.58,460.63,429.79,398.77,366.49,336.56,306.39,279.53,254.95,233.50,212.74,193.61,175.01,160.59,146.12,131.85 ;double x090;double y090;double min=1000000;doublea20=1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5, 2.6,2.7,2.8,2.9;doubleb20=0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.011,0.012,0.013,0.014,0.015,0.016,0.017,0.018,0.019,0.02;doublec20=0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1;doubled20=0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.011,0.012,0.013,0.014,0.015,0.016,0.01