1

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?

lueckem
  • 31
  • 4
  • 2
    You can create a list outside your right-hand side function and append all evaluated t to it. After integration you can check if that list has the same length as set(list). – Peter Meisrimel Dec 16 '20 at 15:14
  • 2
    In my sporadic tests of it, I never encountered a failed step. I'm not sure if a failed-step mechanism is even included in the code references of Shampine or Hairer. // There is a single additional Euler step performed in the first step, presumably to compute a sensible initial step size. See https://stackoverflow.com/questions/62343906/understanding-scipy-integrates-internal-behavior – Lutz Lehmann Dec 16 '20 at 15:52

1 Answers1

2

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).

lueckem
  • 31
  • 4