掠食者—猎物方程,MATLAB

掠食者—猎物方程,MATLAB
考虑一个简单的生态系统,包含无限食物供给的兔群和依靠捕食兔子为生的狐狸群体.这个系统可以通过一对非线性一阶微分方程来建模：
dr/dt=2r-arf,r(0)=r0
df/dt=-f=arf,f(0)=f0
其中t为时间,r(t)为兔子的数量,f(t)为狐狸的数量,a为一个正常数.如果a>0,这两个种群没有关系,兔子很好的生存,而狐狸因
饥饿而死亡.如果a>0,狐狸捕食和自身数量成比例的兔子,这种捕食会造成兔子数量的减少,和狐狸数量的增加（原因不像减少那么显而易见）.
这个非线性系统的解无法写成已知的其他函数的组合,只能通过数值求解.可以证明该系统的解有周期性,其周期取决于初始条件.也就是说,对任意的r(0)和f(0)值,都存在一个时间t=tp,此时这两个种的数量都等于初始值.因些对所的和t有
r(t+tp)=r(t),f(t+tp)=f(t)
1.计算系统一思想的解,取r0=300,f0=150以及a=0.01.结果应该有tp接近期.绘制两幅,一个以r和f做为t的函数,一个相位平面图分别以r和t为两个坐标轴.
2.计算并绘制r0=15,f0=22和a=0.01时的解.tp应该接近6.62.
3.计算并绘制r0=102,f0=198和a=0.01时的解.并计算出周期tp,可以通过误差试算或者消息句柄机制.
4.点(r0,f0)=(1/a,2/a)是一个稳定的平衡点,如果初值为此值,那么种群数量不会变化.如果初值不等于但比较接近该平衡点,那么
其数量不会发生在的变化.令u(t)=r(t)-1/a,v(t)=f(t)-2/a,函数u(t)和v(t)满足另外一个非线性微分方程系统,但如果忽
略项,系统变为线性.试问这个线性系统为何?其周期解的周期是多少?
最好有截图

说明几点：
1、第二个方程错误,应为df/dt=-f+arf
2、周期可以放大曲线后观察两个峰值间的距离
3、这道题里面错误太多,除上面的公式错误外,还有
3.1“如果a>0,这两个种群没有关系”==> 这里应为a=0
3.2“对所的和t”==> “和”字去掉
3.3“相位平面图分别以r和t为两个坐标轴”==> 应该是“r和f”
3.4 下面这几句,完全就猜不出什么意思：“计算系统一思想的解”、“结果应该有tp接近期”、
“通过误差试算或者消息句柄机制”、“忽略项”,等等
4、希望以后提问能稍微认真一点,是对别人的尊重,也有助于你的问题尽快得到解决

 
给几个截图,代码见后.



% 1.取r0=300,f0=150以及a=0.01
% 纠错：从形式上判断,第二个方程应为df/dt=-f+arf,仿真结果也证实了这个判断
Y0=[300; 150];
dY=inline('[2*Y(1)-0.01*Y(1)*Y(2); -Y(2)+0.01*Y(1)*Y(2)]','t','Y');
[t,Y] = ode45(dY,10,Y0,[],a);
figure(101)
plot(t,Y)
legend('Rabbit','Fox');
xlabel('Time')
figure(102)
r=Y(:,1);
f=Y(:,2);
plot(r,f)
xlabel('r')
ylabel('f')
 
% 2.计算并绘制r0=15,f0=22和a=0.01时的解.tp应该接近6.62
Y0=[15; 22];
dY=inline('[2*Y(1)-0.01*Y(1)*Y(2); -Y(2)+0.01*Y(1)*Y(2)]','t','Y');
[t,Y] = ode45(dY,10,Y0,[],a);
figure(103)
plot(t,Y)
legend('Rabbit','Fox');
xlabel('Time')
 
% 3.计算并绘制r0=102,f0=198和a=0.01时的解.并计算出周期tp
% 从图中观察容易看到,周期约为4.5（放大后看两个峰值间的距离）
Y0=[102; 198];
dY=inline('[2*Y(1)-0.01*Y(1)*Y(2); -Y(2)+0.01*Y(1)*Y(2)]','t','Y');
[t,Y] = ode45(dY,10,Y0,[],a);
figure(104)
plot(t,Y)
legend('Rabbit','Fox');
xlabel('Time')
 
% 4.验证：点(r0,f0)=(1/a,2/a)是一个稳定的平衡点
Y0=[100; 200];
dY=inline('[2*Y(1)-0.01*Y(1)*Y(2); -Y(2)+0.01*Y(1)*Y(2)]','t','Y');
[t,Y] = ode45(dY,10,Y0,[],a);
figure(105)
plot(t,Y)
legend('Rabbit','Fox');
xlabel('Time')
 
% 求新的微分方程.使用符号运算推导（手工推导亦可）,得到新的微分方程为:
%    du/dt = -auv-v; dv/dt = auv + 2u
% 但题中错字比较多,有些地方实在猜不到是什么意思,比如这个“忽略项”
syms r f u v a
dr=2*r-a*r*f;
df=-f+a*r*f;
du=simple(subs(dr,{r,f},{u+1/a,v+2/a}))
dv=simple(subs(df,{r,f},{u+1/a,v+2/a}))