实验六 机理模型与平衡原理Word文档下载推荐.docx

实验六 机理模型与平衡原理Word文档下载推荐.docx

《实验六 机理模型与平衡原理Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《实验六 机理模型与平衡原理Word文档下载推荐.docx（14页珍藏版）》请在冰豆网上搜索。

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,:

t,z（2,:

由数值模拟结果可见，兔子数量递增，但是幼兔和成兔在种群中所占比例很快会趋于一个极限。

由线性差分方程的定性理论可知，这个极限就是对应于差分方程系数矩阵A主特征值得归一化特征向量。

因为A是非负矩阵，由矩阵理论知，A主特征值是正实数，是最大的特征值。

[v,d]=eig（a）

v 

=

-0.8507 

0.5257

0.5257 

0.8507

d 

-0.6180 

0 

1.6180

max（max（d））

ans 

v（:

2）/sum（v（:

2））

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）蓝鲸捕食磷虾，蓝鲸平均增长率的增加量正比于单位区域内磷虾数量，磷虾平均增长率的减少量正比于蓝鲸数量

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

其中r1=0.05和r2=0.25分别表示蓝鲸和磷虾种群的内禀增长率，k1=150000（条）和k2=500（吨/英亩）分别表示蓝鲸和磷虾种群环境载量，a12=0.02/500表示每英亩每吨磷虾对蓝鲸平均增长率的改变，a21=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）,'

-'

）,gridon,holdon

plot（t,y（:

2）,'

-.'

）,gridon,xlabel（'

时间'

）,ylabel（'

种群数量'

）;

gtext（'

实线表示蓝鲸，虚线表示磷虾'

由图可见，蓝鲸的数量随着时间的增加而逐渐增加，磷虾的数量随时间的增加而逐渐减少，最后都分别趋近一个稳定值。

（2）由常微分方程理论知，方程组的常数值解为方程的平衡点，由平衡点的稳定性确定了方程解的动态性质，即解的渐近行为。

上一问的数值结果表明，方程组具有一个稳定的正平衡点，也就是说蓝鲸与磷虾可以长期共存。

首先求方程组的平衡点，Matlab指令为：

symsy1y2

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]

接着，考虑方程组在平衡点（N1,N2）=（xi,yi）,i=1,2,3,4附近的局部线性方程零解的稳定性

这些线性方程组零解的稳定性由其系数矩阵的特征值确定。

利用Matlab指令求系数矩阵的特征值

x=[0,150,0,181.034];

y=[0,0,500,258.621];

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）]）

0.0500 

0.2500

-0.0500 

0.1500

-0.2500 

0.0700

-0.0948 

+ 

0.0077i 

- 

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=[];

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（:

t,N（:

）,xlabel（'

）,grid 

on,hold 

on

d=find（N（:

1）>

181.034-0.001,1）;

T=[T 

t（d）];

subplot（1,3,2）,plot（0.05:

0.9,T）,xlabel（'

损伤程度r'

恢复时间T'

gtext（'

）

由图可知，得到的函数关系T=T（r）是不光滑的。

注意到计算机只会计算离散值，尽管ode45采用了四阶、五阶龙格-库塔法，可能是离散的时间步长太大，计算得结果并未反映连续系统的真实规律。

重新规定步长，计算量大些，但能够保证离散结果逼近连续情形的精度。

tspan=linspace（0,150,1000）;

[t,N]=ode45（G,tspan,[181.034*r,258.621]）;

subplot（1,3,3）,plot（0.05:

）,grid

果然得到理想的效果。

在利用计算机进行模型的数值求解时，对模型解的性质先有一个定性分析是非常重要的，以确保离散网格点上的解确实保持原连续系统的性质。

三、离散与连续

上节的例子说明连续系统仿真必须限制离散步长，否则可能会产生不同的结果，甚至导致复杂的动力学行为，这些行为并不是连续系统所具有的。

我们在借助数值模拟研究连续系统时，计算机只会在离散点上计算，不同的离散格式同样会导致不同的动力学行为，有些行为并不是连续系统所具有的。

