使用SciPy解决Lotka-Volterra模型

Question

我有以下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中执行此操作？

我在线阅读本教程(幻灯片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[:])

但情节不正确.更新:我认为我使用了错误的方法.我正在阅读另一个在线教程.我会更多地解决这个问题.在此期间,如果有人知道如何解决它,请告诉我.

Answer 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演变为某些最终值,但两个值都不同.隐式图中的曲线与给定方程不相交.