地中海鲨鱼问题的定量描述.doc

题目: 地中海鲨鱼问题的定量描述内容摘要：以下通过数学方法对地中海鲨鱼问题进行了定量分析。通过建立了相应的微分方程模型，然后求其相应的数值解，以图像的方式，清晰且直观的反映了战争前和战争中鲨鱼种群数量，鲨鱼百分数以及食饵种群数量随时间的变化。然后根据其变化特点呈周期性，基于此基础上，我们将傅里叶级数应用于此。利用傅里叶级数将数值解的图像进行数据拟合做出一定的结果分析。建立了一个相对简单的n次三角多项式来近似代替了复杂的函数等。通过利用图像，公式的结合，既清晰又定量地描绘了这一变化。最终我们得出结论：（1）一般情况下食饵种群数量总是多于鲨鱼种群数量；战争中的鲨鱼百分数多于战前时的。（2）人工捕获量的降低，导致鲨鱼在一段时间内大幅度的增加，乃至超过食饵的数量。（3）鲨鱼种群数量，所占百分数以及食饵种群数量随时间呈周期性变化。关键词：微分方程模型 数值解 n次三角多项式 傅里叶级数 数据拟合 问题重述：意大利生物学家Ancona安科纳曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼等的比例有明显增加（见下表），而供其捕食的食用鱼的百分比却明显下降.显然战争使捕鱼量下降，食用鱼增加，鲨鱼等也随之增加，但为何鲨鱼的比例大幅增加呢？ 年代19141915191619171918百分比11.921.422.121.236.4年代19191920192119221923百分比27.316.015.914.819.7模型假设与符号说明：假设（一）：食饵由于捕食者的存在使增长率降低（种群数量减少），假设降低的程度与捕食者数量成正比；假设（二）：捕食者由于食饵为它提供食物的作用使其死亡率降低或使出生率增加，假定增长的程度与食饵数量成正比。假设（三）：由于资源，空间等有限。因此种群之间存在着相互竞争，相互制约等，使之不能够无限增长下去，它应存在某一最大值，此时增长率为零。假定增长率为种群数量的减函数且成正比例关系。假设（四）：每种种群的内部制约(阻滞作用)是相互独立的，比如食饵种群数量不会与鲨鱼进行竞争。反之亦然。假设（五）：忽略环境等因素的变化对种群数量的影响。假设（六）：以上考虑的各种因素同时存在时对种群数量的总影响可以叠加。相关符号变量的说明：食饵在t时刻的数量； 捕食者在t时刻的数量；（）=食饵独立生存时的增长率；()捕食者独自存在时的死亡率；一一食饵的固有增长率； 一一捕食者的固有增长率；捕食者掠取食饵的能力； 食饵对捕食者的供养能力.一一一定条件下该生态系统所能容纳食饵的最大数量（即环境容纳量）； 一一一定条件下该生态系统所能容纳鲨鱼的最大数量；e捕获能力系数（可以指人类对其的捕获）；一一鲨鱼的百分数随时间变化的函数；模型的分析与建立：由假设（一）和假设（六）知：从时刻有：；由假设（二），假设（三）和假设（四）可知：（；）当时即： 所以有：；所以综合上述所有假设可得如下模型:；因此由上述模型可以针对一组具体的数据用Matlab软件进行计算.设： ；战前e=0.3;战争中e=0.1；将上述数据分别代入得如下两个模型：战前和战争中：；模型的求解：在matlab中输入如下命令得战前的数值解：(时间取0-15)；首先，建立m-文件shier.m如下： function dx=shark1(t,x) dx=zeros(2,1); dx(1)=x(1)*(0.7-0.1*x(2)-x(1)/1000);dx(2)=x(2)*(-0.8+0.02*x(1)+0.5*x(2)/100);其次，建立主程序shark.m如下： t,x=ode45(shark1,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2)求解结果如下：图表 1图表 2图表1反映了战前（即食饵和鲨鱼）随t（时间）的变化；其中蓝*线代表食饵，实线代表鲨鱼。同理输入如下命令（数据仍按上述计算）得战争中的结果：function dx=shark2(t,x) dx=zeros(2,1); dx(1)=x(1)*(0.9-0.1*x(2)-x(1)/1000);dx(2)=x(2)*(-0.6+0.02*x(1)+0.5*x(1)/100);其次，建立主程序shark.m如下： t,x=ode45(shark2,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2)结果如下：图表 3图表 4 输入如下命令得鲨鱼所占鱼总数的百分数随时间的变化情况：t,x=ode45(shark1,0 15,25 2)；plot(t,x(:,2)./(x(:,2)+x(:,1),r)hold ont,x=ode45(shark2,0 15,25 2);plot(t,x(:,2)./(x(:,2)+x(:,1),g)结果如下：图表 5将题目中的表格数据一起与上述图4画成折线图如下：hold ont=0:1:9;x=0.01*11.9 21.4 22.1 21.2 36.4 27.3 16 15.9 14.8 19.7;plot(t,x,g)结果如下：图表 6模型的检验与结果分析：由图表1， 3， 5， 6可以猜测到鲨鱼百分数或者是鲨鱼和食饵种群数量无论是战争中还是战争前均为时间的周期函数，并且从图表6中还可以看出题目当中的实际鲨鱼百分含量随时间变化也有周期性（图表6中的黄线），且变化走势与模型中所设定的相一致。因此可以初步断定只要我们选取适当的参数(比如；等)就可以得到与题目中的实际数据相吻合的曲线图，故该模型可以应用于此问题当中。 上述的图表1，图表3，分别反映了战争中跟战争前鲨鱼种群数量跟食饵种群数量随时间的变化规律。图表5，图表6.，分别反映了战争中跟战争前鲨鱼百分数随时间的变化规律。由于这种变化呈周期性的变化（作为数学上的周期函数），基于此种情况之下，我们可以对其作进一步的数据拟合，建立一个相对简单的多项式函数来计算任意时刻的鲨鱼百分数为多少，下面进行数据拟合构造一个相对简单的函数式来刻画这种变化（即用来计算任意时刻的鲨鱼百分数）。（仍以模型中所设定的数据为例）。由图表5我们可以看出鲨鱼百分数跟时间的函数的图像为一连续曲线，且为周期函数。由高等数学知识我们知道，若一个函数是以为周期的函数，且在上按段光滑，则函数在上收敛于它的傅里叶级数，即：由图像可直接看出周期应该为T=9.647；即： 9.647；因此在一定的截取误差允许的范围内我们可以利用上述级数的其中一部分来近似计算任意时刻的鲨鱼百分数，即用一个n次三角多项式来近似代替鲨鱼百分数随时间变化的函数，以n=10为例。因此由上述可知:令则由线性最小二乘法原理得如下超定方程组：其中A=(RTR)-1RTW.我们只要求出A即可。若以以上战前的鲨鱼百分数数据为例（n值取10）则得如下结果（求解的matlab命令见以下的附录）：A=0.3612 -0.1245 -0.1450 0.0071 0.0646 0.0074 -0.0059 0.0020 0.0031 0.0013 -0.0006 0.0006 -0.0001 0.0001 -0.0003 0.0003 -0.0007 0.0001 -0.0008 -0.0003 -0.0002;令B=则有 .*；由*式即可计算战前任意时刻的鲨鱼百分数为多少。为了验证其计算结果是否可靠，将其计算结果与上述图表6中的红色曲线（即战前的鲨鱼百分数）一并放入同一个坐标当中，并且分别取做数据拟合，结果分别如下（相应的计算程序见以下的附录中）：图表 7A=0.3618 -0.1233 -0.1456 0.0071 0.0646;B=;图表 8A=0.3615 -0.1236 -0.1456 0.0069 0.0646 0.0077 -0.0060; B=;图表 9 A=0.3613 -0.1238 -0.1456 0.0067 0.0646 0.0075 -0.0060 0.0021 0.0030 0.0013 -0.0007 0.0006 -0.0002;图表 10A=0.3612 -0.1245 -0.1450 0.0071 0.0646 0.0074 -0.0059 0.0020 0.0031 0.0013 -0.0006 0.0006 -0.0001 0.0001 -0.0003 0.0003 -0.0007 0.0001 -0.0008 -0.0003 -0.0002;B=;图表6至图表10中的红色曲线为以上模型中所求的数值解（战前的数据为例，同样的方法可以用来分析战争中的相关数据），不同的黑线分别代表了不同的n值所得到的结果，由上述的拟合结果可以看出通过增大n的值可以提高计算结果的精确度，并且n越大，计算结果就越精确。当时，计算结果就几乎可以和已知数据重合（图表10）。利用同样的方法可以分别找出鲨鱼以及食饵种群数量随时间变化的函数的近似表达式。综合以上的所有数据以及结果，可以得出如下结论:在上述所规定的所有假设条件下，（1）在一般情况之下食饵的种群数量总是位于鲨鱼种群数量上方，由上图表1中的结果可以看出。但由于战争的发生，改变了其一般情况之下的环境条件，即战前时和战争中的人工捕获量不同，从而导致鲨鱼的种群数量大幅度增加，乃至存在一段时间超过食饵的数量，如图表3中显示的结果。（2）战争中的鲨鱼百分数多于战前的（如图表6）。（3）食饵，鲨鱼的种群数量以及百分数随时间呈现周期性的变化。模型的优点及推广：在该模型中充分考虑了人工捕获量，鲨鱼对食饵的捕获能力以及食饵对鲨鱼的供养能力。没有把增长率作为一个常数来处理，而是将其看做一个随种群数量变化的减函数（），在此基础上，引入了固有增长率（）和增长率函数以及环境容纳量（）的相关概念。此模型不仅可以应用于地中海鲨鱼问题之中，它可以应用于其它一般情况之下很多存在简单捕食关系的生态系统中（在此不再列举实际的案例来进行讨论）。并且可以准确反映出该生态系统中捕食者和被吃者之间的相互作用关系，以及它们各自的种群数量随时间变化的情况，计算任意时刻捕食者的种群数量或被吃者的种群数量以及它们各自所占的百分数。也可对其相应的种群数量进行相应的未来变化的预测。 在结果分析中所使用的数学方法，可以应用于一切呈周期性变化的数据拟合中。它比起使用一般的n次多项式拟合要优越的多。模型的缺点:（1）该模型中没有考虑到环境，气温，等因素的变化以及种群年龄结构的分布对其种群数量的影响。因此只能适用于一个小范围的生态系统在短时间内的变化情况。（2）仅仅考虑了鲨鱼和食饵这两个种群，并未考虑其它种群的存在而对其产生的影响。并且还假设了各种群的阻滞作用是相互独立的，但在实际的两个存在捕食关系的种群之间，它们不仅仅存在捕食关系，它们也存在一定的竞争。所以对于十分复杂的生态系统来讲，该模型并不成立。（3）对于影响鲨鱼的增长率的因素中，也不仅仅存在人类的捕获。并且更一般情况下的捕获量并非为一常数，它也应该是一个变量，就像上述的增长率一样。参考文献：数学实验与MATLAB第五版；华东大学出版社数学分析下册第三版；华东师范大学数学系主编；高等教育出版社经济数学一一线性代数第二版；吴传生主编；高等教育出版社矩阵理论及其应用数学建模 附录：n=2时的数据拟合MATLAB程序：o=1/2*ones(49,1);t1=0:9.647/48:9.647;t2=(t1);n1=cos(1*pi*t2)/4.8235);m1=sin(1*pi*t2)/4.8235); n2=cos(2*pi*t2)/4.8235);m2=sin(2*pi*t2)/4.8235); R=o n1 m1 n2 m2 ;t,x=ode45(shark4,0 9.647,25 2);w=x(:,2)./(x(:,2)+x(:,1);A=Rw;A1=A;syms tn1=cos(1*pi*t)/4.8235);m1=sin(1*pi*t)/4.8235); n2=cos(2*pi*t)/4.8235);m2=sin(2*pi*t)/4.8235); B=0.5 n1 m1 n2 m2 ; W=dot(A1,B) .,; W =1*1 .,; 6517300416699079/36028797018963968-8885740318467653/72057594037927936*cos(2000/9647*pi*t)-5244734929237683/36028797018963968*sin(2000/9647*pi*t)+4121111377597937/576460752303423488*cos(4000/9647*pi*t)+4652180975160443/72057594037927936*sin(4000/9647*pi*t) f=(t)6517300416699079/36028797018963968-8885740318467653/72057594037927936*cos(2000/9647*pi*t)-5244734929237683/36028797018963968*sin(2000/9647*pi*t)+4121111377597937/576460752303423488*cos(4000/9647*pi*t)+4652180975160443/72057594037927936*sin(4000/9647*pi*t); t=0:9.647/48:9.647; w=f(t); t,x=ode45(shark4,0 9.647,25 2); plot(t,w,ko,t,x(:,2)./(x(:,2)+x(:,1),r)n=3时的数据拟合MATLAB程序：o=1/2*ones(49,1);t1=0:9.647/48:9.647;t2=(t1);n1=cos(1*pi*t2)/4.8235);m1=sin(1*pi*t2)/4.8235); n2=cos(2*pi*t2)/4.8235);m2=sin(2*pi*t2)/4.8235); n3=cos(3*pi*t2)/4.8235);m3=sin(3*pi*t2)/4.8235); R=o n1 m1 n2 m2 n3 m3 ;t,x=ode45(shark4,0 9.647,25 2);w=x(:,2)./(x(:,2)+x(:,1);A=Rw;A1=A;syms tn1=cos(1*pi*t)/4.8235);m1=sin(1*pi*t)/4.8235); n2=cos(2*pi*t)/4.8235);m2=sin(2*pi*t)/4.8235); n3=cos(3*pi*t)/4.8235);m3=sin(3*pi*t)/4.8235); B=0.5 n1 m1 n2 m2 n3 m3 ; W=dot(A1,B) .,; W =1*1 ., 101751143016943/562949953421312-8906649372926567/72057594037927936*cos(2000/9647*pi*t)-5244734929237683/36028797018963968*sin(2000/9647*pi*t)+123557466935207/18014398509481984*cos(4000/9647*pi*t)+4652180975160443/72057594037927936*sin(4000/9647*pi*t)+4432719545289157/576460752303423488*cos(6000/9647*pi*t)-6971811585880695/1152921504606846976*sin(6000/9647*pi*t) f=(t)101751143016943/562949953421312-8906649372926567/72057594037927936*cos(2000/9647*pi*t)-5244734929237683/36028797018963968*sin(2000/9647*pi*t)+123557466935207/18014398509481984*cos(4000/9647*pi*t)+4652180975160443/72057594037927936*sin(4000/9647*pi*t)+4432719545289157/576460752303423488*cos(6000/9647*pi*t)-6971811585880695/1152921504606846976*sin(6000/9647*pi*t); t=0:9.647/48:9.647;w=f(t);t,x=ode45(shark4,0 9.647,25 2);plot(t,w,k*,t,x(:,2)./(x(:,2)+x(:,1),r)n=6时的数据拟合MATLAB程序：o=1/2*ones(49,1);t1=0:9.647/48:9.647;t2=(t1);n1=cos(1*pi*t2)/4.8235);m1=sin(1*pi*t2)/4.8235); n2=cos(2*pi*t2)/4.8235);m2=sin(2*pi*t2)/4.8235); n3=cos(3*pi*t2)/4.8235);m3=sin(3*pi*t2)/4.8235); n4=cos(4*pi*t2)/4.8235);m4=sin(4*pi*t2)/4.8235); n5=cos(5*pi*t2)/4.8235);m5=sin(5*pi*t2)/4.8235); n6=cos(6*pi*t2)/4.8235);m6=sin(6*pi*t2)/4.8235); R=o n1 m1 n2 m2 n3 m3 n4 m4 n5 m5 n6 m6;t,x=ode45(shark4,0 9.647,25 2);w=x(:,2)./(x(:,2)+x(:,1);A=Rw;A1=A;syms tn1=cos(1*pi*t)/4.8235);m1=sin(1*pi*t)/4.8235); n2=cos(2*pi*t)/4.8235);m2=sin(2*pi*t)/4.8235); n3=cos(3*pi*t)/4.8235);m3=sin(3*pi*t)/4.8235); n4=cos(4*pi*t)/4.8235);m4=sin(4*pi*t)/4.8235); n5=cos(5*pi*t)/4.8235);m5=sin(5*pi*t)/4.8235); n6=cos(6*pi*t)/4.8235);m6=sin(6*pi*t)/4.8235); B=0.5 n1 m1 n2 m2 n3 m3 n4 m4 n5 m5 n6 m6; W=dot(A1,B) W =1*1 ., 1627360572103877/9007199254740992-278661650987561/2251799813685248*cos(2000/9647*pi*t)-5244734929237685/36028797018963968*sin(2000/9647*pi*t)+7739302545047061/1152921504606846976*cos(4000/9647*pi*t)+1163045243790111/18014398509481984*sin(4000/9647*pi*t)+8697063751772127/1152921504606846976*cos(6000/9647*pi*t)-3485905792940369/576460752303423488*sin(6000/9647*pi*t)+4788843151125117/2305843009213693952*cos(8000/9647*pi*t)+7027773051795359/2305843009213693952*sin(8000/9647*pi*t)+6080212008601511/4611686018427387904*cos(10000/9647*pi*t)-6349771707463659/9223372036854775808*sin(10000/9647*pi*t)+1431694478914249/2305843009213693952*cos(12000/9647*pi*t)-486902548663185/2305843009213693952*sin(12000/9647*pi*t) f=(t)1627360572103877/9007199254740992-278661650987561/2251799813685248*cos(2000/9647*pi*t)-5244734929237685/36028797018963968*sin(2000/9647*pi*t)+7739302545047061/1152921504606846976*cos(4000/9647*pi*t)+1163045243790111/18014398509481984*sin(4000/9647*pi*t)+8697063751772127/1152921504606846976*cos(6000/9647*pi*t)-3485905792940369/576460752303423488*sin(6000/9647*pi*t)+4788843151125117/2305843009213693952*cos(8000/9647*pi*t)+7027773051795359/2305843009213693952*sin(8000/9647*pi*t)+6080212008601511/4611686018427387904*cos(10000/9647*pi*t)-6349771707463659/9223372036854775808*sin(10000/9647*pi*t)+1431694478914249/2305843009213693952*cos(12000/9647*pi*t)-486902548663185/2305843009213693952*sin(12000/9647*pi*t);t=0:9.647/48:9.647;w=f(t);t,x=ode45(shark4,0 9.647,25 2);plot(t,w,ks,t,x(:,2)./(x(:,2)+x(:,1),r)n=10时的数据拟合MATLAB程序：o=1/2*ones(49,1);t1=0:9.647/48:9.647;t2=(t1);n1=cos(1*pi*t2)/4.8235);m1=sin(1*pi*t2)/4.8235); n2=cos(2*pi*t2)/4.8235);m2=sin(2*pi*t2)/4.8235); n3=cos(3*pi*t2)/4.8235);m3=sin(3*pi*t2)/4.8235); n4=cos(4*pi*t2)/4.8235);m4=sin(4*pi*t2)/4.8235); n5=cos(5*pi*t2)/4.8235);m5=sin(5*pi*t2)/4.8235); n6=cos(6*pi*t2)/4.8235);m6=sin(6*pi*t2)/4.8235); n7=cos(7*pi*t2)/4.8235);m7=sin(7*pi*t2)/4.8235); n8=cos(8*pi*t2)/4.8235);m8=sin(8*pi*t2)/4.8235); n9=cos(9*pi*t2)/4.8235);m9=sin(9*pi*t2)/4.8235); n10=cos(10*pi*t2)/4.8235);m10 =sin(10*pi*t2)/4.8235);R=o n1 m1 n2 m2 n3 m3 n4 m4 n5 m5 n6 m6 n7 m7 n8 m8 n9 m9 n10 m10;t,x=ode45(shark4,0 9.647,25 2);w=x(:,2)./(x(:,2)+x(:,1);A=Rw;A1=A;syms tn1=cos(1*pi*t)/4.8235);m1=sin(1*pi*t)/4.8235); n2=cos(2*pi*t)/4.8235);m2=sin(2*pi*t)/4.8235); n3=cos(3*pi*t)/4.8235);m3=sin(3*pi*t)/4.8235); n4=cos(4*pi*t)/4.8235);m4=sin(4*pi*t)/4.8235); n5=cos(5*pi*t)/4.8235);m5=sin(5*pi*t)/4.8235); n6=cos(6*pi*t)/4.8235);m6=sin(6*pi*t)/4.8235); n7=cos(7*pi*t)/4.8235);m7=sin(7*pi*t)/4.8235); n8=cos(8*pi*t)/4.8235);m8=sin(8*pi*t)/4.8235); n9=cos(9*pi*t)/4.8235);m9=sin(9*pi*t)/4.8235); n10=cos(10*pi*t)/4.8235);m10 =sin(10*pi*t)/4.8235); B=0.5 n1 m1 n2 m2 n3 m3 n4 m4 n5 m5 n6 m6 n7 m7 n8 m8 n9 m9 n10 m10; W=dot(A1,B) .,; W =1*1 .,; 6509329302005301/36028797018963968-8917624777242777/72057594037927936*cos(2000/9647*pi*t)-5244734929237683/36028797018963968*sin(2000/9647*pi*t)+3866035707396947/576460752303423488*cos(4000/9647*pi*t)+4652180975160441/72057594037927936*sin(4000/9647*pi*t)+8689832621518959/1152921504606846976*cos(6000/9647*pi*t)-3485905792940353/576460752303423488*sin(6000/9647*pi*t)+4774380890618715/2305843009213693952*cos(8000/9647*pi*t)+7027773051795459/2305843009213693952*sin(8000/9647*pi*t)+6051287487588919/4611686018427387904*cos(10000/9647*pi*t)-3174885853731663/4611686018427387904*sin(10000/9647*pi*t)+5668928873631551/9223372036854775808*cos(12000/9647*pi*t)-7790440778612307/36893488147419103232*sin(12000/9647*pi*t)+2319307828488339/18446744073709551616*cos(14000/9647*pi*t)-1692162567446817/4611686018427387904*sin(14000/9647*pi*t)+5529304472617079/18446744073709551616*cos(16000/9647*pi*t)-3671878516630713/4611686018427387904*sin(16000/9647*pi*t)+3808957445717893/36893488147419103232*cos(18000/9647*pi*t)-3853959857305949/4611686018427387904*sin(18000/9647*pi*t)-6224299460416903/18446744073709551616*cos(20000/9647*pi*t)-1337409325405315/461168601842738790