西农 建模实验六

实验六 机理模型与平衡原理

实验目的

如果对所研究的问题了解的比较深入，知道产生现象的内在的机理，那么依

据机理建模，则模型具有更好的可靠性和广泛性。不考虑随机因素，假设每一时刻是确定的如果对系统状态的观测和描述只在离散的时间点上，则构成差分方程模型；如果考虑系统随时间连续变化，则是微分方程模型。本节主要以这两类方程为例，介绍用MATLAB 软件求解机理模型的基本方法。

差分方程模型

一、实验题目

由一对兔子开始，一年可以繁殖出多少只兔子？如果一对兔子每个月可以生

一对小兔子，兔子在出生两个月后就具有繁殖能力，由一对刚出生一个月的兔子开始，一年内兔子种群数量如何变化。求这个种群的稳定分布和固有增长率。

二、实验内容

解 假设

（a ）兔子每经过一个月底就增加一个月龄； （b ）月龄大于等于2的兔子都具有繁殖能力；

（c ）具有繁殖能力的兔子每一个月一定生一对兔子； （d ）兔子不离开群体（不考虑死亡）

记第n 个月初的幼兔（一月龄兔）数量为a 0(n ）,成兔（月龄大于等于2）数量为a 1(n ），则兔子总数为a(n)= a 0(n ）+a 1(n ）,平衡关系为：

?

??

+==上月初幼兔数量上月初成兔数量本月初成兔数量上月初成兔数量本月初幼兔数量 建立模型：

??

?

??==-+-=-=0)1(,1)1()1()1()()1()(1010110a a n a n a n a n a n a 这个一阶差分方程的矩阵表达式为

)1()(-=n Aa n a

其中

??

? ??=)()()(10n a n a n a ， ?

?? ??=1110A

利用迭代方法求数值解，也就是按时间步长法仿真种群增长的动态过程，模拟幼兔和成兔占整体比例随时间的变化。

>> a=[0 1;1 1];

x=[1 0]';

for k=2:12

y=a*x(:,k-1);

x=[x y];

end

zz=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.8507 0.5257

0.5257 0.8507

d =

-0.6180 0

0 1.6180

>> max(max(d))

ans =

1.6180

>> v(:,2)/sum(v(:,2))

ans =

0.3820 0.6180

由数值计算可知，系数矩阵模A 最大的特征值r=1.618，生物上称之为种群的内禀增长率，是个大于一的实数。因此种群数量随时间递增。相应的归一化的特征向量的两个分量0.382和0.618正是幼兔和成兔在种群中所占比例趋近的稳定值，生物上称之为种群的稳定分布。 从这个例子的讨论可见，数值模拟能帮我们对系统的变化有更直观的认识。但是只有通过计算方程组系数矩阵的特征值和特征向量，运用差分方程定性分析才能对解得渐进性质给出确定的结论。

微分方程模型

一、实验题目

蓝鲸的內禀增长率每年估计为5%，估计蓝鲸的最大环境载量为150000条，磷

虾是蓝鲸喜欢的一种食物。磷虾的最大饱和种群为500吨/英亩*

。当缺少捕食者，环境不拥挤时，磷虾种群以每年25%的速率增长。磷虾500吨/英亩可以提高蓝鲸2%的年增长率，同时150000条蓝鲸将减少磷虾10%的年增长率。

（1）组建一个蓝鲸和磷虾的动态模型，模拟两个种群随时间的变化。假设初始状态为蓝鲸5000条，磷虾750吨/英亩； （2）确定蓝鲸与磷虾是否可以长期共存；

（3）假设捕捞使得蓝鲸只剩下它的平衡态的5%，而磷虾保持平衡态的数量。描述一旦停止捕捞将发生什么情况。蓝鲸恢复需要多长时间？磷虾群体将发生什么变化？进一步，给出蓝鲸种群恢复时间对它所受伤害程度的依赖关系。

二、实验内容

解 ：

（1）假设

（a ）蓝鲸和磷虾单独存在时，两种群都依靠有限的自然资源增长，即遵循Logistic 模型；

（b ）蓝鲸捕食磷虾，蓝鲸平均增长率的增加量正比于单位区域内磷虾数量，磷虾平均增长率的减少量正比于蓝鲸数量

记N 1(t)为蓝鲸在t 时刻的种群数量（条），N 2(t)为磷虾在t 时刻的种群质量（吨/英亩），于是，依据假设，建立蓝鲸和磷虾两个种群的动态模型

???????+-=+-=2

1212222212121

1111

)1()1(N N a k N N r dt

dN N N a k N N r dt dN 其中r 1=0.05和r 2=0.25分别表示蓝鲸和磷虾种群的内禀增长率，k 1=150000(条）

和k 2=500（吨/英亩）分别表示蓝鲸和磷虾种群环境载量，a 12=0.02/500表示每英亩每吨磷虾对蓝鲸平均增长率的改变，a 21=0.10/150000表示每条蓝鲸对磷虾平均增长率的改变。

下一步，为了便于数值模拟，保持数量级的协调，把鲸鱼的单位改为以千条为单位。当初始状态为蓝鲸5千条、磷虾750吨/英亩时，动态模型

??

??

????

?==+-=+-=750)0(,5)0(1501.0)5001(25.050002.0)1501(05.02121222

12111

N N N N N N dt

dN N N N N dt dN 的数值模拟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); >> disp('平衡点是：'); 平衡点是：

>> disp([vpa(y1steady,6) vpa(y2steady,6)]); [ 0, 0] [ 150.0, 0] [ 0, 500.0] [ 181.034, 258.621]

接着，考虑方程组在平衡点（N 1,N 2）=(x i ,y i ),i=1,2,3,4附近的局部线性方程零解的稳定性

?????++-=++-=u

y v x y dt

dv v x u y x dt du

i i i i i i 1501.0)1501.05005.025.0(50002.0)50002.01501.005.0( 这些线性方程组零解的稳定性由其系数矩阵的特征值确定。 利用Matlab 指令求系数矩阵的特征值

>> x=[0,150,0,181.034];y=[0,0,500,258.621]; >> for i=1:4

A=[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)]) end

0.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/100×181.034,N2(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,hold on

d=find(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('实线表示蓝鲸，虚线表示磷虾')

由图可知，得到的函数关系T=T(r)是不光滑的。注意到计算机只会计算离散值，尽管ode45采用了四阶、五阶龙格-库塔法，可能是离散的时间步长太大，计算得结果并未反映连续系统的真实规律。重新规定步长，计算量大些，但能够保证离散结果逼近连续情形的精度。

>> tspan=linspace(0,150,1000);

T=[];

for r=0.05:0.05:0.9

[t,N]=ode45(G,tspan,[181.034*r,258.621]);

d=find(N(:,1)>181.034-0.001,1);T=[T t(d)];

end

>> subplot(1,3,3),plot(0.05:0.05:0.9,T),xlabel('损伤程度

r'),ylabel('恢复时间T'),grid

果然得到理想的效果。在利用计算机进行模型的数值求解时，对模型解的性质先

有一个定性分析是非常重要的，以确保离散网格点上的解确实保持原连续系统的性质。

离散与连续

一、实验题目

上节的例子说明连续系统仿真必须限制离散步长，否则可能会产生不同的结果，甚至导致复杂的动力学行为，这些行为并不是连续系统所具有的。我们在借助数值模拟研究连续系统时，计算机只会在离散点上计算，不同的离散格式同样会导致不同的动力学行为，有些行为并不是连续系统所具有的。我们通过下面的例子说明这一点。

利用数值模拟，讨论离散的和连续的Logistic模型的动力学性质。

二、实验内容

解：

先考察连续Logistic模型

dt

=

dN-

rN

)

1(

/N

的动力学行为。运用微分方程定性理论分析，得知对任意给定的不同参数r>0,该方程从任意初值出发的轨迹线都将趋近于平衡点N=1.运用数值模拟方法,可以得到图:

Matlab指令为：

>> r=[1 3.6];

for i=1:2

funlog=@(t,x)r(i)*x*(1-x);

[tt,xx]=ode45(funlog,[0,10],0.1);

subplot(2,2,1),axis([0,10,0.1,1.5]),plot(tt,xx,'k'),hold on

end

>> n=16;t=linspace(0,10,n);

for i=1:2

rr=r(i);

x1=dislog1(rr,0.1,n);

subplot(2,2,2),axis([0,10,0.1,1.5]),plot(t,x1,'--k'),hold on x2=dislog2(rr,0.1,n);

subplot(2,2,3),axis([0,10,0.1,1.5]),plot(t,x2,'k'),hold on x3=dislog3(rr,0.1,n);

subplot(2,2,4),axis([0,10,0.1,1.5]),plot(t,x3,'k'),hold on end

求解子程序的函数文件为： 文件一：

function N=dislog1(r,N0,n)

N=[N0 zeros(1,n-1)];x=[(exp(r*10/n)-1)*N0 zeros(1,n-1)]; for i=2:n

x(i)=exp(r*10/n)*x(i-1)/(1+x(i-1)); N(i)=x(i)/(exp(r*10/n)-1); end

文件二：

function N=dislog2(r,N0,n)

N=[N0 zeros(1,n-1)];x=[r*10/n*N0 zeros(1,n-1)]; for i=2:n

x(i)=exp(r*10/n)*x(i-1)*exp(-x(i-1)); N(i)=x(i)/(r*10/n); end

文件三：

function N=dislog3(r,N0,n)

N=[N0 zeros(1,n-1)];x=[r*10/n*N0/(r*10/n+1) zeros(1,n-1)]; for i=2:n

x(i)=(r*10/n+1)*x(i-1)*(1-x(i-1)); N(i)=x(i)*(r*10/n+1)/(r*10/n); end

由左上图可以看到，对不同的参数值r=1，r=3.6，微分方程解轨线的变化特征是一样的，都会从原点附近离开，单调增加趋于1.

接着，以三种不同格式将微分方程离散化，为了比较离散方程的动力学行为，统一取离散步长为：

16

10=?t

在图中，每个子图横轴变量为t ，纵轴变量为N=N （t ）. （1）如果在区间[t ，t+Δt]上，对

)1(N rN dt

dN

-= 积分，分离变量得

rdt N

dN

dt dN =-+1 由此得到

t r e t N t t N t N t t N ?=-?+-?+)

(1)

(1/)()(

记

)()1()(,t N t X e t r -==??

λλ

则得到差分方程

)()

(1)(t X t X t t X +=

?+λ

因为在转化的过程中没有附加任何限制，这个离散模型仅仅描述了连续的Logistic 模型在给定的离散时间点上的动态，所以它应该具有与连续模型同样的动态行为。数值实验验证了这一点，即右上图。 （2）假设在[t ，t+Δt]上，单位变化率保持不变，即

]),[,)(1(1

t t t r dr t N r dN N

?+∈-= 在[t ，t+Δt]上积分，得到

))(1()()(t N t r e t N t t N -?=?+

记

)()(,t tN r t X e t r ?==??

λ

得到微分方程

)

()()(t X e

t X t t X -=?+λ

当Δt=1 时，称为Ricker 模型。该模型常被用于描述鱼群的变化。由于这个离散模型在小区间上部分改变了连续模型，当参数r 变化时，它会表现出不同于连续系统的性质，如左下图。

（3）假设在[t ，t+Δt]上，变化率dN/dt 保持不变，即按最简单的欧拉差分格式离散化，

)())(1()()(t N t N tr t N t t N -?+=?+

记

)()

1()(,1t N t r t

r t X t r +??=+?=μ

得到差分方程

))(1)(()(t X t X t t X -=?+μ

当Δt=1时，称为离散的Logistic 模型。这个离散模型在小区间上完全改变了连续模型。当参数r 变化时，它也会表现出不同于连续系统的性质，如上面截图右下图。

在上面的程序中，我们取模拟步数n=16,相当于将模拟时间区间[0,10]分成步长为Δt=10/16的时间间隔，结果对英语较小的参数值r=1，不同形式的离散格式没有导致不同的变化；而对于较大的参数值r=3.6，不同形式的离散格式导致了不同的变化。如果取模拟步数n=100，时间步长缩短为Δt=0.1，则3种离散格式对于不同的参数值r=1和r=3.6都能保持原来系统的动态变化性质。可见只要限制离散步长，采用不同格式的离散化的系统仿真都能保持连续系统的特征。 最后，考虑离散Logistic 模型

))1(1)(1()(---=t X t X t X μ

的动力学行为。运用差分方程定性理论做些初步的分析，得到这个差分方程的两个平衡点，μ/11,0**

-==X X

当0<μ<1时，非负平衡解只有X *=0，在其附近的局部线性化方程为

)1()(-=t u t u μ

此时随着t →∞，u(t)→0，从而x(t)→0，所以零平衡点是渐近稳定的。

当μ>1时，X *=0不再是稳定的平衡解，此时出现正平衡解x *=1-1/μ，在其附近的局部线性化方程为

)1()2()(--=t u t u μ

当1<μ<3时，∣2-μ∣<1,随着t →∞，u(t)→0,从而x(t)→1-1/μ，所以非零平衡解X *=1-1/μ是渐近稳定的。因为μ=r Δt+1，仅当0

Matlab命令为：

>> rate=[0.8 1.1 2 2.75 3.25 3.5 3.55 3.75]; for k=1:8

r=rate(k);

x=zeros(1,101);

x(1)=0.1;

for i=1:100

x(i+1)=r*x(i)*(1-x(i));

end

subplot(4,2,k),plot(0:100,x,'*');

End

随着参数μ在[0,4]的连续变化，离散Logistic模型的解从趋于一个稳定的平衡态，变化为2,4,8周期振荡，还出现混乱无规则游荡，直接画出解的极限态随着参数变化的分岔图为：

Matlab命令为：

>> mup=1:0.001:4;

>> m=length(mup);

>> np=zeros(m,50);

>> p=1;

>> for mu=1:0.001:4

n=0.6;

for i=1:200

n=mu*(1-n)*n;

end

for i=1:50

n=mu*(1-n)*n;

np(p,i)=n;

end

p=p+1;

end

>> for i=1:50

plot(mup,np(:,i),'.','MarkerSize',0.5);hold on;

end