I am trying to understand the behavior of Diffeq.jl in a toy example, with the eventual goal of implementing networks of integrate-and-fire neurons. These questions are probably very basic, but any clarification would be much appreciated! The vector field of interest (parameters are pretty much random) is

function vf!(du, u,p,t)
du .= (2 .- u) ./ 3.0 .+ 0.2
end

I am comparing the behavior of two Euler stepping functions, one written by hand and one using the Diffeq.jl integrator interface.

function mystep!(du,u)
vf!(du,u,nothing,nothing)
u .= u .+ 0.01*du
end

When I run @benchmark mystep!(du0,u0) and @benchmark step!(integ), I get mean results 21.107 Ī¼s Ā± 42.136 Ī¼s and 27.882 Ī¼s Ā± 2.594 Ī¼s respectively. The increased variability in the first case seems to be attributable to one step in every 500 or so taking on the order of ~1 ms. This canāt be chalked up to any function compilation, as far as I understand.

Question 1: why does the by-hand step function occasionally take much longer, and how does the integrator in Diffeq.jl avoid this?

Question 2: what exactly is the additional overhead introduced by the integrator interface such that each step! takes longer on average than each mystep!?

Ah, thanks! That makes sense. However, this seems to result in the difference between step! and mystep! becoming even more dramatic (30 microsecs per step for step! vs 10 microsecs for mystep!). So my remaining question is: what is Diffeq.jl doing with those extra 20 microsecs?

Thatās simplified a bit, but you catch the drift. So first of all, this computation is so cheap that any(isnan,u) is the most expensive thing in the list . You can always turn that off:

but you can guess why a production ODE solver has defaults for exiting on NaN. But thereās 2 other extra O(n) kernel calls, why? This is explained in a very old video:

Basically whatās going on is that if you do this inverted form, update then calculate f, then you can have the derivative at both sides which defines a smooth interpolation over the solve. Note that the video explains how for most RK methods that are used there is a separate trick that actually makes it free. Itās really just Euler (and a few other low order less used methods), that require this trick, which then results in 2 extra copyto! operations.

So yes, if your f is only some trivial O(n) calculation, you can handwrite an Euler method to be 2n operations per step, while the default here is 4n operations + error checking.

So oh no, why is this generally not a problem? Letās take a look at the stepping code for the standard āsimpleā method: Tsit5

In this code you see that a single step of the 5th order RK method, with eigenvalue checks for stiffness detection, takes 13 O(n) operations and 6 f calls (6 instead of 7 due to FSAL). The FSAL gives one more O(n) operation by cutting out 1 f call. Then you have the O(n) operation for error checking / unstable check. That is how the other pieces usually get masked: the more complicated the method, and the complicated the f, the less the other details matter. Also (as mentioned in the video), here the FSAL handling gets rid of an f call so itās cheaper than not doing it (as opposed to an extra O(n) operation), and the interpolant construction is free by design. Thus if you think about it like this, the worst case scenario is to minimize the number of internal stages and minimize the complexity of f, which is precisely Eulerās method and your kind of example.

But really whatās going on here is that it doesnāt make sense for the simplest of simplest cases to use all of the machinery, which is why for this edge case we created SimpleDiffEq.jl

Now to finish off, other interesting things can come into play if n is really small, and thatās its own can of worms which is interesting tooā¦

Itās sometimes fun to see if you can get the fully generic version down to match the most basic implementation when u is a scalar and f is identity (@Oscar_Smith do you know if thereās some way to do a faster isnan?). But yeah, I hope that gives a bit of insight into how the simplest ODEs are somewhat of a corner case and why SimpleDiffEq.jl exists.