I want to measure the performance of my own ODE integrator against SciPy RK45. Thus, I need to know exactly the number of right hand side function evaluations that RK45 uses.
Does anyone know if the number sol.nfev is accurate, i.e. without repititions? For example, if RK45 rejects a step size and repeats the step, is the initial evaluation f(t,x) of the step counted multiple times?
Also, the Dormand-Prince RK method used by RK45 has 7 stages, but in reality only uses 6 evaluations per step because of the "First Same As Last" property (the last stage is evaluated at the same point as the first stage of the next step). Is that accounted for in sol.nfev?
After some experiments I can confirm that RK45 counts accurately, i.e. repititions are not evaluated and not counted. For anyone wondering how the number sol.nfev can be explained:
Denote the number of RK steps reported by RK45 as num_steps = len(sol.t) - 1 and the number of rejected steps by num_rej. Then it holds
sol.nfev = (num_steps + num_rej) * 6 + 2
as every step needs 6 function evaluations, except the first one, which needs 7. Additionaly, there is one Euler step in the beginning to generate an initial step size (ergo the +2). Thanks @LutzLehmann for pointing that out.
Here is an example:
class Fun:
def __init__(self):
self.t = []
self.x = []
def __call__(self, t, x):
self.t.append(t)
self.x.append(x)
if t <= 1:
return - x ** 2
return 10 # discontinuity to provoke step rejection
t0 = 0
t1 = 1.01
x0 = [1]
sol = solve_ivp(Fun(), (t0, t1), x0)
print(sol.nfev)
print(len(f.t))
num_steps = len(sol.t) - 1
print(num_steps)
print(num_steps * 6 + 2)
From the output
56
56
7
44
we can see that there are 7 RK steps in the final solution sol but because of the discontinuity in the function, there also were 2 rejected steps (44 + 2 * 6 = 56).