如何使用scipy绘制两物种互利共生模型的相平面/相图

核心实现步骤 #

• 生成N1、N2的二维网格作为相平面坐标基底，覆盖种群的合理取值范围
• 计算每个网格点对应的dN1/dt、dN2/dt导数值，用于绘制向量场
• 求解dN1/dt=0、dN2/dt=0的零倾线，标记模型的平衡点位置
• 代入多组初始值求解ODE，绘制不同初始条件下的种群演化轨迹

完整可运行代码 #

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

# 原有模型参数
r1 = 1.0
r2 = 0.5
e1 = 1
e2 = 0.75
a12 = 0.25
a21 = 0.25

# 模型方程（保留原有定义）
def mutualism(X,t):
    N1, N2 = X
    dX = np.zeros(2)
    dX[0] = N1 * (r1 - (e1 * N1) + (a12 * N2))
    dX[1] = N2 * (r2 - (e2 * N2) + (a21 * N1))
    return dX

# --------------------------
# 原有时间序列绘制（可保留）
# --------------------------
X0 = [1, 1]
t = np.linspace(0, 100, 1000)
X = odeint(mutualism, X0, t)
N1, N2 = X[:,0], X[:,1]

plt.figure(figsize=(10,4))
plt.subplot(121)
plt.plot(t, N1, 'r-', label='物种1')
plt.plot(t, N2, 'b-', label='物种2')
plt.grid()
plt.legend(loc='best')
plt.xlabel('时间')
plt.ylabel('种群数量')
plt.title('互利共生模型时间序列')

# --------------------------
# 新增相平面绘制部分
# --------------------------
plt.subplot(122)
# 1. 生成相平面网格
n1_grid = np.linspace(0, 3, 20)
n2_grid = np.linspace(0, 3, 20)
N1_grid, N2_grid = np.meshgrid(n1_grid, n2_grid)

# 2. 计算网格点的导数
dN1_dt, dN2_dt = np.zeros(N1_grid.shape), np.zeros(N2_grid.shape)
for i in range(N1_grid.shape[0]):
    for j in range(N1_grid.shape[1]):
        dX = mutualism([N1_grid[i,j], N2_grid[i,j]], 0)
        dN1_dt[i,j] = dX[0]
        dN2_dt[i,j] = dX[1]

# 归一化导数，让箭头长度一致只显示方向
norm = np.sqrt(dN1_dt**2 + dN2_dt**2)
dN1_dt_norm = dN1_dt / norm
dN2_dt_norm = dN2_dt / norm

# 绘制向量场
plt.quiver(N1_grid, N2_grid, dN1_dt_norm, dN2_dt_norm, color='gray', alpha=0.6)

# 3. 绘制零倾线
n1_line = np.linspace(0, 3, 100)
# dN1/dt=0的零倾线（N1≠0时）
n2_null1 = (e1 * n1_line - r1)/a12
plt.plot(n1_line, n2_null1, 'r--', label='物种1零倾线')
# dN2/dt=0的零倾线（N2≠0时）
n2_null2 = (r2 + a21 * n1_line)/e2
plt.plot(n1_line, n2_null2, 'b--', label='物种2零倾线')

# 4. 绘制多组初值的轨迹
init_list = [[1,1], [0.2, 0.5], [2, 0.5], [0.5, 2]]
t_phase = np.linspace(0, 100, 1000)
for init in init_list:
    X_phase = odeint(mutualism, init, t_phase)
    plt.plot(X_phase[:,0], X_phase[:,1], linewidth=1.5)
    plt.scatter(init[0], init[1], c='black', s=20) # 标记初始点

# 相平面样式设置
plt.xlim(0, 3)
plt.ylim(0, 3)
plt.xlabel('物种1种群数量')
plt.ylabel('物种2种群数量')
plt.legend()
plt.grid(alpha=0.3)
plt.title('互利共生模型相平面')

plt.tight_layout()
plt.show()

小提示：可以根据你的参数调整网格范围、向量密度、初值列表，来获得更符合你需求的可视化效果。如果不需要保留时间序列图，删掉对应子图的代码即可。

内容的提问来源于stack exchange，提问作者thepajama