Python SciPy integrate.solve_ivp用法及代码示例

本文简要介绍 python 语言中 scipy.integrate.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)#

求解 ODE 系统的初始值问题。

该函数对给定初始值的常微分方程组进行数值积分：

dy / dt = f(t, y)
y(t0) = y0

这里 t 是一维自变量(时间)，y(t) 是N-D vector-valued 函数(状态)，N-D vector-valued 函数 f(t, y) 确定微分方程。目标是在给定初始值 y(t0)=y0 的情况下找到近似满足微分方程的 y(t)。

一些求解器支持复数域中的积分，但请注意，对于刚性 ODE 求解器，右侧必须是 complex-differentiable(满足 Cauchy-Riemann 方程 [11])。要解决复杂域中的问题，请传递具有复杂数据类型的 y0。另一个始终可用的选择是分别重写实部和虚部的问题。

参数 ：：

fun： 可调用的

系统右侧：状态的时间导数y在某个时间t。调用签名是fun(t, y)，其中t是一个标量并且y是一个 ndarraylen(y) = len(y0)。需要传递额外的参数，如果args使用(参见文档args争论)。fun必须返回与以下形状相同的数组y.看矢量化了解更多信息。

t_span： 2元序列

积分间隔(t0，tf)。求解器从 t=t0 开始进行积分，直到达到 t=tf。 t0 和 tf 都必须是浮点数或可由浮点转换函数解释的值。

y0： 数组, 形状 (n,)

初始状态。对于复杂域中的问题，使用复杂数据类型传递 y0(即使初始值是纯实数)。

method： 字符串或 OdeSolver ，可选

使用的积分方法：

• ‘RK45’ (default): Explicit Runge-Kutta method of order 5(4) [1]. The error is controlled assuming accuracy of the fourth-order method, but steps are taken using the fifth-order accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [2]. Can be applied in the complex domain.

• ‘RK23’: Explicit Runge-Kutta method of order 3(2) [3]. The error is controlled assuming accuracy of the second-order method, but steps are taken using the third-order accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.

• ‘DOP853’: Explicit Runge-Kutta method of order 8 [13]. Python implementation of the “DOP853” algorithm originally written in Fortran [14]. A 7-th order interpolation polynomial accurate to 7-th order is used for the dense output. Can be applied in the complex domain.

• ‘Radau’: Implicit Runge-Kutta method of the Radau IIA family of order 5 [4]. The error is controlled with a third-order accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output.

• ‘BDF’: Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation [5]. The implementation follows the one described in [6]. A quasi-constant step scheme is used and accuracy is enhanced using the NDF modification. Can be applied in the complex domain.

• ‘LSODA’: Adams/BDF method with automatic stiffness detection and switching [7], [8]. This is a wrapper of the Fortran solver from ODEPACK.

显式 Runge-Kutta 方法(‘RK23’、‘RK45’、‘DOP853’)应用于非刚性问题，隐式方法(‘Radau’、‘BDF’)应用于刚性问题[9].在Runge-Kutta 方法中，推荐使用‘DOP853’进行高精度求解(低值rol和环礁)。

如果不确定，请先尝试运行“RK45”。如果它进行异常多次迭代、发散或失败，您的问题可能会很僵硬，您应该使用“Radau”或“BDF”。 “LSODA”也可以是一个很好的通用选择，但使用它可能不太方便，因为它包装了旧的 Fortran 代码。

您还可以传递从实现求解器的 OdeSolver 派生的任意类。

t_eval： 数组 或 None，可选

存储计算解决方案的时间必须排序并位于 t_span 内。如果无(默认)，则使用求解器选择的点。

dense_output： 布尔型，可选

是否计算连续解。默认为假。

events： 可调用的，或可调用的列表，可选

要跟踪的事件。如果 None(默认)，则不会跟踪任何事件。每个事件都发生在时间和状态的连续函数的零点处。每个函数都必须有签名event(t, y)如果必须传递附加参数args使用(参见文档args争论)。每个函数必须返回一个浮点数。求解器将找到一个准确的值t在哪个event(t, y(t)) = 0使用root-finding 算法。默认情况下，将找到所有零。求解器会在每个步骤中寻找符号变化，因此如果在一个步骤中发生多个过零，则可能会错过事件。另外每个事件函数可能具有以下属性：

terminal: bool, optional

Whether to terminate integration if this event occurs. Implicitly False if not assigned.

direction: float, optional

Direction of a zero crossing. If direction is positive, event will only trigger when going from negative to positive, and vice versa if direction is negative. If 0, then either direction will trigger event. Implicitly 0 if not assigned.

您可以将 event.terminal = True 等属性分配给 Python 中的任何函数。

vectorized： 布尔型，可选

是否可以以矢量化方式调用乐趣。默认值为 False。

如果vectorized是假的，乐趣总是会被调用y形状的(n,)，其中n = len(y0).

如果vectorized是真的，乐趣可以调用y形状的(n, k)，其中k是一个整数。在这种情况下，乐趣必须表现得这样fun(t, y)[:, i] == fun(t, y[:, i])(即返回数组的每一列都是与一列对应的状态的时间导数y)。

设置vectorized=True允许通过方法“Radau”和“BDF”更快地进行雅可比行列式的有限差分近似，但在某些情况下会导致其他方法以及“Radau”和“BDF”的执行速度变慢(例如小len(y0) )。

args： 元组，可选

传递给用户定义函数的附加参数。如果给定，附加参数将传递给所有用户定义的函数。因此，例如，如果乐趣有签名fun(t, y, a, b, c)， 然后江淮(如果给定)并且任何事件函数必须具有相同的签名，并且参数必须是长度为 3 的元组。

