自定义欧拉法求解微分方程组始终发散至无穷的原因分析

问题背景 #

针对如下捕食者-猎物微分方程组：
$$
\begin{cases}
\frac{dh}{dt} = a h - b h l \
\frac{dl}{dt} = d h l - g l
\end{cases}
$$
其中参数 $a=0.63, b=0.02, g=0.61, d=0.015$，初始条件 $h(0)=24.2, l(0)=16.16$。使用自定义显式欧拉法求解时结果趋于无穷，但用scipy.integrate.odeint却能得到正确的周期解，最终通过减小步长解决了发散问题。

一、odeint与显式欧拉法的核心差异 #

• 算法类型与自适应能力：odeint默认采用LSODA算法，这是一种自适应步长的多步法，会根据当前解的误差动态调整步长：解变化平缓时用大步长提升效率，变化剧烈时自动切换为小步长保证稳定性；同时它能在显式与隐式方法间自动切换，适配不同类型的微分方程。
• 显式欧拉法的固有局限：显式欧拉是一阶单步法，精度低且稳定区域有限——只有当步长小于某个阈值时，才能保证数值解不发散。对于非线性系统（如捕食者-猎物模型），这个阈值通常较小。

二、自定义欧拉法发散的原因 #

1. 步长过大超出稳定区域：你的自定义实现中，步长 $h=1$（每次迭代直接推进1年），而该系统下显式欧拉的稳定步长阈值小于1。步长超过阈值后，数值误差会在迭代中不断被放大，导致解偏离真实轨道并最终发散到无穷。
2. 状态变量传递不严谨：调用f(result[i], i)时，传入的第一个参数是包含年份的列表[year, h, l]，虽在f中正确取到了h和l，但这种传递方式容易引发混淆，增加调试难度。

对比来看，odeint的自适应步长会自动选择远小于1的步长，因此能稳定追踪真实的周期解。

三、修正方案 #

1. 减小步长 #

将大步长拆分为多个小步迭代，例如设置步长 $h=0.1$，每年拆分为10个小步计算：

import numpy as np

def f(s, t):
    a = 0.63
    b = 0.02
    g = 0.61
    d = 0.015
    h, l = s  # 仅传递状态变量[h,l]，避免混淆
    dhdt = a * h - b * h * l
    dldt = d * h * l - g * l
    return [dhdt, dldt]

result = [[1890, 24.2, 16.16]]
h = 0.1  # 减小步长
steps_per_year = int(1 / h)

for year in range(40):
    current_h, current_l = result[-1][1], result[-1][2]
    # 每年拆分为多个小步计算
    for _ in range(steps_per_year):
        dhdt, dldt = f([current_h, current_l], 0)
        current_h += h * dhdt
        current_l += h * dldt
    result.append([1890 + 1 + year, current_h, current_l])

# 输出结果
for item in result:
    print(item)

2. 规范状态变量传递 #

调用f时仅传递状态变量[h,l]，而非包含年份的完整列表，避免索引错误。

总结 #

自定义欧拉法发散的核心原因是步长超出显式欧拉的稳定区域，而odeint的自适应步长算法从根本上避免了这个问题；同时规范状态变量的传递方式能降低调试风险。减小步长后，显式欧拉法可以得到与真实解一致的周期结果。

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