scipy的微分方程组求解包solve_ivp使用详解

官方文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

基础：

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)
'''
部分输入参数：
fun：函数

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

y0：初值


'''

进阶模块：

method：
使用的积分方法：

'RK45'（默认）：5(4) 阶显式龙格-库塔方法[1]。假设四阶方法的精度可以控制误差，但使用五阶精确公式（完成局部外推）采取步骤。四次插值多项式用于密集输出[2]。可应用于复杂领域。

'RK23'：3(2) 阶显式龙格-库塔方法[3]。假设二阶方法的精度可以控制误差，但使用三阶精确公式（完成局部外推）采取步骤。三次 Hermite 多项式用于密集输出。可应用于复杂领域。

“DOP853”：8 阶显式龙格-库塔方法[13]。“DOP853”算法的 Python 实现最初是用 Fortran 编写的[14]。精确到七阶的七阶插值多项式用于密集输出。可应用于复杂领域。

“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_eval 

在(t0, tf)时间段中，存在时间发生时间点。

events

使用需要检测的函数。在ode求解过程中，每一步都会执行events对应的函数，函数返回的是需要进行判断该变量情况的数值：比如官方案例，加农炮发射，存在两个变量，高度和速度。因此构建ode方程upward_cannon。假设此时需要在上升阶段高度为2000时，设置新的速度初值。代码改写如下：

原代码

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]

改写后：

new_v = 50
def upward_cannon(t, y): return [y[1], -0.5]
def hit_ground(t, y): 
    if (y[0] - 2000)>=0:
        y[1] = new_v
    return y[0] - 2000
hit_ground.terminal = True
hit_ground.direction = -1
sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)

其中hit_ground是每一步都检测的函数，它返回的数值需要结合terminal&direction来看。direction为-1，就表示y[0]-2000从正值变负，事件发生触发，terminal表示当事件触发了，是否需要结束ode求解。