可在积分步修改y向量的Python ODE积分器选型咨询

需要求解一组常微分方程（ODEs），但无法使用scipy中的现有求解器——因为每次成功完成积分步骤后需要修改求解器中的y值，而scipy的求解器中y是由内部控制的，直接修改solver.y无效。现寻求具备该功能的Python工具包。

以下是基于scipy的无效设想代码：

solver = integrate.LSODA(fun=psr_model, t0=0, y0=y0, t_bound=simulation_time, first_step=1e-10, max_step=1e-3)
while solver.t < simulation_time:
    solver.step()
    # 该语句无效，scipy不允许直接修改内部y值
    solver.y = correct_isentropic(P_meas(solver.t), solver.y)

可行工具包及实现方式 #

1. PyDSTool #

PyDSTool专为动力系统建模设计，支持完全自定义的积分流程，允许在每步积分后手动修改状态变量。你可以通过构建步进循环，在每次积分完成后直接更新状态向量。

示例代码框架：

import PyDSTool as dst

# 定义ODE模型
dsargs = dst.args(name='psr_model')
dsargs.pars = {}  # 传入模型参数
dsargs.varspecs = {
    # 按你的ODE定义变量导数，比如 'y1': 'f(y1, y2, t)', ...
}
dsargs.ics = {'y1': y0[0], 'y2': y0[1]}  # 初始条件
dsargs.tdomain = [0, simulation_time]

# 创建模型对象
model = dst.Generator.Vode_ODEsystem(dsargs)
model.set(pars={})  # 设置参数

# 自定义步进循环
t_current = 0.0
y_current = y0.copy()
while t_current < simulation_time:
    # 积分一步到下一个时间点
    model.set(ics={**{f'y{i+1}': y_current[i] for i in range(len(y0))}, 't': t_current})
    traj = model.compute('step', tspan=[t_current, t_current + 1e-3])  # 这里可以控制步长
    t_next = traj['t'][-1]
    y_next = [traj[f'y{i+1}'][-1] for i in range(len(y0))]
    
    # 修改y值
    y_next = correct_isentropic(P_meas(t_next), y_next)
    
    # 更新当前状态
    t_current = t_next
    y_current = y_next

2. 自定义显式求解器（搭配Numba加速） #

如果追求最高自由度，完全手动实现显式求解器（比如RK4、欧拉法），可以在每步积分后直接修改y值，没有任何限制。用Numba可以把性能拉到接近C语言水平。

示例代码（RK4方法）：

import numba as nb
import numpy as np

# 用Numba加速ODE右侧函数
@nb.jit(nopython=True)
def psr_model_numba(t, y):
    # 实现你的ODE导数计算，返回dy/dt数组
    dy = np.zeros_like(y)
    # dy[0] = ..., dy[1] = ...
    return dy

# 用Numba加速RK4步进函数
@nb.jit(nopython=True)
def rk4_step(t, y, dt, fun):
    k1 = fun(t, y)
    k2 = fun(t + dt/2, y + dt*k1/2)
    k3 = fun(t + dt/2, y + dt*k2/2)
    k4 = fun(t + dt, y + dt*k3)
    return y + dt*(k1 + 2*k2 + 2*k3 + k4)/6

# 自定义积分循环
t = 0.0
y = y0.copy()
dt = 1e-10  # 初始步长
max_dt = 1e-3

while t < simulation_time:
    # 限制步长不超过剩余时间
    current_dt = min(dt, simulation_time - t)
    # RK4步进
    y = rk4_step(t, y, current_dt, psr_model_numba)
    t += current_dt
    
    # 修改y值
    y = correct_isentropic(P_meas(t), y)
    
    # 可以根据需要调整步长（比如误差控制）
    dt = min(dt * 1.1, max_dt)

3. 改造SciPy求解器的使用方式（无需换工具） #

其实不用换工具包，而是把每步积分后的修正作为新的初始条件，重新初始化求解器。虽然效率略低，但能快速适配现有代码：

示例代码：

from scipy import integrate

t_current = 0.0
y_current = y0.copy()
max_step = 1e-3

while t_current < simulation_time:
    # 计算下一步的时间边界
    t_next_bound = min(t_current + max_step, simulation_time)
    # 初始化求解器，从当前状态开始积分到t_next_bound
    solver = integrate.LSODA(
        fun=psr_model,
        t0=t_current,
        y0=y_current,
        t_bound=t_next_bound,
        first_step=1e-10,
        max_step=max_step
    )
    # 执行积分到边界
    while solver.t < t_next_bound:
        solver.step()
    # 获取积分结果
    t_current = solver.t
    y_current = solver.y.copy()
    
    # 修改y值
    y_current = correct_isentropic(P_meas(t_current), y_current)

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