如何使用scipy.derivative处理solve_ivp输出的OdeSolution对象？

问题 #

我需要计算常微分方程数值解的残差，相关代码如下：

import numpy as np
from scipy.integrate import solve_ivp

def f(t, x): 
    return 0.5 * x

t_start = 0
t_end = 2

n = 50

t = np.linspace(t_start, t_end, 50)

x_init = 1

solution = solve_ivp(f, [t_start, t_end], [x_init], dense_output = True)

尝试用scipy的derivative函数在t点对结果做数值微分，直接调用：

from scipy.differentiate import derivative

derivative(solution.sol, t)

出现尺寸不匹配错误。之后定义了：

def g(t):
    return solution.sol(t)[0]

derivative(g, t)

仍无法正常运行。仅想了解如何使用scipy的derivative函数处理solve_ivp的输出。

解决方法 #

核心错误点 #

1. 模块导入错误：不存在scipy.differentiate，正确的derivative函数在scipy.misc模块中。
2. 函数输入输出维度不匹配：solution.sol返回的是二维数组，derivative默认要求函数是标量输入→标量输出的映射，直接传入数组会触发维度问题。
3. 未启用批量计算：derivative默认只处理标量，需要指定参数让它支持数组输入的批量处理。

修正后的代码 #

import numpy as np
from scipy.integrate import solve_ivp
from scipy.misc import derivative  # 修正模块导入路径

def f(t, x): 
    return 0.5 * x

t_start = 0
t_end = 2
n = 50
t = np.linspace(t_start, t_end, n)
x_init = 1

solution = solve_ivp(f, [t_start, t_end], [x_init], dense_output=True)

# 定义标量输入→标量输出的映射函数
def g(t_val):
    # 当输入是标量时，sol返回(1,1)二维数组，提取标量值
    return solution.sol(t_val)[0][0]

# 调用derivative，启用vectorize=True支持数组批量计算
dx_dt = derivative(g, t, vectorize=True)

# 计算残差：数值导数 - 原方程右端项
residual = dx_dt - 0.5 * solution.sol(t)[0]
print(residual)

关键说明 #

• solution.sol(t_val)在输入为标量时，返回形状为(1,1)的二维数组，必须用[0][0]提取标量，确保g函数符合derivative的输入输出要求。
• vectorize=True参数会让derivative自动遍历数组t中的每个元素，逐个计算导数，彻底解决维度不匹配问题。

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