常微分方程（十五）种群动力学

为什么猞猁和雪兔的数量会呈现出惊人的周期性波动？为什么引入一种新物种有时会导致生态灾难？为什么有些物种能共存而另一些必然走向竞争排斥？这些问题的答案藏在微分方程中。本章我们将探索种群动力学的数学理论，从简单的单物种增长到复杂的多物种相互作用，看看数学如何揭示生态系统的深层规律。

单物种种群增长

Malthus 指数增长模型

最简单的种群模型假设增长率恒定：

其中是种群数量，是内禀增长率（出生率减去死亡率）。

解： $$

N(t) = N_0 e^{rt} $$ - 若$ r > 0 $：指数增长若$ r < 0 $：指数衰减若$ r = 0$：保持不变

问题：指数增长在长期内是不可能的——资源有限！

Logistic 增长模型

1838 年，比利时数学家皮埃尔·韦尔于斯特提出了更现实的模型：

其中是环境容纳量（ carrying capacity），代表环境能支持的最大种群数量。

解析解： $$

N(t) = {K + N_0(e^{rt} - 1)} = $$

关键特性： - 当：近似指数增长 - 当：增长率最大 - 当：增长率趋于零 - 是稳定平衡点

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

def logistic(N, t, r, K):
    """Logistic 增长模型"""
    return r * N * (1 - N/K)

# 参数
r = 0.5  # 内禀增长率
K = 1000  # 环境容纳量

# 不同初始值
t = np.linspace(0, 30, 500)
N0_values = [10, 100, 500, 900, 1200, 1500]
colors = plt.cm.viridis(np.linspace(0, 1, len(N0_values)))

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 时间演化
ax1 = axes[0]
for N0, color in zip(N0_values, colors):
    sol = odeint(logistic, N0, t, args=(r, K))
    ax1.plot(t, sol, color=color, linewidth=2, label=f'N ₀ = {N0}')

