计算科学 - 在 Python 中求解和绘制互惠模型 - 吾爱随笔录

在 Python 中求解和绘制互惠模型

计算科学 Python 颂数值建模

2021-12-12 00:54:41

我是编程初学者。我需要在 python 中编写两个物种的共生模型，该模型将使用以下方程求解和绘制图形：

\frac{d N_{1}}{d t} = N_{1} (r_{1} - e_{1} N_{1} + α_{12} N_{2})

$\frac{dN_1}{dt} = N_1(r_1 - e_1N_1 + \alpha _{12} N_2)$ 和

\frac{d N_{2}}{d t} = N_{2} (r_{1} - e_{1} N_{2} + α_{21} N_{1})

$\frac{dN_2}{dt} = N_2(r_1 - e_1N_2 + \alpha _{21} N_1)$

$N_1$ 是物种 1 的种群。的值是物种 1 对物种 2 的影响（反之亦然），并且会在 0-2 之间变化。同时，、、和具有固定值。 $\alpha_{12}$ $r_1$ $r_2$ $e_1$ $e_2$

我尝试编写与 Lotka-Volterra 模型类似的代码，但没有成功。

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

# parameters
r1 = 1.0  #fixed
r2 = 0.5  #fixed
e1 = 1.0  #fixed
e2 = 0.75 #fixed

alpha12 = 1.5 #vary from 0-2
alpha21 = 1 #vary from 0-2

N1_0 = 1
N2_0 = 1

# store initial values in an array
X0 = [N1_0,N2_0]

# The two equations are contained in dX
def mutualism(X,t):
  N1, N2 = X
  dX = np.zeros(2)
  dX[0] = N1 * (r1 - (e1 * N1) + (alpha12 * N2)) # equation for dN1/dt
  dX[1] = N2 * (r2 - (e2 * N2) + (alpha21 * N1)) # equation for dN2/dt
  return dX

# set time length
t = np.linspace(0,100,300)

X = odeint(mutualism,X0,t)

N1 = X[:,0]; N2 = X[:,1]

plt.plot(t,N1, color='blue', lw=3)
plt.plot(t,N2, color='red', lw=3)
plt.show()

谁能告诉我我做错了什么以及如何改进代码？

1个回答

查看 scipy 上的文档，他们建议使用scipy.integrate.solve_ivp而不是scipy.integrate.odeint 在这里看到的。

因此我改变了功能并重新安排了一点。

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

# parameters
r1 = 1.0  #fixed
r2 = 0.5  #fixed
e1 = 1.0  #fixed
e2 = 0.75 #fixed

alpha12 = 1.5 #vary from 0-2
alpha21 = 1 #vary from 0-2

N1_0 = 1
N2_0 = 1

# store initial values in an array
X0 = [N1_0,N2_0]

# The two equations are contained in dX
def mutualism(t,X):
  N1, N2 = X
  dX = np.zeros(2)
  dX[0] = N1 * (r1 - (e1 * N1) + (alpha12 * N2)) # equation for dN1/dt
  dX[1] = N2 * (r2 - (e2 * N2) + (alpha21 * N1)) # equation for dN2/dt
  return dX

# set time length
t = np.linspace(0,1.6,100)

X = solve_ivp(mutualism,t[[0,-1]],X0,t_eval=t)

N1 = X.y[0]
N2 = X.y[1]
t_fin = X.t

plt.plot(t_fin,N1, color='blue', lw=3)
plt.plot(t_fin,N2, color='red', lw=3)
plt.show()

这给了我以下情节。（不确定这是否是您想要的，因为我并不真正了解共生模型。

如评论中所述，我们在 t~=1.5 处获得了一个趋于无穷大的奇点。因此，代码仅略高于该值。

其它你可能感兴趣的问题

上一篇打算尝试将我的一些 scipy/numpy 计算转移到新的 GPU，如何避免令人失望的结果？下一篇一种用于查找具有给定已知特征值的特征向量数量的方法