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(t0)=y0 的情况下，近似找到满足微分方程的 y(t)。
一些求解器支持在复数域中的积分，但要注意，对于刚性 ODE 求解器，右端必须是复可微的（满足柯西-黎曼方程 [11]）。要解决复数域中的问题，请传递具有复数数据类型的 y0。始终可用的另一个选项是分别重写实部和虚部的求解问题。
- 参数:
- fun可调用对象
系统的右端：状态
y在时间t处的时变导数。调用签名为fun(t, y)，其中t是一个标量，y是一个具有len(y) = len(y0)的 ndarray。如果使用了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]。最初用 Fortran 语言编写的“DOP853”算法的 Python 实现 [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参数的文档）。每个函数都必须返回一个浮点数。求解器将使用一个求根算法找到 t 的精确值，使得event(t, y(t)) = 0。默认情况下，将找到所有零点。求解器在每个步骤中查找符号变化，因此，如果在一步内发生多个零交叉，则可能会错过事件。此外，每个 event 函数可能具有以下属性- terminal: 布尔值或整数，可选
当为布尔值时，如果发生此事件，是否终止积分。当为整数时，在指定次数的此事件发生后终止。如果未分配，则隐式为 False。
- direction: 浮点数，可选
零交叉的方向。如果 direction 为正，则 event 将仅在从负到正时触发，反之亦然，如果 direction 为负。如果为 0，则任一方向都会触发事件。如果未分配，则隐式为 0。
你可以将诸如
event.terminal = True之类的属性分配给 Python 中的任何函数。- vectorized布尔值，可选
是否可以以向量化方式调用 fun。默认为 False。
如果
vectorized为 False，则始终会使用形状为(n,)的y调用 fun，其中n = len(y0)。如果
vectorized为 True，则可以使用形状为(n, k)的y调用 fun，其中k是一个整数。在这种情况下，fun 必须表现得像fun(t, y)[:, i] == fun(t, y[:, i])（即，返回数组的每一列都是与y的一列相对应的状态的时变导数）。设置
vectorized=True允许更快地使用方法 ‘Radau’ 和 ‘BDF’ 对雅可比矩阵进行有限差分近似，但在某些情况下会导致其他方法以及 ‘Radau’ 和 ‘BDF’ 的执行速度变慢（例如，len(y0)较小）。- args元组，可选
要传递给用户定义函数的额外参数。如果给出，则将额外参数传递给所有用户定义的函数。例如，如果 fun 的签名为
fun(t, y, a, b, c)，则 jac（如果给出）和任何事件函数必须具有相同的签名，并且 args 必须是一个长度为 3 的元组。- **options**
传递给所选求解器的选项。以下列出了已实现求解器的所有可用选项。
- first_stepfloat 或 None，可选
初始步长。默认为 None，这意味着算法应自行选择。
- max_stepfloat，可选
允许的最大步长。默认为 np.inf，即步长不受限制，仅由求解器决定。
- rtol，atolfloat 或 array_like，可选
相对和绝对容差。求解器将局部误差估计保持在小于
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,) 的 array_like，可能有利于为不同的分量设置不同的 atol 值。默认值为 rtol 为 1e-3，atol 为 1e-6。- jacarray_like，sparse_matrix，callable 或 None，可选
系统右侧关于 y 的雅可比矩阵，‘Radau’、‘BDF’ 和 ‘LSODA’ 方法需要它。雅可比矩阵的形状为 (n, n)，其元素 (i, j) 等于
d f_i / d y_j。有三种方法可以定义雅可比矩阵如果为 array_like 或 sparse_matrix，则假定雅可比矩阵是常数。‘LSODA’ 不支持。
如果为 callable，则假定雅可比矩阵依赖于 t 和 y；它将在必要时被调用为
jac(t, y)。如果使用args（请参阅args参数的文档），则必须传递其他参数。对于 ‘Radau’ 和 ‘BDF’ 方法，返回值可能是一个稀疏矩阵。如果为 None（默认值），则雅可比矩阵将通过有限差分近似。
通常建议提供雅可比矩阵，而不是依赖有限差分近似。
- jac_sparsityarray_like，sparse matrix 或 None，可选
定义有限差分近似雅可比矩阵的稀疏性结构。其形状必须为 (n, n)。如果 jac 不是 None，则忽略此参数。如果雅可比矩阵的每一行中只有几个非零元素，则提供稀疏性结构将极大地提高计算速度 [10]。零条目表示雅可比矩阵中的对应元素始终为零。如果为 None（默认值），则假定雅可比矩阵是密集的。‘LSODA’ 不支持，请改用 lband 和 uband。
- lband，ubandint 或 None，可选
定义 ‘LSODA’ 方法的雅可比矩阵带宽的参数，即
jac[i, j] != 0 only for 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_stepfloat，可选
‘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。
- nfevint
右侧评估次数。
- njevint
雅可比矩阵评估次数。
- nluint
LU 分解次数。
- statusint
算法终止的原因
-1: 积分步长失败。
0: 求解器成功到达 tspan 的末尾。
1: 发生终止事件。
- messagestring
终止原因的人类可读描述。
- successbool
如果求解器到达区间末尾或发生终止事件（
status >= 0），则为 True。
参考文献
[1]J. R. Dormand、P. J. Prince，“嵌入式龙格-库塔公式族”，计算与应用数学杂志，第 6 卷，第 1 期，第 19-26 页，1980 年。
[2]L. W. Shampine，“一些实用的龙格-库塔公式”，计算数学，第 46 卷，第 173 期，第 135-150 页，1986 年。
[3]P. Bogacki、L.F. Shampine，“一对 3(2) 龙格-库塔公式”，应用数学通讯，第 2 卷，第 4 期。第 321-325 页，1989 年。
[4]E. Hairer、G. Wanner，“常微分方程求解 II：刚性和微分代数问题”，第 IV.8 节。
[5][6]L. F. Shampine、M. W. Reichelt，“MATLAB ODE 套件”，SIAM 科学计算杂志，第 18 卷，第 1 期，第 1-22 页，1997 年 1 月。
[7]A. C. Hindmarsh，“ODEPACK，一个系统化的 ODE 求解器集合”，科学计算 IMACS 交易，第 1 卷，第 55-64 页，1983 年。
[8]L. Petzold，“为解决刚性和非刚性常微分方程组自动选择方法”，SIAM 科学与统计计算杂志，第 4 卷，第 1 期，第 136-148 页，1983 年。
[9][10]A. Curtis、M. J. D. Powell 和 J. Reid，“关于稀疏雅可比矩阵估计”，数学及其应用研究所杂志，13，第 117-120 页，1974 年。
[11][13]E. Hairer、S. P. Norsett G. Wanner，“常微分方程求解 I：非刚性问题”，第 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。顶点没有定义为终止点，因此找到顶点和 hit_ground。在 t=20 时没有信息，因此使用 sol 属性来评估解。sol 属性是通过设置
dense_output=True来返回的。或者，可以使用 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]])]
作为具有附加参数的系统的示例，我们将实现洛特卡-沃尔泰拉方程 [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 解决微分方程
y' = Ay的几个示例，其中A是复矩阵。>>> 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和y作为 3x1 向量解决 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和y作为 3x3 矩阵解决 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]]