实验六 机理模型与平衡原理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;
升级会员 