solve_ivp#
- scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)[source]#
求解常微分方程组的初值问题。
此函数在给定初值的情况下，对常微分方程组进行数值积分。
dy / dt = f(t, y) y(t0) = y0
这里 t 是一个一维自变量（时间），y(t) 是一个 N 维向量值函数（状态），N 维向量值函数 f(t, y) 确定了微分方程。目标是找到近似满足微分方程的 y(t)，给定初值 y(t0)=y0。
一些求解器支持复数域上的积分，但请注意，对于刚性 ODE 求解器，右侧必须是复数可微分的（满足柯西-黎曼方程 [11]）。要在复数域中求解问题，请传递具有复数数据类型的 y0。另一种始终可用的选项是分别重写实部和虚部的问题。
- 参数:
- fun可调用对象
系统的右侧：状态
y在时间t处的导数。调用签名为fun(t, y)，其中t是标量，y是一个 ndarray，其len(y) = len(y0)。如果使用args（请参阅args参数的文档），则需要传递额外参数。fun必须返回与y形状相同的数组。有关更多信息，请参阅 vectorized。- t_span2 成员序列
积分区间 (t0, tf)。求解器从 t=t0 开始积分，直到达到 t=tf。t0 和 tf 都必须是浮点数或可由浮点转换函数解释的值。
- y0类数组, 形状 (n,)
初始状态。对于复数域中的问题，即使初始值为纯实数，也请传递具有复数数据类型的 y0。
- method字符串或
OdeSolver, 可选 要使用的积分方法
‘RK45’ (默认): 5(4) 阶显式龙格-库塔方法 [1]。误差控制假设四阶方法的精度，但使用五阶精确公式进行步进（进行局部外推）。对于密集输出，使用四次插值多项式 [2]。可应用于复数域。
‘RK23’: 3(2) 阶显式龙格-库塔方法 [3]。误差控制假设二阶方法的精度，但使用三阶精确公式进行步进（进行局部外推）。对于密集输出，使用三次 Hermite 多项式。可应用于复数域。
‘DOP853’: 8 阶显式龙格-库塔方法 [13]。Python 实现的“DOP853”算法最初用 Fortran 编写 [14]。对于密集输出，使用精确到 7 阶的 7 阶插值多项式。可应用于复数域。
‘Radau’: 5 阶 Radau IIA 族隐式龙格-库塔方法 [4]。误差通过三阶精确嵌入公式控制。对于密集输出，使用满足搭配条件的三次多项式。
‘BDF’: 基于导数近似的反向微分公式的隐式多步变阶（1 到 5）方法 [5]。实现遵循 [6] 中描述的方法。使用准恒定步长方案，并通过 NDF 修改提高精度。可应用于复数域。
‘LSODA’: 具有自动刚度检测和切换的 Adams/BDF 方法 [7], [8]。这是 ODEPACK 中 Fortran 求解器的包装器。
显式龙格-库塔方法（'RK23'、'RK45'、'DOP853'）应用于非刚性问题，隐式方法（'Radau'、'BDF'）应用于刚性问题 [9]。在龙格-库塔方法中，建议使用 'DOP853' 进行高精度求解（rtol 和 atol 值低）。
如果不确定，首先尝试运行 'RK45'。如果它迭代次数异常多、发散或失败，则您的问题很可能是刚性问题，您应该使用 'Radau' 或 'BDF'。'LSODA' 也可以是一个很好的通用选择，但由于它包装了旧的 Fortran 代码，使用起来可能不那么方便。
您也可以传递任何派生自
OdeSolver并实现求解器的任意类。- t_eval类数组或 None, 可选
存储计算结果的时间点，必须排序并在 t_span 范围内。如果为 None（默认），则使用求解器选择的点。
- dense_output布尔值, 可选
是否计算连续解。默认为 False。
- events可调用对象，或可调用对象列表，可选
要跟踪的事件。如果为 None（默认），则不跟踪任何事件。每个事件都发生在时间连续函数和状态的零点处。每个函数必须具有签名
event(t, y)，如果使用了args（请参阅args参数的文档），则必须传递额外参数。每个函数必须返回一个浮点数。求解器将使用寻根算法找到event(t, y(t)) = 0的精确 t 值。默认情况下，将找到所有零点。求解器在每一步中寻找符号变化，因此如果一步内发生多个零点穿越，则可能会错过事件。此外，每个 event 函数可能具有以下属性：- terminal: 布尔值或整数, 可选
当为布尔值时，表示如果发生此事件是否终止积分。当为整数时，终止发生在此事件指定次数后。如果未指定，则隐式为 False。
- direction: 浮点数, 可选
零点穿越的方向。如果 direction 为正，则 event 仅在从负到正时触发，反之亦然，如果 direction 为负。如果为 0，则任一方向都会触发事件。如果未指定，则隐式为 0。
您可以在 Python 中将
event.terminal = True等属性赋给任何函数。- vectorized布尔值, 可选
是否可以以向量化方式调用 fun。默认为 False。
如果
vectorized为 False，fun 将始终以形状为(n,)的y调用，其中n = len(y0)。如果
vectorized为 True，fun 可能会以形状为(n, k)的y调用，其中k是整数。在这种情况下，fun 必须表现为fun(t, y)[:, i] == fun(t, y[:, i])（即返回数组的每一列是与y的一列对应的状态的时间导数）。设置
vectorized=True允许 'Radau' 和 'BDF' 方法更快地进行有限差分近似雅可比矩阵，但在某些情况下（例如len(y0)较小）会导致其他方法以及 'Radau' 和 'BDF' 的执行速度变慢。- args元组, 可选
要传递给用户定义函数的额外参数。如果给定，则将额外参数传递给所有用户定义函数。因此，例如，如果 fun 的签名为
fun(t, y, a, b, c)，则 jac（如果给定）和任何事件函数都必须具有相同的签名，并且 args 必须是长度为 3 的元组。- **options
传递给所选求解器的选项。下面列出了已实现的求解器可用的所有选项。
- first_step浮点数或 None, 可选
初始步长。默认为 None，这意味着算法应自行选择。
- max_step浮点数, 可选
允许的最大步长。默认为 np.inf，即步长不受限制，仅由求解器决定。
- rtol, atol浮点数或类数组, 可选
相对和绝对容差。求解器将局部误差估计保持在
atol + rtol * abs(y)以下。这里 rtol 控制相对精度（正确位数），而 atol 控制绝对精度（正确小数位数）。为了达到所需的 rtol，请将 atol 设置为小于rtol * abs(y)中可能预期的最小值的数字，以便 rtol 在允许误差中占主导地位。如果 atol 大于rtol * abs(y)，则不保证正确位数。相反，为了达到所需的 atol，请设置 rtol，使rtol * abs(y)始终小于 atol。如果 y 的分量具有不同的比例，则通过为 atol 传递形状为 (n,) 的类数组，为不同分量设置不同的 atol 值可能会有所帮助。默认值是 rtol 为 1e-3，atol 为 1e-6。- jac类数组, 稀疏矩阵, 可调用对象或 None, 可选
系统右侧关于 y 的雅可比矩阵，'Radau'、'BDF' 和 'LSODA' 方法需要。雅可比矩阵的形状为 (n, n)，其元素 (i, j) 等于
d f_i / d y_j。定义雅可比矩阵有三种方式：如果为类数组或稀疏矩阵，则假定雅可比矩阵为常数。'LSODA' 不支持。
如果为可调用对象，则假定雅可比矩阵同时依赖于 t 和 y；它将根据需要以
jac(t, y)的形式调用。如果使用args（请参阅args参数的文档），则必须传递额外参数。对于 'Radau' 和 'BDF' 方法，返回值可以是稀疏矩阵。如果为 None（默认），则雅可比矩阵将通过有限差分近似。
通常建议提供雅可比矩阵，而不是依赖于有限差分近似。
- jac_sparsity类数组, 稀疏矩阵或 None, 可选
定义雅可比矩阵的稀疏结构，用于有限差分近似。其形状必须为 (n, n)。如果 jac 不为 None，则忽略此参数。如果雅可比矩阵在每行中只有少数非零元素，提供稀疏结构将大大加快计算速度 [10]。零条目表示雅可比矩阵中相应的元素始终为零。如果为 None（默认），则假定雅可比矩阵是密集的。'LSODA' 不支持，请改用 lband 和 uband。
- lband, uband整数或 None, 可选
定义 'LSODA' 方法雅可比矩阵带宽的参数，即
jac[i, j] != 0 仅当 i - lband <= j <= i + uband。默认为 None。设置这些参数要求您的 jac 例程以打包格式返回雅可比矩阵：返回的数组必须有n列和uband + lband + 1行，其中写入雅可比矩阵的对角线。具体来说，jac_packed[uband + i - j , j] = jac[i, j]。相同的格式用于scipy.linalg.solve_banded（查看示例）。这些参数也可以与jac=None一起使用，以减少通过有限差分估计的雅可比矩阵元素的数量。- min_step浮点数, 可选
'LSODA' 方法允许的最小步长。默认情况下，min_step 为零。
- 返回:
- 具有以下定义的 Bunch 对象
- tndarray, 形状 (n_points,)
时间点。
- yndarray, 形状 (n, n_points)
解在 t 处的值。
- sol
OdeSolution或 None 找到的解作为
OdeSolution实例；如果 dense_output 设置为 False，则为 None。- t_eventsndarray 列表或 None
对于每种事件类型，包含检测到该类型事件的数组列表。如果 events 为 None，则为 None。
- y_eventsndarray 列表或 None
对于 t_events 的每个值，解的对应值。如果 events 为 None，则为 None。
- nfev整数
右侧的评估次数。
- njev整数
雅可比矩阵的评估次数。
- nlu整数
LU 分解的次数。
- status整数
算法终止原因
-1: 积分步骤失败。
0: 求解器成功到达 tspan 的末尾。
1: 发生终止事件。
- message字符串
终止原因的人类可读描述。
- success布尔值
如果求解器到达区间末尾或发生终止事件（
status >= 0），则为 True。
参考文献
[1]J. R. Dormand, P. J. Prince, “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.
[2]L. W. Shampine, “Some Practical Runge-Kutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
[3]P. Bogacki, L.F. Shampine, “A 3(2) Pair of Runge-Kutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
[4]E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”, Sec. IV.8.
[5][6]L. F. Shampine, M. W. Reichelt, “THE MATLAB ODE SUITE”, SIAM J. SCI. COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
[7]A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 55-64, 1983.
[8]L. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations”, SIAM Journal on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148, 1983.
[9][10]A. Curtis, M. J. D. Powell, and J. Reid, “On the estimation of sparse Jacobian matrices”, Journal of the Institute of Mathematics and its Applications, 13, pp. 117-120, 1974.
[11][13]E. Hairer, S. P. Norsett G. Wanner, “Solving Ordinary Differential Equations I: Nonstiff Problems”, Sec. II.
示例
显示自动选择时间点的基本指数衰减。
>>> import numpy as np >>> from scipy.integrate import solve_ivp >>> def exponential_decay(t, y): return -0.5 * y >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8]) >>> print(sol.t) [ 0. 0.11487653 1.26364188 3.06061781 4.81611105 6.57445806 8.33328988 10. ] >>> print(sol.y) [[2. 1.88836035 1.06327177 0.43319312 0.18017253 0.07483045 0.03107158 0.01350781] [4. 3.7767207 2.12654355 0.86638624 0.36034507 0.14966091 0.06214316 0.02701561] [8. 7.5534414 4.25308709 1.73277247 0.72069014 0.29932181 0.12428631 0.05403123]]
指定需要解的点。
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8], ... t_eval=[0, 1, 2, 4, 10]) >>> print(sol.t) [ 0 1 2 4 10] >>> print(sol.y) [[2. 1.21305369 0.73534021 0.27066736 0.01350938] [4. 2.42610739 1.47068043 0.54133472 0.02701876] [8. 4.85221478 2.94136085 1.08266944 0.05403753]]
炮弹向上发射，并在撞击时发生终止事件。事件的
terminal和direction字段通过猴子补丁函数应用。这里y[0]是位置，y[1]是速度。炮弹从位置 0 开始，速度为 +10。请注意，积分永远不会达到 t=100，因为事件是终止的。>>> def upward_cannon(t, y): return [y[1], -0.5] >>> def hit_ground(t, y): return y[0] >>> hit_ground.terminal = True >>> hit_ground.direction = -1 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground) >>> print(sol.t_events) [array([40.])] >>> print(sol.t) [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
使用 dense_output 和 events 找到炮弹轨迹顶点处的位置（即 100）。顶点未定义为终止，因此同时找到顶点和着地事件。t=20 处没有信息，因此使用 sol 属性评估解。通过设置
dense_output=True返回 sol 属性。或者，可以使用 y_events 属性访问事件发生时的解。>>> def apex(t, y): return y[1] >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], ... events=(hit_ground, apex), dense_output=True) >>> print(sol.t_events) [array([40.]), array([20.])] >>> print(sol.t) [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01] >>> print(sol.sol(sol.t_events[1][0])) [100. 0.] >>> print(sol.y_events) [array([[-5.68434189e-14, -1.00000000e+01]]), array([[1.00000000e+02, 1.77635684e-15]])]
作为带有额外参数的系统示例，我们将实现 Lotka-Volterra 方程 [12]。
>>> def lotkavolterra(t, z, a, b, c, d): ... x, y = z ... return [a*x - b*x*y, -c*y + d*x*y] ...
我们使用 args 参数传入参数值 a=1.5, b=1, c=3 和 d=1。
>>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1), ... dense_output=True)
计算密集解并绘制。
>>> t = np.linspace(0, 15, 300) >>> z = sol.sol(t) >>> import matplotlib.pyplot as plt >>> plt.plot(t, z.T) >>> plt.xlabel('t') >>> plt.legend(['x', 'y'], shadow=True) >>> plt.title('Lotka-Volterra System') >>> plt.show()
使用 solve_ivp 求解带有复数矩阵
A的微分方程y' = Ay的几个示例。>>> A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j], ... [0.25 + 0.58j, -0.2 + 0.14j, 0], ... [0, 0.2 + 0.4j, -0.1 + 0.97j]])
使用上述
A和 3x1 向量y求解 IVP>>> def deriv_vec(t, y): ... return A @ y >>> result = solve_ivp(deriv_vec, [0, 25], ... np.array([10 + 0j, 20 + 0j, 30 + 0j]), ... t_eval=np.linspace(0, 25, 101)) >>> print(result.y[:, 0]) [10.+0.j 20.+0.j 30.+0.j] >>> print(result.y[:, -1]) [18.46291039+45.25653651j 10.01569306+36.23293216j -4.98662741+80.07360388j]
使用上述
A和 3x3 矩阵y求解 IVP>>> def deriv_mat(t, y): ... return (A @ y.reshape(3, 3)).flatten() >>> y0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j], ... [5 + 0j, 6 + 0j, 7 + 0j], ... [9 + 0j, 34 + 0j, 78 + 0j]])
>>> result = solve_ivp(deriv_mat, [0, 25], y0.flatten(), ... t_eval=np.linspace(0, 25, 101)) >>> print(result.y[:, 0].reshape(3, 3)) [[ 2.+0.j 3.+0.j 4.+0.j] [ 5.+0.j 6.+0.j 7.+0.j] [ 9.+0.j 34.+0.j 78.+0.j]] >>> print(result.y[:, -1].reshape(3, 3)) [[ 5.67451179 +12.07938445j 17.2888073 +31.03278837j 37.83405768 +63.25138759j] [ 3.39949503 +11.82123994j 21.32530996 +44.88668871j 53.17531184+103.80400411j] [ -2.26105874 +22.19277664j -15.1255713 +70.19616341j -38.34616845+153.29039931j]]