使用 scipy.integrate 的 solve_ivp 的意外行为

Question

我想用一个开关函数来集成一个简单的ODE。切换函数将确保对于某些 t 值梯度为零，并使梯度对于任何其他 t 值成为数值常数。

我可以通过使用虚拟状态获得所需的结果，但是当我对单个状态重复计算时，solve_ivp 的输出是不同的。我想明白为什么会这样。

这是重现结果的代码：

#==============================================================================
# import modules
#==============================================================================

from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
import numpy as np

#==============================================================================
# user defined functions
#==============================================================================
    
def ode_multi(t, y):
    
   dy1_dt = 0 if t < 3 else 1
   dy2_dt = 0 if t < 4 else 1
   dy3_dt = dy1_dt - dy2_dt
    
   return [dy1_dt, dy2_dt, dy3_dt]

def ode_single(t, y):
    
   dy1_dt = 0 if t < 3 else 1
   dy2_dt = 0 if t < 4 else 1
   dy3_dt = dy1_dt - dy2_dt
    
   return [dy3_dt]

#==============================================================================
# solve ode
#==============================================================================

sol_multi = solve_ivp(ode_multi, [0,6], [0,0,0], method='LSODA', t_eval = np.linspace(0,6,100))
sol_single = solve_ivp(ode_single, [0,6], [0], method='LSODA', t_eval = np.linspace(0,6,100))

plt.figure(1)
plt.plot(sol_multi.t, sol_multi.y[0], label='dy1_dt')
plt.plot(sol_multi.t, sol_multi.y[1], label='dy2_dt')
plt.plot(sol_multi.t, sol_multi.y[2], label='dy3_dt')
plt.legend(); plt.grid()

plt.figure(2)
plt.plot(sol_single.t, sol_single.y[0], label='dy3_dt')
plt.legend(); plt.grid()

Ode_multi 输出：

Ode_single输出

我尝试使用 numpy.heaviside 而不是切换函数的 if 语句来确保梯度是可微的。输出是相同的，所以我在代码示例中使用了 if 语句来简化解释。

感谢您的意见！

Answer 1

将步长限制为“奇异”事件之间差异的一半，max_step=0.5 应该适用于这种情况。降低误差容限可以产生类似的效果。

好像积分器跳过了导数不为零的区间。就求解器而言，当内部步长适应导数函数时，就会发生这种情况。现在在初始零段中，求解器假设导数函数是零函数，并相应地增加步长。幸运的是，区间 [3,4] 中的任何一点都没有被评估，所以没有任何东西使假设无效。

第一个示例不是这种情况，因为前两个组件确实发生了变化，因此迫使求解器仔细处理跳跃点。

使用 scipy.integrate 的 solve_ivp 的意外行为

Unexpected behavior using scipy.integrate's solve_ivp

python

scipy

ode