西北农林科技大学实验报告学院名称：理学院专业年级：2013级信计1班课程：数学模型与数学建模报告日期：2015年12月15日实验六机理模型与平衡原理实验目的如果对所研究的问题了解的比较深入，知道产生现象的内在的机理，那么依据机理建模，则模型具有更好的可靠性和广泛性。
不考虑随机因素，假设每一时刻是确定的如果对系统状态的观测和描述只在离散的时间点上，则构成差分方程模型；如果考虑系统随时间连续变化，则是微分方程模型。
本节主要以这两类方程为例，介绍用MATLAB^件求解机理模型的基本方法。
一、差分方程模型1.实验题目由一对兔子开始，一年可以繁殖出多少只兔子？如果一对兔子每个月可以生一对小兔子，兔子在出生两个月后就具有繁殖能力，由一对刚出生一个月的兔子开始，一年内兔子种群数量如何变化。
求这个种群的稳定分布和固有增长率。
2.实验内容解假设(a)兔子每经过一个月底就增加一个月龄；(b)月龄大于等于2的兔子都具有繁殖能力；(c)具有繁殖能力的兔子每一个月一定生一对兔子；(d)兔子不离开群体(不考虑死亡)记第n个月初的幼兔(一月龄兔)数量为ao(n),成兔(月龄大于等于2)数量为a1(n )，则兔子总数为a(n)= a o(n ) +a(n),平衡关系为：' 本月初幼兔数量 =上月初成兔数量〔本月初成兔数量 =上月初成兔数量+上月初幼兔数量建立模型：a o(n) =a“n T)g(n) =a°(n -1) +印(n -1)a。
⑴=g(1) =0 这个一阶差分方程的矩阵表达式为a(n) = Aa(n T)其中利用迭代方法求数值解，也就是按时间步长法仿真种群增长的动态过程， 模拟幼兔和成兔占整体比例随时间的变化。
>> a=[01;1 1];x=[10]';for k=2:12 y=a*x(:,k-1); x=[x y]; endzz=repmat(sum(x),[2 1]); z=x./zz; t=1:12;>> plot(t,x(1,:),'r*',t,x(2,:),'b*'),grid; >> plot(t,z(1,:),'r*',t,z(2,:),'b*'),grid;由数值模拟结果可见，兔子数量递增，但是幼兔和成兔在种群中所占比例很快会 趋于一个极限。
由线性差分方程的定性理论可知，这个极限就是对应于差分方程 系数矩阵A 主特征值得归一化特征向量。
因为 A 是非负矩阵，由矩阵理论知，A 主特征值是正实数,是最大的特征值。
>> [v,d]=eig(a)v-0.85070.52570.52570.8507d-0.61800 1.6180>>max(max(d))ans1.6180a(n)=勺。
(n)]®(n)丿，01Q JL 畔 11 a IIFit Edit ViewInsert Fix^i D M Id opWinder Help■M□ JI•丄$ia 11 BI ■4 £i SQ -aunt iFite- Edit View Insert Toda DMklap Wirdiwr HetpimrHi>> v(:,2)/sum(v(:,2)) ans0.3820 0.6180由数值计算可知，系数矩阵模 A 最大的特征值r=1.618，生物上称之为种群的内 禀增长率，是个大于一的实数。
因此种群数量随时间递增。
相应的归一化的特征 向量的两个分量0.382和0.618正是幼兔和成兔在种群中所占比例趋近的稳定值， 生物上称之为种群的稳定分布。
从这个例子的讨论可见，数值模拟能帮我们对系统的变化有更直观的认识。
但 是只有通过计算方程组系数矩阵的特征值和特征向量，运用差分方程定性分析才 能对解得渐进性质给出确定的结论。
一、微分方程模型1. 实验题目蓝鲸的內禀增长率每年估计为5%,估计蓝鲸的最大环境载量为150000条，磷 虾是蓝鲸喜欢的一种食物。
磷虾的最大饱和种群为500吨/英亩*。
当缺少捕食者， 环境不拥挤时，磷虾种群以每年 25%的速率增长。
磷虾500吨/英亩可以提高蓝 鲸2%的年增长率，同时150000条蓝鲸将减少磷虾10%的年增长率。
（1） 组建一个蓝鲸和磷虾的动态模型，模拟两个种群随时间的变化。
假设初始 状态为蓝鲸5000条，磷虾750吨/英亩；（2） 确定蓝鲸与磷虾是否可以长期共存；（3） 假设捕捞使得蓝鲸只剩下它的平衡态的 5%，而磷虾保持平衡态的数量。
描 述一旦停止捕捞将发生什么情况。
蓝鲸恢复需要多长时间？磷虾群体将发生什么 变化？进一步，给出蓝鲸种群恢复时间对它所受伤害程度的依赖关系。
2. 实验内容解（1）假设（a ） 蓝鲸和磷虾单独存在时，两种群都依靠有限的自然资源增长，即遵循 Logistic 模型； （b ） 蓝鲸捕食磷虾，蓝鲸平均增长率的增加量正比于单位区域内磷虾数量，磷 虾平均增长率的减少量正比于蓝鲸数量记N （t ）为蓝鲸在t 时刻的种群数量（条），2（t ）为磷虾在t 时刻的种群质量 （吨/英亩），于是，依据假设，建立蓝鲸和磷虾两个种群的动态模型+ a ?1 N 1N 2其中「1=0.05和r2=0.25分别表示蓝鲸和磷虾种群的内禀增长率，^=150000（条）N 1 W%) MN , dt dN 2 dt和k2=500 (吨/英亩)分别表示蓝鲸和磷虾种群环境载量， ai2=0.02/500表示每 英亩每吨磷虾对蓝鲸平均增长率的改变， a2i=0.10/150000表示每条蓝鲸对磷虾 平均增长率的改变。
下一步，为了便于数值模拟，保持数量级的协调，把鲸鱼的单位改为以千条为 单位。
当初始状态为蓝鲸5千条、磷虾750吨/英亩时，动态模型的数值模拟MATLAB!令为>> vG=@(t,y)[0.05*y(1)*(1-y(1)/150)+(0.02/500)*y (2) *y(1);0.25* y(2) *(1-y( 2)/500)-(0.1/150)*y(1)*y(2)];>> [t,y]=ode45(vG,[0,150],[5,750]);plot(t,y(:,1),'-'),grid on,hold on >> plot(t,y(:,2),'-.'),grid on,xlabel(' 时间'),ylabel(' 种群数量');>> gtext('实线表示蓝鲸，虚线表示磷虾');由图可见，蓝鲸的数量随着时间的增加而逐渐增加， 磷虾的数量随时间的增加而逐渐减少，最后都分别趋近一个稳定值。
(2)由常微分方程理论知，方程组的常数值解为方程的平衡点，由平衡点的稳 定性确定了方程解的动态性质，即解的渐近行为。
上一问的数值结果表明，方程 组具有一个稳定的正平衡点，也就是说蓝鲸与磷虾可以长期共存。
首先求方程组的平衡点，Matlab 指令为：>> syms y1 y2>> f1=0.05*y1*(1-y1/150)+(0.02/500)*y2*y1; >> f2=0.25*y2*(1-y2/500)-(0.1/150)*y1*y2; >> [y1steady,y2steady]=solve(f1,f2);<dN 1N 10.02— 1= 0.052(1—dt 150 500 dN 2N 2 0.1= 0.252(1 - 2) + dt 500 150□ (0)=5, 2(0) = 750N 2N ,N 1N 2>> disp('平衡点是：’)；平衡点是：>> disp([vpa(y1steady,6) vpa(y2steady,6)]);[0, 0][150.0, 0][ 0, 500.0][181.034, 258.621]接着，考虑方程组在平衡点(N,N2) =(Xi,y i),i=1,2,3,4 附近的局部线性方程零解的稳定性du0.1 0.02 0.02—二(0.05 - X i +------ yju X i Vdt 150 500 500dv 0.5 0.1 、0.1一二(0.25 - y K )v y i udt 500 150 150这些线性方程组零解的稳定性由其系数矩阵的特征值确定。
利用Matlab指令求系数矩阵的特征值>> x=[0,150,0,181.034];y=[0,0,500,258.621];>> for i=1:4A=[0.05-0.1*x(i)/150+0.02*y(i)/500,0.02*x(i)/500;-0.1*y(i)/150 0.25- 0.5*y(i)/500-0.1*x(i)/150,];b=eig(A);disp([b(1) b(2)])end0.0500 0.2500-0.0500 0.1500-0.2500 0.0700-0.0948 + 0.0077i -0.0948 0.0077i得到的结果表明，在正平衡点(181.034, 258.621),相应的两个特征值的实部都 是负的。
按照微分方程定性理论可知，方程组正平衡点 (181.034, 258.621)是渐近稳定的，即从任意初值点出发，解轨线都会趋近该点。
因此，可以断言，只要 停止捕捞，蓝鲸与磷虾可以长期共存。
(3)为了确定当蓝鲸数量为平衡态数量的r%，磷虾数量为平衡态时，停止捕捞，蓝鲸恢复到平衡态需要的时间，只有通过数值模拟。
对不同的初值N1(0)=r/100X 181.034,2(0)=258.621在一个充分长的时间区间上求方程组的数值解，注意到蓝鲸数量会递增趋于平衡态，可以N1(T) > 181.034-0.001确定的时间T 近似表示 种群恢复所需的时间，得到对应的函数关系T=T(r).Matlab 指令如下：G=@(t,N)[0.05*N(1)*(1-N(1)/150)+(0.02/500)* N(2) *N(1);0.25*N(2)*(1-N(2 )/500)-(0.1/150)*N(1)* N(2)]; T=[]；for r=0.05:0.05:0.9[t,N]=ode45(G,[0,200],[181.034*r,258.621]); subplot(1,3,1),plot(t,N(:,1),'-',t,N(:,2),'-.'),xlabel(' 时间'),ylabel(' 种群数量'),grid on,holdon d=fi nd(N(:,1)>181.034-0.001,1);T=[T t(d)];end>> subplot(1,3,2),plot(0.05:0.05:0.9,T),xlabel(' 损伤程度r'),ylabel('恢复时间 T'),grid>> gtext('实线表示蓝鲸，虚线表示磷虾')由图可知，得到的函数关系 尽管ode45采用了四阶、五阶龙格-库塔法，可能是离散的时间步长太大，计算 得结果并未反映连续系统的真实规律。