[FIXED] first_step parameter in solve_ivp doesn't work when large enough

Issue

I ran into a wird case where i tested the first_step parameter in the ode solver solve_ivp in scipy an when it crossed a into being 2 (ie first_step=2) it was just being ignored by the solver.

So the code that i tested is

from scipy.integrate import solve_ivp

def exponential_decay(t, y): return -0.5 * y

and when i'm calling the solver function with first_step=1.9, it works as expected:

>>> sol = solve_ivp(exponential_decay, [0, 10], [2], first_step=1.9)
>>> print(sol.t)
[ 0.          1.9         3.64913379  5.40785745  7.16661775  8.92690194
 10.        ]

Where it's correct and the second value in the list is 1.9, but when I run it with first_step=2 or for any number larger than 2 for that matter, it gives me the same result which seems to be contradictory to the documantation where it states "first_step: float or None, optional. Initial step size. Default is None which means that the algorithm should choose."

>>> sol = solve_ivp(exponential_decay, [0, 10], [2], first_step=2)
>>> print(sol.t)
[ 0.          1.743044    3.486088    5.24515667  7.00381091  8.76389907
 10.        ]

It's always 1.743044 for some reason and not 2.

Does anyone know why this happens and if there is a way to solve this problem?

Solution

You can change this behavior by adjusting the error tolerances atol and rtol. If you make the tolerances too loose, the whole integration process may break down.

Adaptive solvers take as step size the (guessed) optimal step size of the last step. Obviously at the start there is no last step, so the first step has to be even more heuristically guessed or passed as parameter.

What you observe is that at some point the large first step does not pass the error tolerances. The local error is allowed to be somewhat smaller than the tolerance, so small first step sizes pass through. But for larger errors the step is rejected and recomputed with the optimal guess. This guess tends to be largely independent of the too large first step parameter.

Change your test ODE to

def exponential_decay(t, y): print(t); return -0.5 * y

to see where the function is actually evaluated, note that the RK45 method has 6+1 FSAL stages.

Answered By - Lutz Lehmann