如何用Scipy odeint返回普通变量？将DeSolve的R代码转Python

Scipy的odeint和R的deSolve在返回额外变量的机制上有差异：deSolve允许导数函数同时返回微分和额外变量，而odeint仅专注于求解状态变量的积分结果。要获取Total、rate、Conserved这类派生变量，可按以下两种方式实现：

方法1：积分后逐点计算额外变量 #

如果所有额外变量都是状态变量和时间的函数，最直接的方式是先求解状态变量，再通过辅助函数批量计算每个时间点的派生变量。

示例代码 #

假设你的R模型逻辑如下（模拟）：

model <- function(t, y, parms) {
  A <- y[1]
  B <- y[2]
  k1 <- parms[1]
  k2 <- parms[2]
  
  rate1 <- k1 * A
  rate2 <- k2 * B
  dA <- -rate1 + rate2
  dB <- rate1 - rate2
  
  Total <- A + B
  Conserved <- A + 2*B  # 示例守恒量
  return(list(c(dA, dB), Total=Total, rate=rate1, Conserved=Conserved))
}

对应的Python实现：

import numpy as np
from scipy.integrate import odeint

# ODE导数函数：仅返回状态变量的微分
def dydt(y, t, k1, k2):
    A, B = y
    rate1 = k1 * A
    rate2 = k2 * B
    dA = -rate1 + rate2
    dB = rate1 - rate2
    return [dA, dB]

# 辅助函数：计算指定状态下的额外变量
def get_extra_vars(y, k1, k2):
    A, B = y
    rate1 = k1 * A
    rate2 = k2 * B
    return {
        "Total": A + B,
        "rate": np.array([rate1, rate2]),
        "Conserved": A + 2*B  # 完全照搬R中Conserved的计算逻辑
    }

# 初始化参数
y0 = [10, 0]  # 初始状态
params = (0.1, 0.05)  # 模型参数
t = np.linspace(0, 50, 100)  # 时间点

# 求解ODE
y_sol = odeint(dydt, y0, t, args=params)

# 批量计算所有时间点的额外变量
Total = np.array([get_extra_vars(y, *params)["Total"] for y in y_sol])
Conserved = np.array([get_extra_vars(y, *params)["Conserved"] for y in y_sol])
rate = np.array([get_extra_vars(y, *params)["rate"] for y in y_sol])

方法2：使用scipy.integrate.ode实时记录变量 #

如果需要在积分过程中实时跟踪额外变量（而非事后批量计算），可以使用更灵活的ode类，在导数函数中直接计算并存储派生变量。

示例代码 #

import numpy as np
from scipy.integrate import ode

# 初始化存储额外变量的容器
extra_vars = {"Total": [], "Conserved": [], "rate": []}

def dydt(t, y, k1, k2):
    A, B = y
    rate1 = k1 * A
    rate2 = k2 * B
    dA = -rate1 + rate2
    dB = rate1 - rate2
    
    # 计算并存储额外变量
    extra_vars["Total"].append(A + B)
    extra_vars["Conserved"].append(A + 2*B)
    extra_vars["rate"].append([rate1, rate2])
    
    return [dA, dB]

# 初始化求解器
solver = ode(dydt).set_integrator("vode", method="bdf")
solver.set_initial_value(y0, t[0])
solver.set_f_params(*params)

# 逐步积分
for ti in t[1:]:
    solver.integrate(ti)
    if not solver.successful():
        break

# 转换为NumPy数组方便后续处理
extra_vars["Total"] = np.array(extra_vars["Total"])
extra_vars["Conserved"] = np.array(extra_vars["Conserved"])
extra_vars["rate"] = np.array(extra_vars["rate"])
# 获取状态变量的解
y_sol = np.vstack([y0, solver.y])

关于Conserved变量的关键注意事项 #

Conserved作为守恒量，其计算逻辑完全由你的模型决定：

• 如果是状态变量的线性/非线性组合（比如sum(y)或2*y[0] + 3*y[1]），直接照搬R代码中的表达式即可；
• 如果是速率的累积值（比如对某个速率积分），可以用np.trapz在时间轴上对速率变量积分得到；
• 守恒量理论上在整个积分过程中数值保持稳定（仅受浮点误差影响），可通过绘制Conserved随时间的变化曲线验证正确性。

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