鱼类物种的数量模型.doc

鱼类物种的数量模型摘要本文针对种鱼以及其捕食者种鱼各自随时间的数量变化，估计出种鱼和种鱼之间的相互影响，在此基础上，依据问题中给出的微分方程组，进行数据拟合以及修改，从而最终得出两物种的数量模型。针对问题一，根据已给数据，运用绘制图形，发现种鱼的数量符合阻滞增长模型（模型），由此建立与的数学关系为，于是运用拟合，计算出种鱼的自然增长率以及环境资源容许的种群最大数量，总体拟合程度较高。针对问题二，将数值以及初始值带入，运用计算出微分方程的数值解和并作图。针对问题三，根据已知数据，首先分别绘制及的图像，再将问题二中数值带入，按照问题二的方法计算出两组数据, ,并与原数据比较，发现图形走势基本一致，再根据数据初步拟合的结果显示，拟合后的图像与原图像横坐标相差倍，故通过遍历()邻域内的值优化求解。然后通过控制变量，采用最小二乘拟合法循环逐个拟合四个参数，用平方和残差值的最小值来评价参数的优劣程度，从而筛选出最优解，并代入微分方程组，将再次计算出的数据与原数据比较，发现拟合程度比较理想。拟合过程中还采用微元的思想，将整个时间区间分成无限个小区间进行积分求和运算，得到的数值解为 ，最后根据求出两个初值分别为。针对问题四，通过观察问题三绘制出的拟合图像，发现计算出的数据比原数据略大，说明原微分方程组是理想模型，没有考虑起物种自身阻滞作用的项，于是在原有的微分方程基础上增加了项，提高了数据的拟合程度。针对问题五，分析了微分方程组中的各参数代表的含义，得出不依赖于数据的结论，即无论y取何种初值，图像的形状不变，只是在原有基础上进行平移，基于这个结论， 的数值可解，但对于y的初值无法给定。本文亮点首先在于问题三中的数据拟合过程，通过对拟定的未优化数据与原数据各自形成的图像特征的比较，确定准确值存在的大致范围，运用最小二乘拟合法逐个拟合优化，最终取得最优解。然后在问题五中，清楚地分析了与、的关系中最重要的一点就是和各时刻对应的数值无关，而是和图像的形状、周期有关。关键词：阻滞作用、最小二乘拟合、一、 问题重述在某水域生长的一种鱼类（A种鱼），该水域出现了一种以A种鱼为食的另一种鱼类（B种鱼），记时刻A种鱼的数量为，B种鱼的数量为。A种鱼与B种鱼的数量相互影响，本题要求通过建立数学模型解决以下五个问题：1 根据表格中记录的一段时间内所观测到的A种鱼在若干观测观测的数量，描述A种鱼的数量随时间变化的情况。2 如果这两种鱼的数量变化符合以下的微分方程组：根据题中给出的及初始值，计算t=10,20,100时两种鱼的数量。3 根据在一段时间内观测到的两种鱼的数量，如果两种鱼的数量变化符合问题二中的微分方程组，求其中的常数以及初始时刻两种鱼的数量。4 问题二中的微分方程组能否很好的反映问题三数据？请你修改微分方程组的形式，使得改进后的微分方程组更好的反映问题三中的数据。5 如果在问题三中，无法观测到B种鱼的数据（），只能观测到A种鱼的数据（），是否可以求解与问题三与问题四相应的问题。二、 模型假设1. A种鱼的数量变化符合阻滞增长模型；2. 除了食饵、捕食者关系以及自身阻滞作用，其他因素对两种鱼的数量无影响；三、 符号说明N：环境资源所能容纳的最大种群数量；r：种群数量的自然增长率；：A种鱼在初始时刻的种群数量；：B种鱼在初始时刻的种群数量；：平方和残差值的最小值；：表示的是单位数量A种鱼提供的供养B种鱼的食物量为单位数量的B种鱼消耗的供养B的食物量的倍； ：表示单位数量B种鱼消耗A种鱼作为食物量为单位数量的A种鱼自身减少数量的倍。四、 模型的建立与求解4.1问题一4.1.1问题一的分析当水域中只有A种鱼时便构成了单种群模型，为了清楚的认识问题一中A种鱼的数量变化规律，根据A种鱼在从0到19之间不同时刻的数量，首先运用绘制出A种鱼数量随时间变化的图像如下（如图4.1.1-1），它大体上描述了该种群数量随时间变化的规律，从而判断出A种鱼符合的模型，并根据数据拟合出模型中的未知参数。 图4.1.1-14.1.2模型的建立根据图4.1.1-1，不难看出A种鱼的增长规律符合阻滞增长模型（logistic模型）。将t时刻的种群数量记为，并将视为连续、可微函数，记初始时刻的数量为,假设A种鱼的自然增长率为常数r，可以注意到r随着数量的增加而下降，即表现为阻滞越来越大。若将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模型的结果分析由图像的拟合程度可以直观地判断出，当水域中只有A种鱼时，该种群的数量变化满足模型，即自身增长阻滞因素的存在是肯定的。4.2问题二4.2.1问题二的分析问题二中出现了以A种鱼为食的B种鱼，要求通过问题中给出的微分方程组，相关系数的值以及两种鱼数量的初始值计算出两种鱼在不同时刻的数量，这属于中微分方程的初值解问题，故直接根据其中数据编写程序求解即可。4.2.2模型的求解已知微分方程组的形式为，根据问题二中的数据，a=2，b=0.01，c=1，d=0.01，另在初始时刻t=0时，有。求解如下(精确到小数点后13位)：106.69561062293209242.8269687425382043.4135062005784555.56049274398930402.368367290486124.6846053281384096.307582411949135.17813481510405019.302277240457459.7529087274689606.50808159252891150.1411739936237012.2164610147701396.89782586901080322.124619684216432.70427018826290182.57493739181940.505424275731910035.104638631399344.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倍，因此，记问题二中参数为(),考虑将问题二中缩小10倍后得()，以时的数值为初值，利用微元法，将分成无数个极小的区间，并在每个小区间内对微分方程进行积分处理，通过编程还原拟合后的数值，即时的数值，所得数据结果如下(此处只列出前20组数据，精确到小数点后13位)：t x(t) y(t) 1 45.7900000000000 41.40000000000002 53.7163569677204 39.36669859669533 63.1286142321031 37.76109952492154 74.2949284523748 36.59263557351635 87.5208046082887 35.89010074330116 103.145459762021 35.71215140919517 121.541823058173 36.14730632021508 143.061062641112 37.35879432925819 168.167999347508 39.495438265894510 197.213300867108 42.843281425403811 230.256935482831 47.950682837306512 267.068172070766 55.628317265363613 306.905500686171 66.978236135377514 347.801978090903 84.121052912449215 385.749575134144 110.17318125515616 414.364741995805 148.96946198322017 424.887801816950 205.06275076590118 407.479644798270 282.04770131633919 358.590063197010 375.25628833583720 286.887513576624 469.433817327986另外，利用问题二中的参数得到拟合曲线，如图4.3.1-2：图4.3.1-2作出关于比较图像，如图4.3.1-3：图4.3.1-3作出关于的比较图像，如图4.3.1-4：图4.3.1-4观察图像,可以判断以()为参数拟合的方程所得到的数据与问题三的数据有着较高的拟合度，因此可以在()各参数的相应邻域内考察更优解。4.3.2模型的求解根据以上分析，可以通过遍历()邻域内的值求得最优解，由于分析中认为问题三数据较好地坐落在函数曲线上，因此，任意数据作为初值计算而得的曲线具有相同的参数与相似的表达式，仅相位有差异，不妨选取处的值,计算相应参数值对应的时的值，并规定用平方和残差值的最小值来评价参数的优劣程度，其中n=90为拟合数据长度。利用以下式子：以及MATLAB软件进行最小二乘拟合2。若对于各遍历n个数据，则时间复杂度为，超出了可接受范围，考虑到已经拥有拟合度比较高的()为，不妨在考虑1个参数时令其余个不变，得到较优参数后立即替代原有参数，并逐个多次进行拟合修正，令拟合结果不断逼近真实最优结果，则时间复杂度将为。利用上述算法，求得较优参数，并得到优化参数拟合序列中的90组数据如下（精确到小数点后13位）：tt1145.790000000000041.40000000000005615.138861457279969.00309568748901254.360795919177839.37864667504955717.422343520928163.45918163106331364.672210450977637.81424460432735820.169980662240758.50728347922591477.066277838805836.72247234591745923.482194752630154.10524540527241591.933719213754536.14393592777826027.481263677354550.217510783100216109.70421118358836.15935520908276132.311321335517846.815122187207917130.84638583274336.88956388021856238.138357640069943.875721392695618155.77628917853738.55466514941876345.136713722011241.402188592139819185.06346638619441.36381849196196453.531307510573039.402859694033220219.01922459118845.73290922348496563.661063058341337.831650334543821257.39915055788752.45299846131746675.890794672368136.687433354563922299.41248660768962.67737517886776790.607178035897636.018749193867223343.45032550993277.868865866179368108.21875020836735.923805891108124384.385849491488101.69953278647569129.15590962540436.550479083822025414.728642761561138.09323862032570153.87091609883238.096312008425426425.268162794879190.41725398555471182.83789081666440.808515500215627405.249485644033261.34831600190172216.33357034984144.929902318401728351.024016487800346.05565198239273254.66464046030551.204413980456329274.641783901213428.36968861623874297.08027467818861.165672455961830196.751522338243490.10805128961375340.79986744442876.907637653628231132.941595890738522.40912493943776381.009198336474101.0833417589683288.0038873030939527.20529120827677410.860432068282136.9048892342973358.6034216805480512.85359084402478421.472118490313188.1434568187323439.8453802685919487.25240956560279402.615608861006257.7546085486483528.1860215083149455.95261725743480350.596735684413340.7529514802883620.6983113177825422.67457147855181276.336736263480421.9499857999543715.8307206710633389.44584580155882199.352233884202483.9504060728033812.6053131914126357.39485067349183133.700452851551518.7537819472613910.4113961448450327.1058997676898488.2828551636930524.357614303775408.90954296930337298.8399933572798558.9715920647978510.263630920899417.88108718175481272.6731142710608640.2082156612421484.976897532318427.18653888312947248.5816280925678728.5160333628630453.981513400348436.73431499277001226.4932146991138820.9848902487768420.979015245146446.47277107575796206.2949568138268916.0778258606198387.988162216631456.36771240761249187.8629300749419012.8230217278798356.140163089553466.39894786128350171.0714301236519110.6064174083346326.024374711171476.55515937956694155.796320241832929.08660959806996297.907333751208486.83299283564912141.915797341273938.04579547838872271.868764055616497.23594100965142129.313878280482947.34319024571603247.888308198598507.77381022211353117.885882250924956.88671058645735225.896141750932518.46222428117785107.534499072907966.62359151767271205.782313839817529.3224721691988198.1691962642994976.51954866447124187.4243467651255310.381508042742989.7062190405298986.55476218018557170.6979425876145411.671951232588582.0685903145883996.71798829256552155.4805915669265513.240968613729775.1877806758082007.00555086926924141.652416308371同时经过程序计算，t从0时刻到10时刻对应的x（t）与y（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从中也可以得到与的初始值分别为=8.68581，=100.0154.3.3模型的可信度检验根据上述得到的数据，可以得到以下通过优化的拟合参数图，如图4.3.2-1:图4.3.2-1优化拟合后的与原数据的比较图，如图4.3.2-2：图4.3.2-2优化拟合后的与原数据的比较图，如图4.3.2-3：图4.3.2-3同时两组数据的相轨线图像也显示了较为理想的拟合程度，如图4.3.2-4：图4.3.2-4借助这些图像可以直观地表达出，拟合后的参数使得计算结果与所给数据更加吻合。这也进一步证明了模型的可信度比较高，所得的参数值及变量初值与真实值有着较高的拟合度。同时由于相轨线是封闭曲线，可以确定的是 和均满足周期函数的特性。 4.4问题四观察问题三中的图像，发现拟合图像总是比真实图像高出一些，结合问题一中采用的人口环境阻滞增长模型，不难发现，任务问题二中的模型缺少了环境阻滞参数。因此提出以下的优化模型IV： （4.4.1）类比于问题二给出的微分方程组，可以得出： (4.4.2) 而方程组中添加的项分别为与.由于A种鱼为B种鱼提供食物，此处表示的则是单位数量A种鱼（相对于）提供的供养B种鱼的食物量为单位数量的B种鱼（相对于）消耗的供养B的食物量的倍。新的微分方程中出现了新的参数、和、，这两组参数拟合并不困难，且以问题三中的拟合程度来看，和对于最终结果的影响有限，究其结果，可以理解为B种鱼已构成了造成A种鱼无法按J型曲线繁殖的大部分环境因素，此处为除了B种鱼以外其余影响A种鱼生存的环境制约因素，因此制约效果有限，同理也是如此。对于和可以直接通过图像求得，而与则因为可求,故根据（4.4.2）式可求。4.5问题五作为种群数量阻滞增长的离散模型，虽然种群数量一直在随时间变化，但是由前几问的结论可以判断出较为明显的周期性。通过上述几个问题的讨论与求解，我们发现，参数组决定了A种鱼与B种鱼数量的走向趋势，而作为这两个种群的基本特征，不会随着时间的变化而变化，例如a是A种鱼的固有增长率，而c则是B种鱼的死亡率，这两个参数是不可能随时间产生较为显著的改变的。因此这种趋势无关于初值()的选取，例如问题三中的原始数据图，如图4.5.1:图4.5.1假设序列已知，序列未知，的曲线在图中的形状、周期都是确定的，t=11时，B种鱼依赖于A种鱼生存，因此，只需要知道任意一个t对应的数值，即可通过问题三中提到的在多个小区间内求积分的方法求得整个序列，且序列参数一定。不妨假设一个t=11时的值，只要的取值在y的值域之内，则通过平移的图像，使得t=11时，图像与轴交点为，而不影响图像的形状、周期。五、 模型的评价及推广优点：对于问题一的讨论，就指出了种群数量增长的阻滞作用，提出的单种群模型也较为新颖，并较为理想地刻画了A种鱼的数量随时间的变化情况，这为此后问题四的解决提前做好了铺垫；第三题对于模型的参数进行拟合的过程中，多种软件的辅助使得优势互补，很大程度上提高了拟合结果的精确程度。缺点：问题一、三的拟合过程都不可避免地存在误差。推广：典型的具有捕食关系的两个鱼类种群的数量变化模型可以被推广到生物领域的构成食饵捕食者关系的生态系统中，用于观察并检测各类物种数量的阻滞增长模型，从而在更多领域得到应用。六、 参考文献1姜启源，谢金星，叶俊，数学模型（第四版），北京，高等教育出版社，2011.1；2 龙辉平,习胜丰,侯新华. 实验数据的最小二乘拟合算法与分析J. 计算技术与自动化. 2008(03)七、附录软件名称：matlab R2011b 1stOpt 1.0 VC程序源码：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,60.45,65.87,92.82,120.64,149.20,184.54,222.62,257.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(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(f2,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);xlabel(x);ylabel(y);3.count.cpp#include #include 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;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;g11;g+)for(int h=0;h10000;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;i90;i+)f1t=22-(i+11),x=x0i,y=y0iendl;for(int i=0;i90;i+)coutt=22-(i+11),x=x0i,y=y0iendl;system(pause);return 0;3.main.cpp#include using 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()double x90=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,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,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;double y90=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.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.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;double a20=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;double b20=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;double c20=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;double d20=0.001,0.002,0.003,0.004,0