使用 SciPy 求解 Lotka-Volterra 模型答案

【问题标题】：Solve Lotka-Volterra model using SciPy使用 SciPy 求解 Lotka-Volterra 模型
【发布时间】：2015-02-04 21:54:22
【问题描述】：

我有以下 Lotka-Volterra 模型

dN1/dt = N1(1-N1-0.7N2)

dN2/dt = N2(1-N2-0.3N1)

其中 N 旁边的 1 和 2 是下标。

我想使用 SciPy 解决这个问题并将结果可视化。我想在 y 轴上用 N2 和在 N1 上用 N1 绘制一个图。如果在第一个等式中将 N1 设置为零，则得到 N2 = 1/0.7，如果在第二个等式中将 N2 设置为零，则得到 N1 = 0.3/1。假设两条线相交。我如何在 Python 中做到这一点？

我在网上阅读了这个tutorial（幻灯片 6 到 16）。这是我目前所拥有的。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

def derivN1(y,t):
    yprime=np.array([1-0.7y[0]])
    return yprime

def derivN2(y,t):
    yprime=np.array([1-0.3y[0]])
    return yprime

start=0
end=1
numsteps=1000
time=np.linspace(start,end,numsteps)
y0=np.array([10])

yN1=integrate.odeint(derivN1,y0,time)
yN2=integrate.odeint(derivN2,y0,time)

plt.plot(time,yN1[:])
plt.plot(time,yN2[:])

但情节不正确。更新：我认为我使用了错误的方法。我正在阅读另一个在线tutorial。我会再解决这个问题。同时，如果有人知道如何解决，请告诉我。

【问题讨论】：

你能描述一下情节是如何不正确的——你的预期与你实际看到的吗？
我想要一个图，其中 N2 在 y 轴上，N1 在 x 轴上。线 N1=0.3 和 N2=1/0.7 应该相交。在我的代码创建的图中，线条不相交，也不清楚哪个轴是哪个。
你描述它的方式听起来很奇怪：这些方程是耦合的，这意味着在每个时间点，你都有一个元组(N1(t), N2(t))。您可以将该元组随时间绘制为隐式参数，但是您只有 一条曲线，而不是两条。如果您想要两条曲线，一条绘制 N1（因变量）与时间（自变量）的关系，那么您可以为 N2 与时间绘制类似的曲线。但是你的 y 轴将代表 N2 或 N1 的值。我不明白你的第二个公式有两个图表。可以添加图片或来源吗？
另外，您确定您不是在绘制两个隐式图，而只是使用不同的起始条件？
查看 scipy 食谱中的 Lotka-Volterra 示例：wiki.scipy.org/Cookbook/…

标签： python numpy ode differential-equations

【解决方案1】：

@WarrenWeckesser 的评论非常好，你应该从那里开始。我只是尝试强调隐式情节和显式情节之间的区别。

首先，设置：

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

time=np.linspace(0,15,5*1024)

def derivN(N, t):
    """Return the derivative of the vector N, which represents
    the tuple (N1, N2). """

    N1, N2  = N
    return np.array([N1*(1 - N1 - .7*N2), N2*(1 - N2 - .3*N1)])

def coupled(time, init, ax):
    """Visualize the system of coupled equations, by passing a timerange and
    initial conditions for the coupled equations.

    The initical condition is the value that (N1, N2) will assume at the first
    timestep. """

    N = integrate.odeint(derivN, init, time)
    ax[0].plot(N[:,0], N[:,1], label='[{:.1f}, {:.1f}]'.format(*init))  # plots N2 vs N1, with time as an implicit parameter
    l1, = ax[1].plot(time, N[:,0], label='[{:.1f}, {:.1f}]'.format(*init))
    ax[1].plot(time, N[:,1], color=l1.get_color())

重要的是要认识到您的方程是耦合的，您应该向odeint 提交一个返回耦合方程导数的函数。由于您有 2 个方程，因此您需要返回一个长度为 2 的数组，每个项表示根据传入变量的导数（在本例中为数组 N(t) = [N1(t), N2(t)]）。

然后您可以完全绘制它，对 N1 和 N2 使用不同的初始条件：

fh, ax = plt.subplots(1,2)
coupled(time, [.3, 1/.7], ax)
coupled(time, [.4, 1/.7], ax)
coupled(time, [1/.7, .3], ax)
coupled(time, [.5, .5], ax)
coupled(time, [.1, .1], ax)
ax[0].legend()
ax[1].legend()
ax[0].set_xlabel('N1')
ax[0].set_ylabel('N2')
ax[1].set_xlabel('time')
ax[1].set_ylabel(r'$N_i$')
ax[0].set_title('implicit')
ax[1].set_title('explicit (i.e. vs independant variable time)')
plt.show()

您会注意到N1 和N2 都演变为某个最终值，但两个值不同。对于给定的方程，隐式图中的曲线不相交。

【讨论】：