实例2 动物种群的相互竞争与相互依存的模型
在生物的种群关系中,一种生物以另一种生物为食的现象,称为捕食.一般说来,由于捕食关系,当捕食动物数量增长时,被捕食动物数量就逐渐下降,捕食动物由于食物来源短缺,数量也随之下降,而被捕食动物数量却随之上升.这样周而复始,捕食动物与被捕食动物的数量随时间变化形成周期性的震荡.
田鼠及其天敌的田间种群消长动态规律也是如此.实验调查数据表明:无论是田鼠还是其天敌的数量都呈周期性的变化,天鼠与天敌的作用系统随时间序列推移,田鼠密度逐渐增加,其天敌随之增加,但时间上落后一步.由于天敌密度增加,则田鼠密度降低,而田鼠密度的降低,则其天敌密度亦减少,如此往复循环,从而形成一定的周期.试用数学模型来概括这一现象,并总结出其数量变化的近似公式.
一 问题分析及模型的建立
设x(t)和y(t)分别表示t时刻田鼠与其天敌的数量,如果单独生活,田鼠的增长速度正比于当时的数量,即
dx=λx dt
dy=-μy dt而田鼠的天敌由于没有被捕食对象,其数量减少的速率正比于当时的数量,即
现在田鼠与其天敌生活一起,田鼠一部分遭到其天敌的消灭,于是以一定的速率α减少,减少的数量正比于天敌的数量,因此有
dx=(λ-αy)x dt
类似地,田鼠的天敌有了食物,数量减少的速率μ减少β,减少的量正比于田鼠的数量,因此有
dy=-(μ-βx)y dt
上述公式,最后两个方程联合起来称为Volterra-Lot方程,这里α,β,λ,μ均为正数,初始条件为
x(0)=x0,y(0)=y0
现在通过实验调查所得到的数据如表,此数据为每隔两个月田间调查一次,得到的田鼠及其天敌种群数量的记录,数量的单位经过处理.试建立合理的数学模型.
表 田鼠种群数量记录
29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5
表 田鼠天敌种群数量记录
1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3
二 模型的求解
Volterra-Lotok方程的解析解即x,y的显示解难求出,因此公式的参数方程不宜直接用Matlab函数来拟合解,可用如下的方法来求其近似解.
Volterra-Lotok可转化为
⎧dlnx=(λ-αy)dt ⎨dlny=(-μ+βx)dt⎩
在区间[ti-1,ti]上积分,得
lnxi-lnxi-1=λ(ti-ti-1)-αS1i
lnyi-lnyi-1=-μ(ti-ti-1)+βS2i
这里,S1i=⎰ti
ti-1ydt,S2i=⎰xdt, i=1, ,m ti-2ti
于是得到方程组
⎧A1P1=B1 ⎨ AP=B2⎩22
这里
⎛t1-t0 t-tA1= 21
t-t⎝mm-1-S11⎫⎛t1-t0⎪ -S12⎪ t2-t1A= 2 ⎪ ⎪ t-t-Sim⎪m-1⎭⎝m-S⎫⎪-S22⎪ ⎪⎪-S2m⎪⎭
⎛-μ⎫⎛λ⎫ ⎪ P=P1= 2 β⎪⎪ α⎪⎝⎭⎝⎭
B1=(lnxyx1y, ,lnm)T B=(ln1, ,lnm)T x0xm-1y0ym-1
T-1TA2B2 因此方程组参数的最小二乘解为 T-1T P=(AA)A1B1 P=(A2A2)111
由于x(t)和y(t)均为未知,因此S1i,S2用数值积分方法的梯形公式解
S1i=⎰
⎰titi-1ydt≈ti-ti-1(yi+yi-1) 2 S2=ti
ti-1xdt=ti-ti-1(xi+xi-1) 2
这样就可求得参数的近似值.
模型参数求解的程序为
clear all,clc
X=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...
34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];
Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ...
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3];
N=[X;Y];
T=[0:2:60];
for i=1:30
A(i,1)=T(i+1)-T(i);
A(i,[2 3])=((T(i+1)-T(i))/2)*[-(N(1,i+1)+N(1,i)),-(N(2,i+1)+N(2,i))];
B(i,[1 2])=[log(N(1,i+1)/N(1,i)),log(N(2,i+1)/N(2,i))];
end;
A1=A(:,[1 3]);
P1=inv((A1'*A1))*A1'*B(:,1)
A2=A(:,[1 2]);
P2=inv((A2'*A2))*A2'*B(:,2)
上述结果代入Volterra-Lotok方程,用MATLAB函数ode45求方程在时间[0,60]的数值解.作图可看到田鼠及其天敌数量的周期震荡.
求方程Volterra-Lotok的数值解的程序为
定义函数vlok为
[vlok.m]
function dydt=vlok(T,Y)
dydt=[(0.8765-0.5468*Y(2))*Y(1);(-0.1037+0.0010*Y(1))*Y(2)];
clear all, clc
X=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...
34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];
Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ...
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3]; N=[X,Y];
T=[0:2:60];
[t,Y]=ode45(@vlok,[0:0.5:60],[29.7 1.6]);
plot(t,Y(:,1)/100,'k');
hold on;
plot(t,Y(:,2),'-.k');
title('田鼠及其天敌的Volterra-Lotok模型拟合曲线');
xlabel('时间');
ylabel('数量(只/每百)');
gtext('田鼠');
gtext('天敌');
legend('田鼠','天敌');legend('田鼠','天敌');
图 田鼠及其天敌的模拟曲线
实线和虚线分别为田鼠和天敌的实际值,田鼠的数量为y坐标乘以100.
实例2 动物种群的相互竞争与相互依存的模型
在生物的种群关系中,一种生物以另一种生物为食的现象,称为捕食.一般说来,由于捕食关系,当捕食动物数量增长时,被捕食动物数量就逐渐下降,捕食动物由于食物来源短缺,数量也随之下降,而被捕食动物数量却随之上升.这样周而复始,捕食动物与被捕食动物的数量随时间变化形成周期性的震荡.
田鼠及其天敌的田间种群消长动态规律也是如此.实验调查数据表明:无论是田鼠还是其天敌的数量都呈周期性的变化,天鼠与天敌的作用系统随时间序列推移,田鼠密度逐渐增加,其天敌随之增加,但时间上落后一步.由于天敌密度增加,则田鼠密度降低,而田鼠密度的降低,则其天敌密度亦减少,如此往复循环,从而形成一定的周期.试用数学模型来概括这一现象,并总结出其数量变化的近似公式.
一 问题分析及模型的建立
设x(t)和y(t)分别表示t时刻田鼠与其天敌的数量,如果单独生活,田鼠的增长速度正比于当时的数量,即
dx=λx dt
dy=-μy dt而田鼠的天敌由于没有被捕食对象,其数量减少的速率正比于当时的数量,即
现在田鼠与其天敌生活一起,田鼠一部分遭到其天敌的消灭,于是以一定的速率α减少,减少的数量正比于天敌的数量,因此有
dx=(λ-αy)x dt
类似地,田鼠的天敌有了食物,数量减少的速率μ减少β,减少的量正比于田鼠的数量,因此有
dy=-(μ-βx)y dt
上述公式,最后两个方程联合起来称为Volterra-Lot方程,这里α,β,λ,μ均为正数,初始条件为
x(0)=x0,y(0)=y0
现在通过实验调查所得到的数据如表,此数据为每隔两个月田间调查一次,得到的田鼠及其天敌种群数量的记录,数量的单位经过处理.试建立合理的数学模型.
表 田鼠种群数量记录
29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5
表 田鼠天敌种群数量记录
1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3
二 模型的求解
Volterra-Lotok方程的解析解即x,y的显示解难求出,因此公式的参数方程不宜直接用Matlab函数来拟合解,可用如下的方法来求其近似解.
Volterra-Lotok可转化为
⎧dlnx=(λ-αy)dt ⎨dlny=(-μ+βx)dt⎩
在区间[ti-1,ti]上积分,得
lnxi-lnxi-1=λ(ti-ti-1)-αS1i
lnyi-lnyi-1=-μ(ti-ti-1)+βS2i
这里,S1i=⎰ti
ti-1ydt,S2i=⎰xdt, i=1, ,m ti-2ti
于是得到方程组
⎧A1P1=B1 ⎨ AP=B2⎩22
这里
⎛t1-t0 t-tA1= 21
t-t⎝mm-1-S11⎫⎛t1-t0⎪ -S12⎪ t2-t1A= 2 ⎪ ⎪ t-t-Sim⎪m-1⎭⎝m-S⎫⎪-S22⎪ ⎪⎪-S2m⎪⎭
⎛-μ⎫⎛λ⎫ ⎪ P=P1= 2 β⎪⎪ α⎪⎝⎭⎝⎭
B1=(lnxyx1y, ,lnm)T B=(ln1, ,lnm)T x0xm-1y0ym-1
T-1TA2B2 因此方程组参数的最小二乘解为 T-1T P=(AA)A1B1 P=(A2A2)111
由于x(t)和y(t)均为未知,因此S1i,S2用数值积分方法的梯形公式解
S1i=⎰
⎰titi-1ydt≈ti-ti-1(yi+yi-1) 2 S2=ti
ti-1xdt=ti-ti-1(xi+xi-1) 2
这样就可求得参数的近似值.
模型参数求解的程序为
clear all,clc
X=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...
34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];
Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ...
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3];
N=[X;Y];
T=[0:2:60];
for i=1:30
A(i,1)=T(i+1)-T(i);
A(i,[2 3])=((T(i+1)-T(i))/2)*[-(N(1,i+1)+N(1,i)),-(N(2,i+1)+N(2,i))];
B(i,[1 2])=[log(N(1,i+1)/N(1,i)),log(N(2,i+1)/N(2,i))];
end;
A1=A(:,[1 3]);
P1=inv((A1'*A1))*A1'*B(:,1)
A2=A(:,[1 2]);
P2=inv((A2'*A2))*A2'*B(:,2)
上述结果代入Volterra-Lotok方程,用MATLAB函数ode45求方程在时间[0,60]的数值解.作图可看到田鼠及其天敌数量的周期震荡.
求方程Volterra-Lotok的数值解的程序为
定义函数vlok为
[vlok.m]
function dydt=vlok(T,Y)
dydt=[(0.8765-0.5468*Y(2))*Y(1);(-0.1037+0.0010*Y(1))*Y(2)];
clear all, clc
X=[29.7 33.1 32.5 69.1 134.2 236.0 269.6 162.2 69.6 39.8 ...
34.0 20.7 22.0 37.6 57.6 124.6 225.0 272.7 195.7 94.5 41.9 25.7 ... 10.9 22.5 33.5 48.2 92.5 183.3 268.5 230.6 115.5];
Y=[1.6 1.3 1.1 1.2 1.1 1.3 1.8 2.2 2.4 2.2 1.9 1.5 1.5 1.2 0.9 ...
1.1 1.3 1.6 2.3 2.4 2.2 1.7 1.8 1.5 1.2 1.0 0.9 1.1 1.3 1.9 2.3]; N=[X,Y];
T=[0:2:60];
[t,Y]=ode45(@vlok,[0:0.5:60],[29.7 1.6]);
plot(t,Y(:,1)/100,'k');
hold on;
plot(t,Y(:,2),'-.k');
title('田鼠及其天敌的Volterra-Lotok模型拟合曲线');
xlabel('时间');
ylabel('数量(只/每百)');
gtext('田鼠');
gtext('天敌');
legend('田鼠','天敌');legend('田鼠','天敌');
图 田鼠及其天敌的模拟曲线
实线和虚线分别为田鼠和天敌的实际值,田鼠的数量为y坐标乘以100.