HOWTO · Python
在 Python 中模仿 ode45() 函数
本文讨论了 ode45() 函数在 python 中的实现。最初，该函数是在 MatLab 中定义的。但是，本文讨论了在 Python 中使用 scipy 模块实现 ode45() 函数。
本页内容
MATLAB 中的 ode45() 函数使用自适应 Runge-Kutta 方法求解常微分方程（ODE）的初值问题。Python 中没有直接名为 ode45() 的函数，但 SciPy 的 solve_ivp() 涵盖了同样的工作流程：你描述方程，为求解器给出时间区间和初始状态，然后从它返回的对象中读取解。
使用 SciPy 的 solve_ivp() 代替 MATLAB 的 ode45()
scipy.integrate 中的 solve_ivp() 函数是与 ode45() 最接近的实用对应物。它的默认方法 RK45 是一个 5(4) 阶显式 Runge-Kutta 公式，带自动步长控制——与 MATLAB 中 ode45() 使用的算法同属一族，因此求解器在解变化缓慢处采用较大步长，在解变化剧烈处采用较小步长。
它并不是可直接替换的等价物：参数名称不同，结果以单个对象返回，而不是一对 t 和 y 数组。不过整体模式——导数函数、时间跨度和初始状态——可以直接对应过去。如果尚未安装 SciPy，可以用 pip install scipy 安装。
定义 ODE 并调用 solve_ivp()
最基本的调用形式是 solve_ivp(fun, t_span, y0, method="RK45", t_eval=None, args=None, ...)。主要输入为：
fun：可调用对象fun(t, y)，针对当前时间t和状态y返回导数dy/dt。单方程时返回一个值；方程组时返回一个数组，每个状态分量对应一个导数。t_span：二元组(t_start, t_end)，给出积分区间的起点和终点。y0：t_start时刻的状态；其长度决定方程组中方程的数量。method（可选）：积分方案，默认为'RK45'。args（可选）：以fun(t, y, *args)形式传递给函数的额外参数。rtol和atol（可选）：控制自适应步长的相对容差和绝对容差。
下面的例子在区间 [0, 2] 上从 y(0) = 1 开始求解 dy/dt = 2y。调用返回一个结果对象；打印其 t 和 y 属性可以看到求解器选取的时间点以及每个时间点上的状态：
本节中的代码以及显示的标签或输出仍为英文。
from scipy.integrate import solve_ivp
def growth(t, y):
return 2 * y
result = solve_ivp(growth, (0, 2), [1])
print("t:", result.t)
print("y:", result.y[0])
t: [0. 0.07580098 0.56590702 1.18757739 1.85907272 2. ]
y: [ 1. 1.16369695 3.10133569 10.75272756 41.18553629 54.59497095]
报告时间的间距并不均匀，因为 solve_ivp() 会根据解的局部行为调整每一步，而不是按固定增量推进。
读取求解器结果
solve_ivp() 返回一个 scipy.integrate.OdeResult 对象。最常用的属性有：
t：报告时间点的 NumPy 数组。未传t_eval时，这些就是求解器内部的自适应步点。y：形状为(n, len(t))的二维 NumPy 数组；第i行保存状态第i个分量在各报告时间点的值，因此y[:, k]是t[k]时刻的状态。单方程时，y[0]就是解曲线。sol：连续的OdeSolution插值对象，可用result.sol(t)求区间内任意时刻的解。只有传入dense_output=True时才会填充，否则为None。success和message：前者是表示积分是否完成的布尔值，后者是对结束方式的可读描述。status：区间完成时为0，发生终止事件时为1，求解器无法继续推进时为-1。
在使用 y 之前先检查 success（或 status）是个好习惯，因为失败的运行仍可能返回部分数组。
用 t_eval 选择输出时间
自适应步点很少恰好是你要报告的时间。把一组排序后的时间传给 t_eval，solve_ivp() 就会精确地在这些时间返回解，而不是返回内部步点：
本节中的代码以及显示的标签或输出仍为英文。
from scipy.integrate import solve_ivp
def growth(t, y):
return 2 * y
times = [0.5, 1.0, 1.5, 2.0]
result = solve_ivp(growth, (0, 2), [1], t_eval=times)
print("t:", result.t)
print("y:", result.y[0])
t: [0.5 1. 1.5 2. ]
y: [ 2.71823524 7.38683161 20.07910026 54.59497095]
此时 t 只包含所请求的时间，y 保存每个时间点上插值得到的状态。t_eval 中的每个值都必须位于 t_span 之内；请求积分区间之外的时间会使调用报错，提示这些值必须在 t_span 范围内。端点是允许的，因此 t_eval 可以同时包含 t_start 和 t_end。
