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)[源代码]¶
求解常微分方程组的初值问题。
此函数在给定初始值的情况下对常微分方程组进行数值积分：
dy / dt = f(t, y) y(t0) = y0
这里t是一维自变量(时间)，y(T)是N维矢值函数(状态)，N维矢值函数f(t，y)确定微分方程。目标是在给定初始值y(T0)=y0的情况下找到近似满足微分方程的y(T)。
一些求解器支持复数域中的积分，但请注意，对于刚性ODE求解器，右侧必须是复数可微的(满足Cauchy-Riemann方程 [11]) 。要解决复杂域中的问题，请使用复杂数据类型传递y0。另一个总是可用的选择是分别重写实部和虚部的问题。
- 参数
- fun可调用
系统的右侧。调用签名为
fun(t, y)。这里 t 是标量，ndarray有两个选项 y ：它可以是形状(n，)；然后 fun 必须返回形状为(n，)的ARRAY_LIKE。或者，它可以具有形状(n，k)；然后 fun 必须返回形状为(n，k)的array_like，即每列对应于 y 。在这两个选项之间的选择由以下因素决定 vectorized 参数(见下文)。矢量化实现允许通过有限差分更快地逼近雅可比(对于刚性解算器是必需的)。- t_span浮点数的2元组
积分间隔(t0，tf)。求解器从t=t0开始并进行积分，直到达到t=tf。
- y0类似数组，形状(n，)
初始状态。对于复杂领域中的问题，请通过 y0 具有复杂数据类型(即使初始值是纯实数)。
- 方法 ：字符串或
OdeSolver，可选字符串或 要使用的集成方法：
‘RK45’(默认)：5阶显式Runge-Kutta方法(4) [1]. 假设四阶方法的精度，误差是可以控制的，但采取的步骤是使用五阶精度公式(进行局部外推)。四次插值多项式用于密集输出。 [2]. 可应用于复杂领域。
“RK23”：3(2)阶显式Runge-Kutta方法 [3]. 在假设二阶方法精度的前提下，误差是可控制的，但采取的步骤是使用三阶精度公式(进行局部外推)。三次Hermite多项式用于密集输出。可应用于复杂领域。
“DOP853”：8阶显式Runge-Kutta方法 [13]. 最初用Fortran编写的“DOP853”算法的Python实现 [14]. 对于密集输出，采用精确到7阶的7次插值多项式。可应用于复杂领域。
“Radau”：5阶Radau IIA族的隐式Runge-Kutta方法 [4]. 误差由三阶精确的嵌入公式控制。对于密集输出，采用满足配置条件的三次多项式。
“BDF”：基于导数逼近向后微分公式的隐式多步变阶(1-5)方法 [5]. 该实现遵循中描述的实现 [6]. 采用准等步长格式，并通过NDF修正提高了精度。可应用于复杂领域。
‘LSODA’：具有自动刚度检测和切换功能的ADAMS/BDF方法 [7], [8]. 这是ODEPACK提供的Fortran求解器的包装。
对于非刚性问题应该使用显式Runge-Kutta方法(‘RK23’，‘RK45’，‘DOP853’)，对于刚性问题应该使用隐式方法(‘Radau’，‘BDF’) [9]. 在龙格-库塔(Runge-Kutta)方法中，推荐使用‘DOP853’求解精度较高(数值较低)的方法 rtol 和 atol )。
如果不确定，请先尝试运行“RK45”。如果它进行了异常多次的迭代、发散或失败，那么您的问题很可能很棘手，您应该使用“Radau”或“bdf”。“LSODA”也可以是一个很好的通用选择，但由于它包装了旧的Fortran代码，使用起来可能不太方便。
还可以传递从
OdeSolver它实现了求解器。- t_evalARRAY_LIKE或NONE，可选
必须对存储计算结果的时间进行排序，并将其置于 t_span 。如果为None(默认)，则使用解算器选择的点。
- dense_output布尔值，可选
是否计算连续解。默认值为False。
- events可调用或可调用列表，可选
要跟踪的事件。如果为None(默认值)，则不会跟踪任何事件。每个事件都发生在时间和状态的连续函数的零点。每个函数都必须有签名
event(t, y)并返回一个浮点数。求解器将找到一个精确值 t 在哪一天event(t, y(t)) = 0使用寻根算法。默认情况下，将找到全零。求解器会在每个步骤中查找符号变化，因此如果在一个步骤中出现多个过零事件，则可能会错过事件。另外，每个 event 函数可能具有以下属性：- 终端：布尔型，可选
如果发生此事件，是否终止集成。如果未赋值，则隐式为false。
- 方向：浮点，可选
过零的方向。如果 direction 是积极的， event 仅在从负数变为正数时触发，反之亦然 direction 是阴性的。如果为0，则任一方向都将触发事件。如果未分配，则隐式为0。
您可以指定如下属性
event.terminal = True添加到Python中的任何函数。- vectorized布尔值，可选
是否 fun 是以矢量化的方式实现的。默认值为False。
- args元组，可选
要传递给用户定义函数的其他参数。如果给定，则将附加参数传递给所有用户定义函数。所以，举个例子，如果 fun 有签名的
fun(t, y, a, b, c)，那么 jac (如果指定)，并且任何事件函数都必须具有相同的签名，并且 args 必须是长度为3的元组。- options
传递给选定求解器的选项。下面列出了已实现的解算器的所有可用选项。
- first_step浮动或无，可选
初始步长。默认值为 None 这意味着算法应该选择。
- max_step浮动，可选
允许的最大步长。默认值为np.inf，即步长不受限制，仅由解算器确定。
- RTOL，ATOLFLOAT或ARRAY_LIKE，可选
相对公差和绝对公差。求解器使局部误差估计值小于
atol + rtol * abs(y)。这里 rtol 控制相对精度(正确的位数)。但是如果一个组件 y 大约低于 atol ，误差只需要落在相同的 atol 阈值，并且不保证正确的位数。如果y的分量具有不同的刻度，则设置不同的刻度可能是有益的 atol 通过为传递具有形状(n，)的array_like来获取不同组件的值 atol 。的默认值为1e-3 rtol 1e-6用于 atol 。- jacARRAY_LIKE、SPARSE_MATRIX、Callable或None，可选
系统右侧相对于y的雅可比矩阵，这是‘Radau’、‘BDF’和‘LSODA’方法所要求的。雅可比矩阵具有形状(n，n)，其元素(i，j)等于
d f_i / d y_j。定义雅可比有三种方法：如果为ARRAY_LIKE或SPARSE_MATRATE，则假定雅可比为常量。‘LSODA’不支持。
如果是可调用的，则假定雅可比依赖于t和y；它将被称为
jac(t, y)，根据需要。对于“Radau”和“bdf”方法，返回值可能是稀疏矩阵。如果为None(默认)，则雅可比将由有限差分近似。
通常建议提供雅可比，而不是依赖有限差分近似。
- jac_sparsityARRAY_LIKE、稀疏矩阵或无，可选
定义了用于有限差分近似的雅可比矩阵的稀疏结构。其形状必须为(n，n)。如果出现以下情况，则忽略此参数 jac 不是 None 。如果雅可比函数中只有几个非零元素 each 行，提供稀疏结构将大大加快计算速度 [10]. 零条目表示雅可比中的相应元素始终为零。如果为None(默认)，则假定雅可比是密集的。“LSODA”不支持，请参见 lband 和 uband 取而代之的是。
- 带区，子带整型或无型，可选
定义‘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_step浮动，可选
“LSODA”方法允许的最小步长。默认情况下， min_step 是零。
- 退货
- 捆绑定义了以下字段的对象：
- tndarray，形状(n_point，)
时间点。
- yndarray，形状(n，n_点)
位于以下位置的解的值 t 。
- sol ：
OdeSolution或无OdeSolution或无 找到的解决方案为
OdeSolution实例；如果 dense_output 设置为false。- t_eventsndarray或None列表
包含每个事件类型的数组列表，在这些数组中检测到该类型的事件。如果有，则无 events 什么都没有。
- y_eventsndarray或None列表
对于每个值 t_events ，则为解的对应值。如果有，则无 events 什么都没有。
- nfev集成
右侧的求值次数。
- njev集成
雅可比的求值次数。
- nlu集成
逻辑单元分解数。
- status集成
算法终止原因：
-1：集成步骤失败。
0：求解器已成功到达 tspan 。
1：发生终止事件。
- message字符串
终止原因的人类可读描述。
- success布尔尔
如果求解程序达到间隔结束或发生终止事件，则为True (
status >= 0)。
参考文献
- 1
J.R.Dormand，P.J.Prince，“一族嵌入式Runge-Kutta公式”，“计算与应用数学杂志”，第6卷，第1期，第19-26页，1980年。
- 2
夏姆平，“一些实用的龙格-库塔公式”，“计算数学”，第46卷，第173期，第135-150页，1986。
- 3
P.Bogacki，L.F.Shampine，“A 3(2)对Runge-Kutta公式”，Appl.数学课。让我们来吧。第2卷，第4期，第321-325页，1989年。
- 4
海尔，G.Wanner，“解常微分方程II：刚性和微分-代数问题”，SEC。IV.8.
- 5
Backward Differentiation Formula 在维基百科上。
- 6
书名/作者：The MATLAB Ode Suite(MATLAB颂歌组曲)，SIAM J.SCI。计算，第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
Stiff equation 在维基百科上。
- 10
A.Curtis，M.J.D.Powell和J.Reid，“关于稀疏雅可比矩阵的估计”，“数学研究所及其应用期刊”，第13期，第117-120页，1974年。
- 11
Cauchy-Riemann equations 在维基百科上。
- 12
Lotka-Volterra equations 在维基百科上。
- 13
海尔，S.P.Norsett G.Wanner，“解常微分方程I：非刚性问题”，美国科学出版社。二、二、
- 14
Page with original Fortran code of DOP853 _.
示例
基本指数衰减，显示自动选择的时间点。
>>> 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]就是速度。弹丸以+10的速度从位置0开始。请注意，积分永远不会达到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]])]
作为具有附加参数的系统的示例，我们将实现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] ...
我们将参数值a=1.5、b=1、c=3和d=1与 args 论点。
>>> 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()