我们通过下面的例子说明这一点。

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

解：

先考察连续Logistic模型

的动力学行为。

运用微分方程定性理论分析，得知对任意给定的不同参数r>

0,该方程从任意初值出发的轨迹线都将趋近于平衡点N=1.运用数值模拟方法,可以得到图:

Matlab指令为：

r=[1 

3.6];

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 

n=16;

t=linspace（0,10,n）;

rr=r（i）;

x1=dislog1（rr,0.1,n）;

subplot（2,2,2）,axis（[0,10,0.1,1.5]）,plot（t,x1,'

--k'

x2=dislog2（rr,0.1,n）;

subplot（2,2,3）,axis（[0,10,0.1,1.5]）,plot（t,x2,'

x3=dislog3（rr,0.1,n）;

subplot（2,2,4）,axis（[0,10,0.1,1.5]）,plot（t,x3,'

end

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

文件一：

functionN=dislog1（r,N0,n）

N=[N0zeros（1,n-1）];

x=[（exp（r*10/n）-1）*N0zeros（1,n-1）];

fori=2:

n

x（i）=exp（r*10/n）*x（i-1）/（1+x（i-1））;

N（i）=x（i）/（exp（r*10/n）-1）;

文件二：

functionN=dislog2（r,N0,n）

x=[r*10/n*N0zeros（1,n-1）];

x（i）=exp（r*10/n）*x（i-1）*exp（-x（i-1））;

N（i）=x（i）/（r*10/n）;

文件三：

functionN=dislog3（r,N0,n）

x=[r*10/n*N0/（r*10/n+1）zeros（1,n-1）];

x（i）=（r*10/n+1）*x（i-1）*（1-x（i-1））;

N（i）=x（i）*（r*10/n+1）/（r*10/n）;

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

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

在图中，每个子图横轴变量为t，纵轴变量为N=N（t）.

（1）如果在区间[t，t+Δt]上，对

积分，分离变量得

由此得到

记

则得到差分方程

因为在转化的过程中没有附加任何限制，这个离散模型仅仅描述了连续的Logistic模型在给定的离散时间点上的动态，所以它应该具有与连续模型同样的动态行为。

数值实验验证了这一点，即右上图。

（2）假设在[t，t+Δt]上，单位变化率保持不变，即

在[t，t+Δt]上积分，得到

得到微分方程

当Δt=1时，称为Ricker模型。

该模型常被用于描述鱼群的变化。

由于这个离散模型在小区间上部分改变了连续模型，当参数r变化时，它会表现出不同于连续系统的性质，如左下图。

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

得到差分方程

当Δ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模型

运用差分方程定性理论做些初步的分析，得到这个差分方程的两个平衡点，

当0<

μ<

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

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

当μ>

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

当1<

3时，∣2-μ∣<

1,随着t→∞，u（t）→0,从而x（t）→1-1/μ，所以非零平衡解X*=1-1/μ是渐近稳定的。

因为μ=rΔt+1，仅当0<

rΔt<

2时，差分方程解与原来微分方程解具有相同的性质。

当3<

μ时，X*=1-1/μ不再是稳定的平衡解，进一步的理论分析可以参考有关混沌动力学的教材，下面通过数值模拟，给出参数μ在[0,4]中取值时，相应差分方程动态性质的直观描述为：

程序：

rate=[0.81.122.753.253.53.553.75];

fork=1:

8

r=rate（k）;

x=zeros（1,101）;

x

（1）=0.1;

fori=1:

100

x（i+1）=r*x（i）*（1-x（i））;

subplot（4,2,k）,plot（0:

100,x,'

.'

End

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

程序为：

mup=1:

0.001:

4;

m=length（mup）;

np=zeros（m,50）;

p=1;

mu=1:

n=0.6;

200

n=mu*（1-n）*n;

50

np（p,i）=n;

p=p+1;

plot（mup,np（:

i）,'

'

MarkerSize'

0.5）;

hold 

on;