使用 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] 区间内的任何点都不会被评估，所以没有什么可以使假设无效。

在第一个示例中情况并非如此，因为前两个组件确实发生了变化，因此迫使求解器仔细处理跳转点。