用Python模拟兔子和羊的‘地盘争夺战’：从零实现Lotka-Volterra竞争模型

用Python模拟兔子和羊的‘地盘争夺战’从零实现Lotka-Volterra竞争模型想象你有一片草地同时放养了兔子和羊。起初它们相安无事但随着种群增长争夺草资源的矛盾逐渐显现——有的结局是兔子灭绝有的是羊消失偶尔也会出现两者共存的奇妙平衡。这种动态博弈可以用Lotka-Volterra竞争模型精确描述而今天我们将用Python让它活起来。1. 环境准备与模型原理1.1 工具链配置推荐使用Anaconda创建独立环境避免依赖冲突conda create -n ecology python3.9 conda activate ecology pip install numpy matplotlib scipy jupyter1.2 竞争模型核心方程Lotka-Volterra模型通过微分方程描述两个物种的竞争关系$$ \begin{cases} \frac{dN_1}{dt} r_1N_1\left(1 - \frac{N_1 \alpha_{12}N_2}{K_1}\right) \ \frac{dN_2}{dt} r_2N_2\left(1 - \frac{N_2 \alpha_{21}N_1}{K_2}\right) \end{cases} $$参数说明$N_1, N_2$兔子和羊的种群数量$r_1, r_2$固有增长率$K_1, K_2$环境承载力$\alpha_{12}$羊对兔子的竞争系数$\alpha_{21}$兔子对羊的竞争系数关键提示当α1表示竞争激烈1则表示影响较弱2. Python实现步骤2.1 微分方程求解器使用Scipy的odeint进行数值求解import numpy as np from scipy.integrate import odeint def competition_model(y, t, r1, r2, K1, K2, alpha12, alpha21): N1, N2 y dN1dt r1 * N1 * (1 - (N1 alpha12 * N2) / K1) dN2dt r2 * N2 * (1 - (N2 alpha21 * N1) / K2) return [dN1dt, dN2dt]2.2 参数配置与模拟设置四组典型参数演示不同竞争结果情景r1r2K1K2α12α21预期结果10.80.65004001.20.8兔子胜出20.70.93005000.71.3羊胜出31.01.04004000.50.5稳定共存40.90.93503501.11.1不确定# 模拟参数 t np.linspace(0, 50, 1000) initial_pop [100, 100] # 初始种群数量 # 情景1兔子胜出 params_case1 {r1:0.8, r2:0.6, K1:500, K2:400, alpha12:1.2, alpha21:0.8} solution odeint(competition_model, initial_pop, t, argstuple(params_case1.values()))3. 结果可视化与分析3.1 动态过程绘图使用Matplotlib制作专业级图表import matplotlib.pyplot as plt def plot_results(t, N1, N2, case_name): plt.figure(figsize(12, 5)) plt.plot(t, N1, r-, labelRabbit Population) plt.plot(t, N2, b-, labelSheep Population) plt.title(fCompetition Outcome: {case_name}) plt.xlabel(Time (days)) plt.ylabel(Population Size) plt.legend() plt.grid(True) plt.show()3.2 相空间分析通过零增长线dN/dt0预测长期趋势def plot_phase_portrait(K1, K2, alpha12, alpha21): # 计算零增长线交点 N1_range np.linspace(0, K1*1.5, 100) N2_isocline (K1 - N1_range) / alpha12 N2_range np.linspace(0, K2*1.5, 100) N1_isocline (K2 - N2_range) / alpha21 plt.plot(N1_range, N2_isocline, r-, labelRabbit Isocline) plt.plot(N1_isocline, N2_range, b-, labelSheep Isocline) plt.axis([0, max(K1, K2)*1.2, 0, max(K1/alpha12, K2)*1.2]) plt.xlabel(Rabbit Population) plt.ylabel(Sheep Population) plt.legend()4. 模型扩展与应用4.1 参数敏感性分析使用SALib库进行全局敏感性分析from SALib.analyze import sobol problem { num_vars: 6, names: [r1, r2, K1, K2, alpha12, alpha21], bounds: [[0.5, 1.5], [0.5, 1.5], [300, 700], [300, 700], [0.3, 1.7], [0.3, 1.7]] } # 生成参数样本并进行模拟后... Si sobol.analyze(problem, simulation_results)4.2 商业竞争模拟将生态模型迁移到市场分析def market_competition(y, t, params): prod1, prod2 y growth1 params[r1] * prod1 * (1 - (prod1 params[alpha12]*prod2)/params[K1]) growth2 params[r2] * prod2 * (1 - (prod2 params[alpha21]*prod1)/params[K2]) return [growth1, growth2]关键参数对应关系种群数量 → 产品市场份额环境承载力 → 市场总容量竞争系数 → 产品替代性在最近的一个咨询项目中我用这个模型成功预测了两款同类APP的用户增长拐点。当把α12调整为0.8表示功能重叠度时模型准确复现了实际观测到的老二逆袭现象。