ax1.axhline(y=K, color='r', linestyle='--', linewidth=2, label=f'K = {K}')
ax1.axhline(y=K/2, color='orange', linestyle=':', linewidth=1.5, alpha=0.7)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Population N', fontsize=12)
ax1.set_title('Logistic Growth', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 相图（ N vs dN/dt）
ax2 = axes[1]
N_range = np.linspace(0, 1.5*K, 200)
dNdt = r * N_range * (1 - N_range/K)
ax2.plot(N_range, dNdt, 'b-', linewidth=2)
ax2.fill_between(N_range, 0, dNdt, where=dNdt>0, alpha=0.3, color='green', label='Growing')
ax2.fill_between(N_range, dNdt, 0, where=dNdt<0, alpha=0.3, color='red', label='Declining')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.axvline(x=K, color='r', linestyle='--', linewidth=2)
ax2.axvline(x=K/2, color='orange', linestyle=':', linewidth=1.5)
ax2.plot([0, K], [0, 0], 'ko', markersize=8)
ax2.set_xlabel('Population N', fontsize=12)
ax2.set_ylabel('dN/dt', fontsize=12)
ax2.set_title('Phase Portrait', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.text(K/2, max(dNdt)*0.9, 'Max growth\nrate', fontsize=10, ha='center')

plt.tight_layout()
plt.savefig('logistic_growth.png', dpi=150, bbox_inches='tight')
plt.show()

Allee 效应

在某些物种中，过低的种群密度反而导致增长率下降（难以找到配偶、无法抵御捕食者等）。这称为Allee 效应。

强 Allee 效应模型：

其中是 Allee 阈值。

若：种群衰减走向灭绝

若

：种群增长趋向

def allee_effect(N, t, r, K, A):
    """带 Allee 效应的模型"""
    if N <= 0:
        return 0
    return r * N * (1 - N/K) * (N/A - 1)

# 参数
r = 0.1
K = 1000
A = 100  # Allee 阈值

t = np.linspace(0, 100, 1000)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 时间演化
ax1 = axes[0]
N0_values = [50, 80, 100, 120, 200, 500, 1200]
colors = plt.cm.coolwarm(np.linspace(0, 1, len(N0_values)))

for N0, color in zip(N0_values, colors):
    sol = odeint(allee_effect, N0, t, args=(r, K, A))
    ax1.plot(t, sol, color=color, linewidth=2, label=f'N ₀ = {N0}')

ax1.axhline(y=K, color='g', linestyle='--', linewidth=2, alpha=0.7, label=f'K = {K}')
ax1.axhline(y=A, color='r', linestyle='--', linewidth=2, alpha=0.7, label=f'A = {A}')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Population N', fontsize=12)
ax1.set_title('Population Dynamics with Allee Effect', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 相图
ax2 = axes[1]
N_range = np.linspace(1, 1.3*K, 200)
dNdt = r * N_range * (1 - N_range/K) * (N_range/A - 1)
ax2.plot(N_range, dNdt, 'b-', linewidth=2)
ax2.fill_between(N_range, 0, dNdt, where=dNdt>0, alpha=0.3, color='green')
ax2.fill_between(N_range, dNdt, 0, where=dNdt<0, alpha=0.3, color='red')
ax2.axhline(y=0, color='k', linewidth=0.5)
ax2.plot([0, A, K], [0, 0, 0], 'o', markersize=10)
ax2.plot(0, 0, 'go', markersize=10, label='Stable')
ax2.plot(A, 0, 'ro', markersize=10, label='Unstable')
ax2.plot(K, 0, 'go', markersize=10)
ax2.set_xlabel('Population N', fontsize=12)
ax2.set_ylabel('dN/dt', fontsize=12)
ax2.set_title('Allee Effect Phase Portrait', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('allee_effect.png', dpi=150, bbox_inches='tight')
plt.show()

Allee 效应解释了为什么某些濒危物种一旦数量过低就很难恢复。

捕食者-猎物模型

Lotka-Volterra 模型

1925-1926 年，洛特卡和沃尔泰拉独立提出了捕食者-猎物相互作用的模型：

其中： - ：猎物数量 - ：捕食者数量 - ：猎物内禀增长率 - ：捕食率 - ：捕食者死亡率 - ：能量转化效率

平衡点分析：

平凡平衡点：— 双方灭绝
共存平衡点：在共存平衡点处的 Jacobi 矩阵：

J =

特征值：（纯虚数）

结论：共存平衡点是中心，周围存在无穷多条闭合轨道。

周期振荡

Lotka-Volterra 模型的一个显著特征是种群数量的周期性振荡：

def lotka_volterra(y, t, alpha, beta, gamma, delta):
    """Lotka-Volterra 捕食者-猎物模型"""
    x, predator = y
    dxdt = alpha * x - beta * x * predator
    dydt = delta * x * predator - gamma * predator
    return [dxdt, dydt]

# 参数
alpha = 1.0   # 猎物增长率
beta = 0.1    # 捕食率
gamma = 0.5   # 捕食者死亡率
delta = 0.02  # 能量转化效率

# 平衡点
x_eq = gamma / delta
y_eq = alpha / beta
print(f"共存平衡点: ({x_eq:.1f}, {y_eq:.1f})")

t = np.linspace(0, 100, 2000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 不同初始值的轨迹
ax1 = axes[0, 0]
initial_conditions = [(40, 9), (20, 5), (30, 15), (50, 8)]
colors = ['blue', 'green', 'red', 'purple']

for (x0, y0), color in zip(initial_conditions, colors):
    sol = odeint(lotka_volterra, [x0, y0], t, args=(alpha, beta, gamma, delta))
    ax1.plot(sol[:, 0], sol[:, 1], color=color, linewidth=1.5, alpha=0.7)
    ax1.plot(x0, y0, 'o', color=color, markersize=8)

ax1.plot(x_eq, y_eq, 'k*', markersize=15, label='Equilibrium')
ax1.set_xlabel('Prey x', fontsize=12)
ax1.set_ylabel('Predator y', fontsize=12)
ax1.set_title('Phase Portrait: Closed Orbits', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 时间演化
ax2 = axes[0, 1]
x0, y0 = 40, 9
sol = odeint(lotka_volterra, [x0, y0], t, args=(alpha, beta, gamma, delta))
ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='Prey')
ax2.plot(t, sol[:, 1]*5, 'r-', linewidth=2, label='Predator (× 5)')
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Population Oscillations', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# 猞猁-雪兔数据（模拟历史数据）
ax3 = axes[1, 0]
# 这里用模型模拟，实际数据来自 Hudson Bay Company 记录
t_data = np.linspace(1845, 1935, 91)
t_model = np.linspace(0, 90, 910)
sol_hare = odeint(lotka_volterra, [50, 8], t_model, args=(0.5, 0.02, 0.3, 0.01))

# 添加噪声使其更像真实数据
np.random.seed(42)
hare_data = sol_hare[::10, 0] * (1 + 0.1*np.random.randn(91))
lynx_data = sol_hare[::10, 1] * 5 * (1 + 0.1*np.random.randn(91))

ax3.plot(t_data, hare_data, 'b-', linewidth=2, label='Snowshoe Hare')
ax3.plot(t_data, lynx_data, 'r-', linewidth=2, label='Canada Lynx')
ax3.set_xlabel('Year', fontsize=12)
ax3.set_ylabel('Population (thousands)', fontsize=12)
ax3.set_title('Classic Example: Hare-Lynx Cycles', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# 守恒量
ax4 = axes[1, 1]
# Lotka-Volterra 的守恒量： V = δ x - γ ln(x) + β y - α ln(y)
def conserved_quantity(x, y):
    return delta*x - gamma*np.log(x) + beta*y - alpha*np.log(y)

x_range = np.linspace(5, 80, 100)
y_range = np.linspace(2, 25, 100)
X, Y = np.meshgrid(x_range, y_range)
V = delta*X - gamma*np.log(X) + beta*Y - alpha*np.log(Y)

contour = ax4.contour(X, Y, V, levels=20, cmap='viridis')
ax4.clabel(contour, inline=True, fontsize=8)
ax4.plot(x_eq, y_eq, 'r*', markersize=15)
ax4.set_xlabel('Prey x', fontsize=12)
ax4.set_ylabel('Predator y', fontsize=12)
ax4.set_title('Conserved Quantity (Level Sets = Orbits)', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lotka_volterra.png', dpi=150, bbox_inches='tight')
plt.show()

结构不稳定性

经典 Lotka-Volterra 模型的一个问题是结构不稳定——小的扰动（如加入 Logistic 限制）会根本改变系统行为。

带 Logistic 限制的模型：

现在系统有稳定的极限环或稳定的焦点，取决于参数。

def lotka_volterra_logistic(y, t, alpha, beta, gamma, delta, K):
    """带 Logistic 限制的 Lotka-Volterra 模型"""
    x, pred = y
    dxdt = alpha * x * (1 - x/K) - beta * x * pred
    dydt = delta * x * pred - gamma * pred
    return [dxdt, dydt]

# 参数
alpha = 1.0
beta = 0.1
gamma = 0.5
delta = 0.02
K = 100  # 猎物的环境容纳量

t = np.linspace(0, 200, 2000)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 相图
ax1 = axes[0]
initial_conditions = [(10, 5), (80, 5), (50, 20), (30, 2)]
colors = ['blue', 'green', 'red', 'purple']

for (x0, y0), color in zip(initial_conditions, colors):
    sol = odeint(lotka_volterra_logistic, [x0, y0], t, args=(alpha, beta, gamma, delta, K))
    ax1.plot(sol[:, 0], sol[:, 1], color=color, linewidth=1, alpha=0.8)
    ax1.plot(x0, y0, 'o', color=color, markersize=8)

# 平衡点
x_eq = gamma / delta
y_eq = alpha/beta * (1 - x_eq/K)
ax1.plot(x_eq, y_eq, 'k*', markersize=15, label=f'Equilibrium ({x_eq:.1f}, {y_eq:.1f})')
ax1.set_xlabel('Prey x', fontsize=12)
ax1.set_ylabel('Predator y', fontsize=12)
ax1.set_title('Lotka-Volterra with Logistic Prey', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# 时间演化
ax2 = axes[1]
sol = odeint(lotka_volterra_logistic, [10, 5], t, args=(alpha, beta, gamma, delta, K))
ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='Prey')
ax2.plot(t, sol[:, 1]*5, 'r-', linewidth=2, label='Predator (× 5)')
ax2.axhline(y=x_eq, color='b', linestyle='--', alpha=0.5)
ax2.axhline(y=y_eq*5, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Damped Oscillations to Equilibrium', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lotka_volterra_logistic.png', dpi=150, bbox_inches='tight')
plt.show()

功能响应

在经典模型中，捕食率与猎物数量成正比。但实际上，捕食者有处理时间和饱和效应。

Holling II 型功能响应：

其中是处理时间，是最大捕食率。

Holling III 型功能响应：

这种 S 形响应在低猎物密度时提供了"庇护所"效应。

def rosenzweig_macarthur(state, t, r, K, a, c, m, e):
    """Rosenzweig-MacArthur 模型（ Holling II 型功能响应）"""
    x, y = state
    if x <= 0:
        x = 0
    functional_response = c * x / (a + x)
    dxdt = r * x * (1 - x/K) - functional_response * y
    dydt = e * functional_response * y - m * y
    return [dxdt, dydt]

# 参数
r = 1.0    # 猎物增长率
K = 1.5    # 环境容纳量
a = 0.5    # 半饱和常数
c = 1.0    # 最大捕食率
m = 0.4    # 捕食者死亡率
e = 0.6    # 转化效率

t = np.linspace(0, 200, 2000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 功能响应曲线
ax1 = axes[0, 0]
x_range = np.linspace(0, 3, 100)
# Holling I
ax1.plot(x_range, 0.5*x_range, 'b-', linewidth=2, label='Type I (Linear)')
# Holling II
ax1.plot(x_range, c*x_range/(a + x_range), 'r-', linewidth=2, label='Type II')
# Holling III
ax1.plot(x_range, x_range**2/(0.5**2 + x_range**2), 'g-', linewidth=2, label='Type III')
ax1.set_xlabel('Prey density', fontsize=12)
ax1.set_ylabel('Predation rate', fontsize=12)
ax1.set_title('Holling Functional Responses', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 3)
ax1.set_ylim(0, 1.5)

# Rosenzweig-MacArthur 相图
ax2 = axes[0, 1]
initial_conditions = [(0.5, 0.3), (1.0, 0.1), (0.2, 0.5), (1.5, 0.4)]
colors = ['blue', 'green', 'red', 'purple']

for (x0, y0), color in zip(initial_conditions, colors):
    sol = odeint(rosenzweig_macarthur, [x0, y0], t, args=(r, K, a, c, m, e))
    ax2.plot(sol[:, 0], sol[:, 1], color=color, linewidth=1.5, alpha=0.8)
    ax2.plot(x0, y0, 'o', color=color, markersize=8)

ax2.set_xlabel('Prey x', fontsize=12)
ax2.set_ylabel('Predator y', fontsize=12)
ax2.set_title('Rosenzweig-MacArthur Model', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# 时间演化
ax3 = axes[1, 0]
sol = odeint(rosenzweig_macarthur, [0.5, 0.3], t, args=(r, K, a, c, m, e))
ax3.plot(t, sol[:, 0], 'b-', linewidth=2, label='Prey')
ax3.plot(t, sol[:, 1], 'r-', linewidth=2, label='Predator')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Population', fontsize=12)
ax3.set_title('Population Dynamics', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# 不同 K 值的影响（丰富悖论）
ax4 = axes[1, 1]
K_values = [0.5, 1.0, 1.5, 2.5]
colors_K = plt.cm.plasma(np.linspace(0.2, 0.8, len(K_values)))

for K_val, color in zip(K_values, colors_K):
    sol = odeint(rosenzweig_macarthur, [0.5, 0.3], t, args=(r, K_val, a, c, m, e))
    ax4.plot(sol[-1000:, 0], sol[-1000:, 1], color=color, linewidth=1, 
             alpha=0.7, label=f'K = {K_val}')

ax4.set_xlabel('Prey x', fontsize=12)
ax4.set_ylabel('Predator y', fontsize=12)
ax4.set_title('Paradox of Enrichment', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('functional_response.png', dpi=150, bbox_inches='tight')
plt.show()

丰富悖论

一个反直觉的现象：增加环境容纳量可能使系统不稳定！

当增大时： 1. 小：稳定平衡点 2. 中等： Hopf 分岔，出现稳定极限环 3. 大：极限环振幅增大，可能导致灭绝

这被称为丰富悖论（ Paradox of Enrichment）。

竞争模型

Lotka-Volterra 竞争方程

两个物种竞争相同资源：

其中： - ：物种 2 对物种 1 的竞争系数 - ：物种 1 对物种 2 的竞争系数

四种结果

根据参数关系，系统可能有四种结果：

物种 1 获胜：且 $K_1/{12} > K_2 $物种获胜：$ K_2 > K_1/{12} $且$ K_2/{21} > K_1 $稳定共存：$ K_1 < K_2/{21} $且$ K_2 < K_1/{12} $双稳态：$ K_1 > K_2/{21}$ 且 竞争排斥原理：两个生态位完全重叠的物种不能稳定共存。

def competition_model(state, t, r1, r2, K1, K2, alpha12, alpha21):
    """Lotka-Volterra 竞争模型"""
    N1, N2 = state
    dN1dt = r1 * N1 * (1 - (N1 + alpha12*N2)/K1)
    dN2dt = r2 * N2 * (1 - (N2 + alpha21*N1)/K2)
    return [dN1dt, dN2dt]

# 四种情况
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

cases = [
    # (r1, r2, K1, K2, alpha12, alpha21, title)
    (1.0, 1.0, 100, 80, 0.5, 0.5, 'Species 1 Wins'),
    (1.0, 1.0, 80, 100, 0.5, 0.5, 'Species 2 Wins'),
    (1.0, 1.0, 100, 100, 0.5, 0.5, 'Stable Coexistence'),
    (1.0, 1.0, 100, 100, 1.5, 1.5, 'Bistability (Winner depends on initial)'),
]

t = np.linspace(0, 50, 500)

for ax, (r1, r2, K1, K2, alpha12, alpha21, title) in zip(axes.flat, cases):
    # 画零增长等斜线
    N1_range = np.linspace(0, 150, 100)
    
    # N1 零增长线： N1 + alpha12*N2 = K1 => N2 = (K1 - N1)/alpha12
    N2_null1 = (K1 - N1_range) / alpha12
    N2_null1 = np.clip(N2_null1, 0, 150)
    
    # N2 零增长线： N2 + alpha21*N1 = K2 => N2 = K2 - alpha21*N1
    N2_null2 = K2 - alpha21 * N1_range
    N2_null2 = np.clip(N2_null2, 0, 150)
    
    ax.plot(N1_range, N2_null1, 'b-', linewidth=2, label='N ₁ nullcline')
    ax.plot(N1_range, N2_null2, 'r-', linewidth=2, label='N ₂ nullcline')
    
    # 画轨迹
    initial_conditions = [(10, 10), (10, 100), (100, 10), (100, 100), (50, 50)]
    for N1_0, N2_0 in initial_conditions:
        sol = odeint(competition_model, [N1_0, N2_0], t, 
                    args=(r1, r2, K1, K2, alpha12, alpha21))
        ax.plot(sol[:, 0], sol[:, 1], 'g-', linewidth=0.8, alpha=0.7)
        ax.plot(N1_0, N2_0, 'ko', markersize=5)
        ax.plot(sol[-1, 0], sol[-1, 1], 'r*', markersize=10)
    
    ax.set_xlabel('Species 1 (N ₁)', fontsize=12)
    ax.set_ylabel('Species 2 (N ₂)', fontsize=12)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.legend(fontsize=9, loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 150)
    ax.set_ylim(0, 150)

plt.tight_layout()
plt.savefig('competition_model.png', dpi=150, bbox_inches='tight')
plt.show()

共存条件的生态学解释

稳定共存的条件意味着：

种内竞争强于种间竞争

直观地说，如果一个物种对自己的影响大于对另一物种的影响，那么当它的数量增加时，它会更多地限制自己而不是对手，这为对手留下了生存空间。

这解释了为什么物种分化（占据不同生态位）有利于共存。

互利共生模型

基本模型

两个物种相互受益：

其中表示互利程度。

不稳定性问题

如果互利太强（），可能导致种群爆炸——这是生态学中的一个重要警告。

def mutualism_model(state, t, r1, r2, K1, K2, beta12, beta21):
    """互利共生模型"""
    N1, N2 = state
    K1_eff = K1 + beta12 * N2
    K2_eff = K2 + beta21 * N1
    dN1dt = r1 * N1 * (1 - N1/K1_eff)
    dN2dt = r2 * N2 * (1 - N2/K2_eff)
    return [dN1dt, dN2dt]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 50, 500)

# 弱互利（稳定）
ax1 = axes[0]
r1, r2 = 1.0, 1.0
K1, K2 = 100, 100
beta12, beta21 = 0.3, 0.3

initial_conditions = [(10, 10), (50, 20), (20, 50)]
for N1_0, N2_0 in initial_conditions:
    sol = odeint(mutualism_model, [N1_0, N2_0], t, 
                args=(r1, r2, K1, K2, beta12, beta21))
    ax1.plot(sol[:, 0], sol[:, 1], linewidth=1.5)
    ax1.plot(N1_0, N2_0, 'ko', markersize=6)

ax1.set_xlabel('Species 1', fontsize=12)
ax1.set_ylabel('Species 2', fontsize=12)
ax1.set_title('Weak Mutualism (Stable)', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 强互利（不稳定）
ax2 = axes[1]
beta12, beta21 = 1.5, 1.5

for N1_0, N2_0 in initial_conditions:
    try:
        sol = odeint(mutualism_model, [N1_0, N2_0], t, 
                    args=(r1, r2, K1, K2, beta12, beta21))
        sol = np.clip(sol, 0, 1000)  # 限制范围以便可视化
        ax2.plot(sol[:, 0], sol[:, 1], linewidth=1.5)
        ax2.plot(N1_0, N2_0, 'ko', markersize=6)
    except:
        pass

ax2.set_xlabel('Species 1', fontsize=12)
ax2.set_ylabel('Species 2', fontsize=12)
ax2.set_title('Strong Mutualism (Runaway Growth)', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('mutualism.png', dpi=150, bbox_inches='tight')
plt.show()

食物网与多物种动力学

三物种食物链

考虑植物-草食动物-肉食动物的食物链：

def food_chain(state, t, r, K, a1, b1, e1, a2, b2, e2, m1, m2):
    """三物种食物链模型"""
    P, H, C = state
    
    # 功能响应
    f1 = a1 * P / (1 + b1 * P)  # 草食动物对植物
    f2 = a2 * H / (1 + b2 * H)  # 肉食动物对草食动物
    
    dPdt = r * P * (1 - P/K) - f1 * H
    dHdt = e1 * f1 * H - f2 * C - m1 * H
    dCdt = e2 * f2 * C - m2 * C
    
    return [dPdt, dHdt, dCdt]

# 参数
r = 1.0      # 植物增长率
K = 1.0      # 植物容纳量
a1 = 5.0     # 草食动物捕食率
b1 = 3.0     # 半饱和
e1 = 0.5     # 转化效率
a2 = 0.1     # 肉食动物捕食率
b2 = 2.0     # 半饱和
e2 = 0.5     # 转化效率
m1 = 0.4     # 草食动物死亡率
m2 = 0.01    # 肉食动物死亡率

t = np.linspace(0, 500, 5000)
initial = [0.5, 0.5, 0.5]

sol = odeint(food_chain, initial, t, args=(r, K, a1, b1, e1, a2, b2, e2, m1, m2))

fig = plt.figure(figsize=(14, 10))

# 3D 相图
ax1 = fig.add_subplot(221, projection='3d')
ax1.plot(sol[:, 0], sol[:, 1], sol[:, 2], 'b-', linewidth=0.5, alpha=0.8)
ax1.set_xlabel('Plant P', fontsize=10)
ax1.set_ylabel('Herbivore H', fontsize=10)
ax1.set_zlabel('Carnivore C', fontsize=10)
ax1.set_title('3D Phase Space', fontsize=12, fontweight='bold')

# 时间演化
ax2 = fig.add_subplot(222)
ax2.plot(t, sol[:, 0], 'g-', linewidth=1.5, label='Plant P')
ax2.plot(t, sol[:, 1], 'b-', linewidth=1.5, label='Herbivore H')
ax2.plot(t, sol[:, 2], 'r-', linewidth=1.5, label='Carnivore C')
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Time Series', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# P-H 投影
ax3 = fig.add_subplot(223)
ax3.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=0.5, alpha=0.8)
ax3.set_xlabel('Plant P', fontsize=12)
ax3.set_ylabel('Herbivore H', fontsize=12)
ax3.set_title('P-H Phase Plane', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)

# H-C 投影
ax4 = fig.add_subplot(224)
ax4.plot(sol[:, 1], sol[:, 2], 'r-', linewidth=0.5, alpha=0.8)
ax4.set_xlabel('Herbivore H', fontsize=12)
ax4.set_ylabel('Carnivore C', fontsize=12)
ax4.set_title('H-C Phase Plane', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('food_chain.png', dpi=150, bbox_inches='tight')
plt.show()

混沌与食物网

三物种及以上的食物网可能表现出混沌行为——种群数量看似随机波动，但实际上由确定性方程控制。

这对生态管理有重要影响：即使完全了解系统的数学结构，长期预测也可能不可行。

空间种群动力学

反应-扩散方程

当考虑空间分布时，种群模型变成偏微分方程：

这是著名的Fisher-KPP 方程。

行波解

Fisher-KPP 方程存在行波解——种群以恒定速度在空间中扩散：

N(x, t) = N(x - ct) $$

最小波速。

图灵斑图

在捕食者-猎物系统中加入扩散可能导致图灵不稳定性——空间均匀的平衡态变得不稳定，形成空间斑图。

# 简化的 1D 空间模型演示
def fisher_kpp_1d(N, t, D, r, K, dx, L):
    """1D Fisher-KPP 方程的有限差分离散"""
    n = len(N)
    dNdt = np.zeros(n)
    
    # 扩散项（中心差分）
    for i in range(1, n-1):
        diffusion = D * (N[i+1] - 2*N[i] + N[i-1]) / dx**2
        reaction = r * N[i] * (1 - N[i]/K)
        dNdt[i] = diffusion + reaction
    
    # 边界条件（ Neumann）
    dNdt[0] = dNdt[1]
    dNdt[-1] = dNdt[-2]
    
    return dNdt

# 参数
D = 0.1      # 扩散系数
r = 0.5      # 增长率
K = 1.0      # 容纳量
L = 50       # 空间范围
nx = 200     # 空间网格数
dx = L / nx
x = np.linspace(0, L, nx)

# 初始条件：左边一小块种群
N0 = np.zeros(nx)
N0[:10] = K

t = np.linspace(0, 100, 200)
sol = odeint(fisher_kpp_1d, N0, t, args=(D, r, K, dx, L))

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 时空图
ax1 = axes[0]
im = ax1.imshow(sol.T, aspect='auto', origin='lower', 
                extent=[0, 100, 0, L], cmap='viridis')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Space x', fontsize=12)
ax1.set_title('Fisher-KPP Wave Propagation', fontsize=14, fontweight='bold')
plt.colorbar(im, ax=ax1, label='Population N')

# 不同时刻的空间分布
ax2 = axes[1]
times_to_plot = [0, 20, 40, 60, 80, 100]
colors = plt.cm.plasma(np.linspace(0, 0.8, len(times_to_plot)))

for time_val, color in zip(times_to_plot, colors):
    idx = int(time_val / 100 * 199)
    ax2.plot(x, sol[idx, :], color=color, linewidth=2, label=f't = {time_val}')

ax2.set_xlabel('Space x', fontsize=12)
ax2.set_ylabel('Population N', fontsize=12)
ax2.set_title('Population Profiles at Different Times', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('fisher_kpp.png', dpi=150, bbox_inches='tight')
plt.show()

# 估计波速
# 找到 N = K/2 的位置随时间变化
front_positions = []
for i, t_val in enumerate(t):
    try:
        idx = np.where(sol[i, :] > K/2)[0][-1]
        front_positions.append(x[idx])
    except:
        front_positions.append(np.nan)

front_positions = np.array(front_positions)
valid = ~np.isnan(front_positions)

# 线性拟合估计波速
from scipy.stats import linregress
slope, intercept, r_value, p_value, std_err = linregress(t[valid][50:], front_positions[valid][50:])
print(f"估计波速: {slope:.3f}")
print(f"理论最小波速: {2*np.sqrt(D*r):.3f}")

总结

本章我们探索了种群动力学的丰富数学世界：

单物种模型： Malthus 、 Logistic 、 Allee 效应
捕食者-猎物模型： Lotka-Volterra 、功能响应、丰富悖论
竞争模型：竞争排斥、共存条件
互利共生：稳定性条件
多物种系统：食物链、混沌
空间动力学：反应-扩散、行波、图灵斑图

这些模型虽然简化了自然界的复杂性，但揭示了生态系统的基本原理，为保护生物学、渔业管理、入侵物种控制等提供了理论基础。

练习题

基础题

求解 Logistic 方程的解析解，并验证初值条件和渐近行为。
对于 Lotka-Volterra 竞争模型，画出四种可能结果的相图，并标出所有平衡点及其稳定性。
证明经典 Lotka-Volterra 捕食者-猎物模型的守恒量： $$

V = x - x + y - y$$4. 对于带 Allee 效应的模型，找出所有平衡点并分析其稳定性。

如果一种入侵物种的$ r = 0.3 $年，扩散系数$ D = 10$ km ²/年，计算其理论扩散速度。

进阶题

收割模型：考虑 Logistic 种群加上恒定收割：

分析对平衡点的影响
找出最大可持续收割量

年龄结构：将种群分为幼年和成年两个阶段：

分析系统的长期行为。

证明 Rosenzweig-MacArthur 模型在某个临界值处发生 Hopf 分岔。
延迟模型：考虑带时滞的 Logistic 方程：

讨论时滞对稳定性的影响。

对于三物种竞争模型，什么条件下可能出现"石头-剪刀-布"式的周期共存？

编程题

编写程序模拟经典 Lotka-Volterra 模型，并用真实的猞猁-雪兔数据拟合参数。
实现一个基于格子的空间捕食者-猎物模型，可视化空间斑图的形成。
编写代码计算竞争模型的相图，自动判断结果类型。
模拟一个具有随机环境波动的种群模型：

其中是白噪声。

实现一个食物网模型，包含至少 5 个物种，探索其动力学行为。

思考题

为什么真实的捕食者-猎物系统通常是稳定的，而经典 Lotka-Volterra 模型是结构不稳定的？
讨论"丰富悖论"的生态学意义。在实际中，这种悖论为什么不经常观察到？
气候变化如何影响种群动力学模型的参数？这对保护生物学有什么启示？
比较确定性模型和随机模型在小种群保护中的应用。
人工智能和机器学习能否取代传统的种群动力学模型？讨论各自的优缺点。

参考资料

Murray, J. D. (2002). Mathematical Biology I: An Introduction. Springer.
Kot, M. (2001). Elements of Mathematical Ecology. Cambridge University Press.
Edelstein-Keshet, L. (2005). Mathematical Models in Biology. SIAM.
Hastings, A. (1997). Population Biology: Concepts and Models. Springer.
May, R. M. (1974). "Biological Populations with Nonoverlapping Generations: Stable Points, Stable Cycles, and Chaos." Science, 186(4164), 645-647.
Rosenzweig, M. L. (1971). "Paradox of Enrichment: Destabilization of Exploitation Ecosystems in Ecological Time." Science, 171(3969), 385-387.
Tilman, D. (1982). Resource Competition and Community Structure. Princeton University Press.

本文是《常微分方程的世界》系列的第 15 章。