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)[源代码]#
求解常微分方程组的初值问题。
此函数以数值方式积分给定的初值的常微分方程组。
dy / dt = f(t, y) y(t0) = y0
其中 t 是 1 维独立变量（时间），y(t) 是 N 维向量值函数（状态），N 维向量值函数 f(t, y) 确定微分方程。目标是找到近似满足微分方程的 y(t)，给定初值 y(t0)=y0。
一些求解器支持复数域中的积分，但请注意，对于刚性 ODE 求解器，右侧必须是复数可微的（满足柯西-黎曼方程 [11]）。要解决复数域中的问题，请传递具有复数数据类型的 y0。另一个始终可用的选项是将您的问题分别重写为实部和虚部。
- 参数:
- funcallable
系统的右侧：时间
t时刻状态y的时间导数。调用签名是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 都必须是浮点数或可以由浮点数转换函数解释的值。
- y0array_like, shape (n,)
初始状态。对于复数域中的问题，请传递具有复数数据类型的 y0（即使初始值是纯实数）。
- methodstring 或
OdeSolver，可选 要使用的积分方法
‘RK45’ (默认)：5(4) 阶显式龙格-库塔方法 [1]。误差控制假设四阶方法的精度，但使用五阶精确公式进行步进（进行局部外推）。四次插值多项式用于密集输出 [2]。可以在复数域中应用。
‘RK23’：3(2) 阶显式龙格-库塔方法 [3]。误差控制假设二阶方法的精度，但使用三阶精确公式进行步进（进行局部外推）。三次埃尔米特多项式用于密集输出。可以在复数域中应用。
‘DOP853’：8 阶显式龙格-库塔方法 [13]。最初用 Fortran 编写的“DOP853”算法的 Python 实现 [14]。7 阶插值多项式，精确到 7 阶，用于密集输出。可以在复数域中应用。
‘Radau’：Radau IIA 系列的 5 阶隐式龙格-库塔方法 [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_evalarray_like 或 None，可选
存储计算解的时间，必须排序并在 t_span 内。如果为 None（默认），则使用求解器选择的点。
- dense_outputbool，可选
是否计算连续解。默认值为 False。
- eventscallable，或 callable 列表，可选
要跟踪的事件。如果为 None（默认），则不会跟踪任何事件。每个事件都发生在时间和状态的连续函数的零点。每个函数都必须具有签名
event(t, y)，如果使用了args，则必须传递其他参数（请参阅args参数的文档）。每个函数都必须返回一个浮点数。求解器将使用寻根算法找到event(t, y(t)) = 0的精确值 t。默认情况下，将找到所有零点。求解器会在每个步骤中查找符号变化，因此如果一个步骤中发生多个零交叉，则可能会错过事件。此外，每个 event 函数可能具有以下属性- terminal: bool 或 int，可选
当为布尔值时，是否在此事件发生时终止积分。当为整数时，在指定次数的此事件发生后终止。如果未分配，则默认为 False。
- direction: float，可选
零交叉的方向。如果 direction 为正，则 event 仅在从负变为正时触发，如果 direction 为负，则反之亦然。如果为 0，则任何方向都会触发事件。如果未分配，则默认为 0。
您可以在 Python 中将诸如
event.terminal = True之类的属性分配给任何函数。- vectorizedbool，可选
是否可以以向量化方式调用 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’ 的雅可比矩阵的有限差分逼近，但对于其他方法以及在某些情况下（例如，小len(y0)），对于 ‘Radau’ 和 ‘BDF’ 将导致执行速度变慢。- argstuple，可选
传递给用户自定义函数的附加参数。如果给定，则附加参数将传递给所有用户自定义函数。例如，如果 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、稀疏矩阵或 None，可选
定义用于有限差分近似的雅可比矩阵的稀疏结构。其形状必须为 (n, n)。如果 jac 不是 None，则忽略此参数。如果雅可比矩阵的每行只有几个非零元素，提供稀疏结构将大大加快计算速度 [10]。零条目表示雅可比矩阵中的相应元素始终为零。如果为 None（默认值），则假定雅可比矩阵为稠密矩阵。'LSODA' 不支持，请改用 lband 和 uband。
- lband, ubandint 或 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_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：发生终止事件。
- message字符串
终止原因的人类可读描述。
- successbool
如果求解器到达区间末尾或发生终止事件（
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 属性来评估解。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]])]
作为带有附加参数的系统示例，我们将实现 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 求解微分方程
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和 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]]