**options：

传递给选定求解器的选项。下面列出了已实现求解器的所有可用选项。

first_step： 浮点数或无，可选

初始步长。默认为无，这意味着算法应该选择。

max_step： 浮点数，可选

最大允许步长。默认值为 np.inf，即步长不受限制，仅由求解器确定。

rtol, atol： float 或 数组，可选

相对和绝对公差。求解器保持局部误差估计小于atol + rtol * abs(y).这里rol控制相对精度(正确位数)，而环礁控制绝对精度(正确的小数位数)。达到预期rol， 放环礁小于可以预期的最小值rtol * abs(y)以便rol支配允许误差。如果环礁大于rtol * abs(y)不保证正确位数。反之，达到预期环礁setrol这样rtol * abs(y)总是小于环礁.如果 y 的分量具有不同的尺度，则设置不同的尺度可能是有益的环礁通过将 数组 与 shape (n,) 传递给不同组件的值环礁.默认值为 1e-3rol和 1e-6 为环礁.

jac： 数组，sparse_matrix，可调用或无，可选

“Radau”、“BDF”和“LSODA”方法所需的系统右侧关于 y 的雅可比矩阵。雅可比矩阵的形状为 (n, n)，其元素 (i, j) 等于 d f_i / d y_j 。定义雅可比行列式有以下三种方法：

• If array_like or sparse_matrix, the Jacobian is assumed to be constant. Not supported by ‘LSODA’.

• If callable, the Jacobian is assumed to depend on both t and y; it will be called as jac(t, y), as necessary. Additional arguments have to be passed if args is used (see documentation of args argument). For ‘Radau’ and ‘BDF’ methods, the return value might be a sparse matrix.

• If None (default), the Jacobian will be approximated by finite differences.

通常建议提供雅可比行列式而不是依赖有限差分近似。

jac_sparsity： 数组，稀疏矩阵或无，可选

为有限差分近似定义雅可比矩阵的稀疏结构。它的形状必须是 (n, n)。如果江淮不是None.如果雅可比矩阵中只有少数非零元素每个行，提供稀疏结构将大大加快计算速度[10].零条目意味着雅可比矩阵中的相应元素始终为零。如果无(默认)，则假定雅可比是密集的。 “LSODA”不支持，请参阅带和乌班反而。

lband, uband： int 或无，可选

定义“LSODA”方法的雅可比带宽的参数，即 jac[i, j] != 0 only for i - lband <= j <= i + uband 。默认为无。设置这些要求您的 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 对象：

t： ndarray，形状(n_points，)

时间点。

y： ndarray，形状(n，n_points)

t 处解的值。

sol： OdeSolution 或无

找到解决方案为OdeSolution实例;没有如果dense_output被设置为假。

t_events： ndarray 或无列表

包含每个事件类型的数组列表，在该列表中检测到该类型事件的事件。如果事件为无，则无。

y_events： ndarray 或无列表

对于 t_events 的每个值，解的对应值。如果事件为无，则无。

nfev： int

右侧的评估次数。

njev： int

雅可比的评估次数。

nlu： int

LU 分解的数量。

status： int

算法终止原因：

• -1: Integration step failed.

• 0: The solver successfully reached the end of tspan.

• 1: A termination event occurred.

message： string

终止原因的人类可读说明。

success： bool

如果求解器达到间隔结束或发生终止事件 (status >= 0)，则为真。

参考：

[1]

J. R. Dormand, P. J. Prince，“一系列嵌入式Runge-Kutta 公式”，计算与应用数学杂志，卷。 6，第 1 期，第 19-26 页，1980 年。

[2]

L. W. Shampine，“一些实用的Runge-Kutta 公式”，计算数学，卷。 46，第 173 期，第 135-150 页，1986 年。

[3]

P. Bogacki, L.F. Shampine，“A 3(2) Pair of Runge-Kutta Formulas”，Appl。数学。莱特。卷。 2，第 4 期。第 321-325 页，1989 年。

[4]

E. Hairer, G. Wanner，“求解常微分方程 II：刚性和 Differential-Algebraic 问题”，秒。 IV.8.

[5]

Backward Differentiation Formula 在维基百科上。

[6]

L. F. Shampine，M. W. Reichelt，“MATLAB ODE 套件”，SIAM J. SCI。计算，卷。 18，第 1 期，第 1-22 页，1997 年 1 月。

[7]

A. C. Hindmarsh，“ODEPACK，ODE 求解器的系统化集合”，IMACS Transactions on Scientific Computation，第 1 卷，第 55-64 页，1983 年。

[8]

L. Petzold，“自动选择求解常微分方程的刚性和非刚性系统的方法”，SIAM 科学与统计计算杂志，卷。 4，第 1 期，第 136-148 页，1983 年。

[9]

Stiff equation 在维基百科上。

[10]

A. Curtis、M. J. D. Powell 和 J. Reid，“关于稀疏雅可比矩阵的估计”，数学研究所及其应用杂志，13，第 117-120 页，1974 年。

[11]

Cauchy-Riemann equations 在维基百科上。

[12]

Lotka-Volterra equations 在维基百科上。

[13]

E. Hairer, S. P. Norsett G. Wanner，“求解常微分方程 I：非刚性问题”，秒。二、

[14]

Page with original Fortran code of DOP853 。

例子：

显示自动选择的时间点的基本 index 衰减。

>>> 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和事件找到炮弹轨迹顶点的位置，即 100。 Apex 未定义为终端，因此找到了 apex 和 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]])]

作为具有附加参数的系统的示例，我们将实现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 和 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 求解 IVP，并将 y 作为 3x3 矩阵：

>>> 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]]