使用sunode独立求解常微分方程：以Lotka-Volterra模型为例-优快云博客

本文链接：https://blog.youkuaiyun.com/gitblog_00424/article/details/148682384

使用sunode独立求解常微分方程：以Lotka-Volterra模型为例

sunode Solve ODEs fast, with support for PyMC 项目地址: https://gitcode.com/gh_mirrors/su/sunode

项目概述

sunode是一个用于求解常微分方程(ODE)的Python库，它提供了强大的符号计算和数值求解能力。本文将重点介绍如何在不依赖PyMC的情况下，使用sunode独立求解常微分方程系统。

Lotka-Volterra模型简介

Lotka-Volterra模型是生态学中描述捕食者-猎物关系的经典模型。它由两个微分方程组成：

野兔(Hares)数量变化：与自身数量成正比增长，被猞猁捕食而减少
猞猁(Lynxes)数量变化：捕食野兔而增长，自然死亡而减少

数学模型表示为：

dH/dt = αH - βLH
dL/dt = δLH - γL

使用sunode求解步骤

1. 定义参数和状态变量

首先需要声明模型中的参数和状态变量及其形状：

params = {
    'alpha': (),  # 野兔自然增长率
    'beta': (),   # 野兔被捕食率
    'gamma': (),  # 猞猁自然死亡率
    'delta': (),  # 猞猁捕食增长率
}

states = {
    'hares': (),  # 野兔数量
    'lynxes': (), # 猞猁数量
}

2. 定义微分方程右侧函数

这是模型的核心部分，需要返回每个状态变量的导数：

def lotka_volterra(t, y, p):
    """Lotka-Volterra方程的右侧函数"""
    return {
        'hares': p.alpha * y.hares - p.beta * y.lynxes * y.hares,
        'lynxes': p.delta * y.hares * y.lynxes - p.gamma * y.lynxes,
    }

3. 对数变换版本

有时需要对变量进行对数变换来改善数值稳定性：

import sympy as sym

def lotka_volterra_log(t, y, p):
    exp = np.vectorize(sym.exp)  # 向量化sympy的exp函数
    
    hares = exp(y.log_hares)
    lynxes = exp(y.log_lynxes)
    
    dhares = p.alpha * hares - p.beta * lynxes * hares
    dlynxes = p.delta * hares * lynxes - p.gamma * lynxes
    return {
        'log_hares': dhares / hares,  # d(logH)/dt = (dH/dt)/H
        'log_lynxes': dlynxes / lynxes,
    }

4. 创建SympyProblem实例

将上述定义封装成问题实例：

problem = sunode.SympyProblem(
    params=params,
    states=states,
    rhs_sympy=lotka_volterra,
    derivative_params=()
)

5. 选择求解器并求解

sunode提供多种求解器选项：

# 创建不计算敏感度的求解器
solver = sunode.solver.Solver(problem, compute_sens=False, solver='BDF')

# 设置初始条件
y0 = np.zeros((), dtype=problem.state_dtype)
y0['hares'] = 1
y0['lynxes'] = 0.1

# 设置时间点
tvals = np.linspace(0, 10)

# 设置参数值
solver.set_params_dict({
    'alpha': 0.1,
    'beta': 0.2,
    'gamma': 0.3,
    'delta': 0.4,
})

# 求解并输出结果
output = solver.make_output_buffers(tvals)
solver.solve(t0=0, tvals=tvals, y0=y0, y_out=output)

6. 结果可视化

可以将结果转换为xarray Dataset或直接访问numpy数组：

# 使用xarray绘图
solver.as_xarray(tvals, output).solution_hares.plot()

# 或直接使用matplotlib绘图
plt.plot(output.view(tvals, problem.state_dtype)['hares'])

注意事项

符号计算特性：右侧函数中的变量是sympy符号而非具体数值，因此不能直接使用Python的条件判断。例如：
```
# 错误示例 - 不会按预期工作
value = 1
if y.some_state > 1:  # y.some_state是符号变量
    value += 1
```
性能考虑：使用结构化数组不会引入运行时开销，同时提高了代码可读性。
求解器选择：对于需要计算梯度的情况，可以使用AdjointSolver。