用Python模拟兔子和羊的生存竞争：从Lotka-Volterra模型到代码实现
生态系统的动态平衡一直是科学家们研究的重点课题。想象一片广袤的草原，兔子和羊作为主要食草动物，它们之间的竞争关系直接影响着整个草场的生态平衡。这种看似简单的生物互动背后，隐藏着精妙的数学规律。本文将带你用Python代码重现这一生态过程，通过Lotka-Volterra模型揭示种群竞争的奥秘。
1. 理解Lotka-Volterra竞争模型
Lotka-Volterra模型是描述两个物种竞争关系的经典数学模型。它由美国数学家Alfred Lotka和意大利数学家Vito Volterra在20世纪20年代独立提出，最初用于解释渔业中掠食者与被掠食者的数量变化，后来被广泛应用于各种生态竞争场景。
模型的核心是两个微分方程：
def lotka_volterra(t, N, r1, r2, K1, K2, alpha12, alpha21): N1, N2 = N dN1_dt = r1 * N1 * (1 - (N1 + alpha12 * N2) / K1) dN2_dt = r2 * N2 * (1 - (N2 + alpha21 * N1) / K2) return [dN1_dt, dN2_dt]其中关键参数含义如下表：
| 参数 | 描述 | 典型取值范围 |
|---|---|---|
| r1, r2 | 兔子和羊的内禀增长率 | 0.1-1.0 |
| K1, K2 | 环境对兔子和羊的承载能力 | 100-1000 |
| alpha12 | 羊对兔子的竞争系数 | 0.1-2.0 |
| alpha21 | 兔子对羊的竞争系数 | 0.1-2.0 |
模型可以预测四种可能的竞争结果：
- 兔子获胜：当兔子能有效抵抗羊的竞争时
- 羊获胜：当羊在竞争中占据绝对优势时
- 稳定共存：当两者达到某种平衡状态时
- 结果不确定：取决于初始种群数量
2. Python环境准备与库安装
要实现这个模型，我们需要以下Python库：
pip install numpy scipy matplotlib ipywidgets各库的作用如下：
- NumPy：处理数值计算和数组操作
- SciPy：求解微分方程
- Matplotlib：可视化模拟结果
- IPyWidgets（可选）：创建交互式控件
建议使用Jupyter Notebook进行实验，可以实时观察代码运行结果。以下是基础环境检查代码：
import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp print("NumPy版本:", np.__version__) print("SciPy版本:", scipy.__version__) print("Matplotlib版本:", matplotlib.__version__)3. 实现竞争模型的核心代码
3.1 定义微分方程
我们需要将数学模型转化为Python函数。使用SciPy的solve_ivp函数可以方便地求解常微分方程：
def competition_model(t, N, r1, r2, K1, K2, alpha12, alpha21): """Lotka-Volterra竞争模型""" N1, N2 = N # 解包当前种群数量 dN1 = r1 * N1 * (1 - (N1 + alpha12 * N2)/K1) dN2 = r2 * N2 * (1 - (N2 + alpha21 * N1)/K2) return [dN1, dN2]3.2 设置模拟参数并求解
选择适当的参数值和初始条件：
# 模型参数 params = { 'r1': 0.5, # 兔子的增长率 'r2': 0.4, # 羊的增长率 'K1': 500, # 兔子的环境承载力 'K2': 600, # 羊的环境承载力 'alpha12': 0.7, # 羊对兔子的竞争影响 'alpha21': 0.5 # 兔子对羊的竞争影响 } # 初始种群数量 N0 = [100, 150] # 初始兔子数量, 初始羊数量 # 时间范围 (0到50个时间单位) t_span = (0, 50) t_eval = np.linspace(*t_span, 500) # 500个时间点求解微分方程：
solution = solve_ivp( competition_model, t_span, N0, args=tuple(params.values()), t_eval=t_eval, method='RK45' # Runge-Kutta方法 )3.3 可视化结果
创建种群数量随时间变化的曲线：
plt.figure(figsize=(10, 6)) plt.plot(solution.t, solution.y[0], label='兔子数量') plt.plot(solution.t, solution.y[1], label='羊数量') plt.xlabel('时间') plt.ylabel('种群数量') plt.title('兔子和羊的种群竞争动态') plt.legend() plt.grid(True) plt.show()还可以创建相位图，展示两个种群数量的关系：
plt.figure(figsize=(8, 8)) plt.plot(solution.y[0], solution.y[1]) plt.xlabel('兔子数量') plt.ylabel('羊数量') plt.title('种群竞争相位图') plt.grid(True) plt.show()4. 探索不同竞争场景
通过调整参数，我们可以模拟四种不同的竞争结果。下面创建交互式可视化，方便探索参数空间。
4.1 兔子获胜的场景
当兔子在竞争中占据绝对优势时：
params_rabbit_win = { 'r1': 0.5, 'r2': 0.4, 'K1': 500, 'K2': 300, 'alpha12': 0.5, 'alpha21': 1.2 } solution = solve_ivp(competition_model, t_span, N0, args=tuple(params_rabbit_win.values()), t_eval=t_eval) plt.figure(figsize=(10, 6)) plt.plot(solution.t, solution.y[0], label='兔子') plt.plot(solution.t, solution.y[1], label='羊') plt.title('兔子获胜的竞争场景') plt.legend(); plt.grid(True) plt.show()4.2 羊获胜的场景
当羊在竞争中占据优势时：
params_sheep_win = { 'r1': 0.4, 'r2': 0.5, 'K1': 300, 'K2': 500, 'alpha12': 1.2, 'alpha21': 0.5 } solution = solve_ivp(competition_model, t_span, N0, args=tuple(params_sheep_win.values()), t_eval=t_eval) plt.figure(figsize=(10, 6)) plt.plot(solution.t, solution.y[0], label='兔子') plt.plot(solution.t, solution.y[1], label='羊') plt.title('羊获胜的竞争场景') plt.legend(); plt.grid(True) plt.show()4.3 稳定共存的场景
当两个物种达到平衡时：
params_coexist = { 'r1': 0.5, 'r2': 0.5, 'K1': 500, 'K2': 500, 'alpha12': 0.5, 'alpha21': 0.5 } solution = solve_ivp(competition_model, t_span, N0, args=tuple(params_coexist.values()), t_eval=t_eval) plt.figure(figsize=(10, 6)) plt.plot(solution.t, solution.y[0], label='兔子') plt.plot(solution.t, solution.y[1], label='羊') plt.title('稳定共存的竞争场景') plt.legend(); plt.grid(True) plt.show()4.4 结果不确定的场景
当结果取决于初始条件时：
params_unstable = { 'r1': 0.5, 'r2': 0.5, 'K1': 500, 'K2': 500, 'alpha12': 0.8, 'alpha21': 0.8 } # 尝试不同的初始条件 N0_case1 = [100, 150] N0_case2 = [150, 100] solution1 = solve_ivp(competition_model, t_span, N0_case1, args=tuple(params_unstable.values()), t_eval=t_eval) solution2 = solve_ivp(competition_model, t_span, N0_case2, args=tuple(params_unstable.values()), t_eval=t_eval) plt.figure(figsize=(10, 6)) plt.plot(solution1.t, solution1.y[0], 'b-', label='兔子 (初始100)') plt.plot(solution1.t, solution1.y[1], 'g-', label='羊 (初始150)') plt.plot(solution2.t, solution2.y[0], 'b--', label='兔子 (初始150)') plt.plot(solution2.t, solution2.y[1], 'g--', label='羊 (初始100)') plt.title('初始条件影响竞争结果') plt.legend(); plt.grid(True) plt.show()5. 高级应用与扩展
5.1 参数敏感性分析
了解哪些参数对结果影响最大：
def analyze_sensitivity(base_params, param_name, values): """分析单个参数的敏感性""" results = [] for value in values: params = base_params.copy() params[param_name] = value sol = solve_ivp(competition_model, t_span, N0, args=tuple(params.values()), t_eval=t_eval) results.append(sol.y[:, -1]) # 记录最终种群数量 # 可视化结果 plt.figure(figsize=(10, 6)) values = np.array(values) results = np.array(results) plt.plot(values, results[:, 0], 'b-', label='兔子最终数量') plt.plot(values, results[:, 1], 'g-', label='羊最终数量') plt.xlabel(param_name) plt.ylabel('最终种群数量') plt.title(f'{param_name}参数敏感性分析') plt.legend(); plt.grid(True) plt.show() # 示例：分析alpha12的影响 analyze_sensitivity(params, 'alpha12', np.linspace(0.1, 2, 20))5.2 三维可视化
展示参数变化如何影响竞争结果：
from mpl_toolkits.mplot3d import Axes3D def plot_3d_phase_diagram(): """创建3D相位图""" fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') # 模拟不同初始条件下的轨迹 for N1_0 in np.linspace(50, 200, 5): for N2_0 in np.linspace(50, 200, 5): sol = solve_ivp(competition_model, t_span, [N1_0, N2_0], args=tuple(params.values()), t_eval=t_eval) ax.plot(sol.y[0], sol.y[1], sol.t, alpha=0.6) ax.set_xlabel('兔子数量') ax.set_ylabel('羊数量') ax.set_zlabel('时间') ax.set_title('种群竞争3D相位图') plt.show() plot_3d_phase_diagram()5.3 交互式探索工具
使用IPyWidgets创建交互式界面：
from ipywidgets import interact, FloatSlider def interactive_simulation(r1, r2, K1, K2, alpha12, alpha21): """交互式模拟函数""" params = locals() sol = solve_ivp(competition_model, t_span, N0, args=tuple(params.values()), t_eval=t_eval) plt.figure(figsize=(10, 6)) plt.plot(sol.t, sol.y[0], label='兔子') plt.plot(sol.t, sol.y[1], label='羊') plt.xlabel('时间'); plt.ylabel('种群数量') plt.title('交互式竞争模拟') plt.legend(); plt.grid(True) plt.show() # 创建交互控件 interact( interactive_simulation, r1=FloatSlider(min=0.1, max=1, step=0.05, value=0.5), r2=FloatSlider(min=0.1, max=1, step=0.05, value=0.4), K1=FloatSlider(min=100, max=1000, step=50, value=500), K2=FloatSlider(min=100, max=1000, step=50, value=600), alpha12=FloatSlider(min=0.1, max=2, step=0.1, value=0.7), alpha21=FloatSlider(min=0.1, max=2, step=0.1, value=0.5) );6. 实际应用中的注意事项
在将模型应用于实际问题时，有几个关键点需要考虑：
参数估计：模型参数需要通过实地观察或实验数据来估计。例如：
- 内禀增长率(r)可以通过种群在理想条件下的增长数据估计
- 环境承载力(K)可以通过资源可用性计算
- 竞争系数(α)需要通过物种间的相互作用实验确定
模型局限性：
- 假设环境条件是恒定的
- 忽略空间异质性
- 不考虑进化适应
- 假设竞争是即时发生的
验证方法：
- 将模型预测与实际观测数据比较
- 进行敏感性分析，了解哪些参数影响最大
- 尝试不同的初始条件，检查结果的稳健性
提示：在实际生态研究中，Lotka-Volterra模型常作为起点，更复杂的模型会考虑更多因素，如空间分布、年龄结构、环境随机性等。
