scipy: BUG: solve_ivp inaccurate for piecewise-smooth system

Describe your issue.

When solving a simple piecewise-smooth system using the default configuration for solve_ivp (specifically rtol=1e-3 and atol=1e-6), the solution is significantly incorrect (on the order of 1e-2). I’ve compared the solution with a manual rk4 solver and as expected the solution quickly converges to the same oscillation, regardless of initial condition. Whereas solve_ivp produces two very different solutions even when the initial conditions differ by 1e-7, and both of these solutions are completely wrong.

The output from solve_ivp can be made accurate by decreasing the tolerances, but I would have thought that the default values should be enough to get a reasonably accurate solution. It seems that the solver is not taking sufficiently small step sizes in order to correctly catch the transition in the ODE. I also tried adding in an event function with the expectation that the step size would correctly shrink close to this transition but it doesn’t seem to effect the accuracy, just outputs the time of the event.

The attached code solves the same ODE with manual rk4 and solve_ivp. The rk4 solutions converge correctly despite different initial conditions, and the solve_ivp solutions do not, despite very similar initial conditions.

Thanks

Reproducing Code Example

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

t = np.linspace(0,300,3001)

def dy(t, y):
    x = 0.02 * np.sin(t) + 0.4
    if x >= y:
        b = 0.2
    else:
        b = 0.02
    return b * (x - y)

def rk4(t,y,h):
    k1 =  h*dy(t, y)
    k2 =  h*dy(t, y + k1/2)
    k3 =  h*dy(t, y + k2/2)
    k4 =  h*dy(t, y + k3)
    return (k1 + 2*k2 + 2*k3 + k4)/6

def solve_rk4(dy,t_span,y0,max_step):
    dt = t_span[1]-t_span[0]
    sol = np.empty_like(t_span)
    sol[0] = y0
    sub_step_n = int(dt/max_step)
    for i,t in enumerate(t_span[:-1]):
        yi = sol[i]
        for j in range(sub_step_n):
            yi += rk4(t,yi,max_step)
        sol[i+1] = yi
    return sol

def event(t,y):
    return 0.02 * np.sin(t) + 0.4 - y

sol1 = solve_ivp(dy, t[[0,-1]], (0.42,), events=event, t_eval=t, rtol=1e-3, atol=1e-6).y.squeeze()
sol2 = solve_ivp(dy, t[[0,-1]], (0.4200001,), events=event, t_eval=t, rtol=1e-3, atol=1e-6).y.squeeze()
sol3 = solve_rk4(dy, t, 0.42, 0.1)
sol4 = solve_rk4(dy, t, 0.43, 0.1)
plt.plot(t,sol1,label='solve_ivp 0.42')
plt.plot(t,sol2,label='solve_ivp 0.4200001')
plt.plot(t,sol3,label='RK4 0.42')
plt.plot(t,sol4,label='RK4 0.43')
plt.legend()
plt.show()

Error message

NA

SciPy/NumPy/Python version information

blas_mkl_info: NOT AVAILABLE blis_info: NOT AVAILABLE openblas_info: libraries = [‘openblas’, ‘openblas’] library_dirs = [‘/usr/local/lib’] language = c define_macros = [(‘HAVE_CBLAS’, None)] runtime_library_dirs = [‘/usr/local/lib’] blas_opt_info: libraries = [‘openblas’, ‘openblas’] library_dirs = [‘/usr/local/lib’] language = c define_macros = [(‘HAVE_CBLAS’, None)] runtime_library_dirs = [‘/usr/local/lib’] lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: libraries = [‘openblas’, ‘openblas’] library_dirs = [‘/usr/local/lib’] language = c define_macros = [(‘HAVE_CBLAS’, None)] runtime_library_dirs = [‘/usr/local/lib’] lapack_opt_info: libraries = [‘openblas’, ‘openblas’] library_dirs = [‘/usr/local/lib’] language = c define_macros = [(‘HAVE_CBLAS’, None)] runtime_library_dirs = [‘/usr/local/lib’] Supported SIMD extensions in this NumPy install: baseline = SSE,SSE2,SSE3 found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2,AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL not found = AVX512_KNL