如何使用scipy.integrate.solve_ivp稠密输出直接评估二维数组上的解

问题原因 #

scipy.integrate.solve_ivp返回的OdeSolution插值对象的调用逻辑默认仅支持标量或一维数组作为时间输入，内部判断输入是否为标量的代码遇到二维数组时会触发数组布尔值判断的歧义错误，就是你遇到的ValueError报错。

最优实现方案 #

你目前使用的ravel展平+形状还原的方案本身就是效率最高的标准实现：

• numpy的ravel对于连续存储的数组仅修改视图属性，不会产生数据拷贝开销
• 后续的reshape操作同样仅修改数组的维度元数据，几乎没有额外性能损耗
• OdeSolution内部本身就是对输入的所有时间点做批量向量化计算，展平操作完全没有浪费插值的并行性能

可以直接封装为通用写法，适配任意维度的时间数组：

# t_arr为任意维度的时间数组，输出解的形状为 (解的维度, *t_arr.shape)
sol_res = res.sol(t_arr.ravel()).reshape(-1, *t_arr.shape)
# 如果是单变量初值问题，只需取第一维的话可以写为：
# sol_res = res.sol(t_arr.ravel()).reshape(t_arr.shape)

没有必要寻找其他更复杂的实现，当前方案的性能已经达到理论上限，同时通用性强，支持二维及更高维度的时间数组输入。

内容的提问来源于stack exchange，提问作者amzon-ex