如何从scipy.odeint中返回非微分类型的额外值？

问题根源 #

odeint采用自适应步长算法，为保证求解精度，会在内部多次调用你的model函数（调用次数远多于你指定的t数组长度），直接在model里追加值自然会导致列表长度远超预期。

实用解决方案 #

方案1：仅在指定输出时间点存储额外值 #

通过判断当前调用model时的t是否属于你定义的输出时间点（注意浮点数比较要加容差），只在这些时刻存储值：

import numpy as np
from scipy.integrate import odeint

def model(y, t, target_times, tol=1e-8):
    dydt = -2 * y
    x = np.random.rand(1,1)
    # 检查当前时间是否在目标输出时间点中（规避浮点数误差）
    if np.any(np.abs(t - target_times) < tol):
        addtional_value.append(x)
    return dydt

y0 = 5
addtional_value = []
t = np.arange(0, 5, 1)

# 把目标时间点传给model函数
y = odeint(model, y0, t, args=(t,))

print(len(y))          # 输出5
print(len(addtional_value))  # 输出5

方案2：将额外值打包进状态变量 #

把不需要微分的额外值作为状态的一部分，设置它们的导数为0，求解后再分离出需要的部分：

import numpy as np
from scipy.integrate import odeint

def model(y, t):
    # y[0]是原状态变量，y[1]是额外值
    dydt = np.zeros_like(y)
    dydt[0] = -2 * y[0]
    dydt[1] = 0  # 额外值的导数设为0，保持状态稳定
    y[1] = np.random.rand()  # 更新当前额外值
    return dydt

y0 = np.array([5, 0])  # 初始状态：原变量+额外值初始占位
t = np.arange(0, 5, 1)

y_full = odeint(model, y0, t)
y = y_full[:, 0]  # 提取原ODE的解
addtional_value = y_full[:, 1]  # 提取额外值

print(len(y))  # 输出5
print(len(addtional_value))  # 输出5

方案3：改用solve_ivp（推荐，更灵活） #

如果可以切换到scipy.integrate.solve_ivp，它支持自定义输出回调，能精准控制何时存储额外值：

import numpy as np
from scipy.integrate import solve_ivp

addtional_value = []

def model(t, y):
    dydt = -2 * y
    return dydt

def store_extra(t, y):
    x = np.random.rand(1,1)
    addtional_value.append(x)
    return 0  # 回调返回值不影响求解逻辑

y0 = [5]
t_span = [0, 4]
t_eval = np.arange(0, 5, 1)

# 用t_eval指定输出点，回调函数在每个t_eval点触发
sol = solve_ivp(model, t_span, y0, t_eval=t_eval, events=store_extra)

print(len(sol.y[0]))  # 输出5
print(len(addtional_value))  # 输出5

方案说明 #

• 方案1适合只需要在指定输出时间点获取额外值的场景，注意容差设置要合理，避免漏判或误判。
• 方案2不需要额外传参，但会增加状态变量维度，适合额外值数量不多的情况。
• 方案3的solve_ivp是scipy官方推荐的新一代求解器，功能更丰富，回调机制更可靠。

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