MATLAB 数学应用 微分方程 常微分方程 求解捕食者-猎物方程

本文说明如何使用 ode23 和 ode45 求解表示捕食者/猎物模型的微分方程。这两个函数用于对使用可变步长大小 Runge-Kutta 积分方法的常微分方程求数值解。ode23 使用一对简单的 2 阶和 3 阶公式实现中等精度，ode45 使用一对 4 阶和 5 阶公式实现更高的精度。

以名为 Lotka-Volterra 方程，也即捕食者-猎物模型的一对一阶常微分方程为例：

$\frac{dx}{dt}=x-\alpha xy$
$\frac{dy}{dt}=-y+\beta xy.$

变量 x 和 y 分别计算猎物和捕食者的数量。二次交叉项表示物种之间的交叉。当没有捕食者时，猎物数量将增加，当猎物匮乏时，捕食者数量将减少。

编写方程代码

为了模拟系统，需要创建一个函数，以返回给定状态和时间值时的状态导数的列向量。在 MATLAB 中，两个变量 x 和 y 可以表示为向量 y 中的前两个值。同样，导数是向量 yp 中的前两个值。函数必须接受 t 和 y 的值，并在 yp 中返回公式生成的值。

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

在此示例中，公式包含在名为 lotka.m 的文件中。此文件使用 α=0.01 和 β=0.02 的参数值。

type lotka % 加载对应lotka.m文件内容

function yp = lotka(t,y)
yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

模拟系统

使用 ode23 在区间 0<t<15 中求解 lotka 中定义的微分方程。使用初始条件 x(0)=y(0)=20，使捕食者和猎物的数量相等。

t0 = 0;
tfinal = 15;
y0 = [20; 20];   
[t,y] = ode23(@lotka,[t0 tfinal],y0);

绘制结果

绘制两个种群数量对时间的图。

plot(t,y)
title('Predator/Prey Populations Over Time')
xlabel('t')
ylabel('Population')
legend('Prey','Predators','Location','North')

现在绘制两个种群数量的相对关系图。生成的相平面图非常清晰地表明了二者数量之间的循环关系。

plot(y(:,1),y(:,2))
title('Phase Plane Plot')
xlabel('Prey Population')
ylabel('Predator Population')

比较不同求解器的结果

使用 ode45 而不是 ode23 再次求解该方程组。ode45 求解器的每一步都需要更长的时间，但它的步长也更大。然而，ode45 的输出是平滑的，因为默认情况下，此求解器使用连续展开公式在每个步长范围内的四个等间距时间点生成输出。（您可以使用 ‘Refine’ 选项调整时间点数。）绘制两个解进行比较。

[T,Y] = ode45(@lotka,[t0 tfinal],y0);

plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'-');
title('Phase Plane Plot')
legend('ode23','ode45')

结果表明，使用不同的数值方法求解微分方程会产生略微不同的答案。