Implement the drift-kick-drift form of the Leapfrog method

[efaulhaber]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

The kick-drift-kick form of the Leapfrog method (currently implemented as VerletLeapfrog) looks like this:

$$\begin{align} v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, f_1(v^0, u^0, t^0) \\\ u^1 &= u^0 + \Delta t\, f_2 \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t\, \right) \\\ v^1 &= v^{1/2} + \frac{1}{2} \Delta t\, f_1(v^{1}, u^{1}, t^0 + \Delta t) \end{align}$$

The evaluation of $f_2$ requires that it does not depend on $u$. If it does, we can add an intermediate step to compute $u^{1/2}$.
The second evaluation of $f_1$ requires that it does not depend on $v$, otherwise the whole thing doesn't work because we need $v^1$ to compute $v^1$.

It can easily be demonstrated that we lose second order accuracy like this:

using OrdinaryDiffEq, DiffEqDevTools, RecursiveArrayTools

function f1!(dv, v, u, p, t)
    dv .= v
    return dv
end

function f2!(du, v, u, p, t)
    du .= v
    return du
end

function analytic(y0, p, t)
    v0, u0 = y0.x
    ArrayPartition(v0 * exp(t), v0 * exp(t) - v0 + u0)
end

v0 = [2.0]
u0 = [0.5]

ff = DynamicalODEFunction(f1!, f2!; analytic = analytic)
prob = DynamicalODEProblem(ff, v0, u0, (0.0, 2.0))

dts = 1 .// 2 .^ (6:-1:3)
sim = test_convergence(dts, prob, VerletLeapfrog())
sim.𝒪est

Dict{Any, Any} with 3 entries:
  :l∞    => 0.931081
  :final => 0.935263
  :l2    => 0.963547

In this case, we can use the drift-kick-drift form, which is commonly used in Smoothed Particle Hydrodynamics, where viscosity makes $f_1$ depend on the velocity.
This form looks like this:

$$\begin{align} u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, f_2(v^0, u^0, t^0) \\\ v^1 &= v^0 + \Delta t\, f_1 \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t\, \right) \\\ u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, f_2(v^{1}, u^{1}, t^0 + \Delta t). \end{align}$$

Now, if $f_1$ depends on $v$, we can just add a half step for $v$ as well:

$$\begin{align} u^{1/2} &= u^0 + \frac{1}{2} \Delta t\, f_2(v^0, u^0, t^0) \\\ v^{1/2} &= v^0 + \frac{1}{2} \Delta t\, f_1(v^0, u^0, t^0) \\\ v^1 &= v^0 + \Delta t\, f_1 \left( v^{1/2}, u^{1/2}, t^0 + \frac{1}{2} \Delta t\, \right) \\\ u^1 &= u^{1/2} + \frac{1}{2} \Delta t\, f_2(v^{1}, u^{1}, t^0 + \Delta t). \end{align}$$

This is exactly what I implemented in this PR.
If $f_2$ does not depend on $u$, this can all be written down cleanly and yields the desired order of convergence:

julia> sim = test_convergence(dts, prob, LeapfrogDriftKickDrift());

julia> sim.𝒪est
Dict{Any, Any} with 3 entries:
  :l∞    => 1.95964
  :final => 1.95356
  :l2    => 1.98679

The cited paper by Verlet does not mention the leapfrog formulation, so I added a paper by Monaghan, which nicely explains the different forms of leapfrog in Section 5.3.

For a Smoothed Particle Hydrodynamics simulation with TrixiParticles.jl, this allows me to use a 5x larger time step than VerletLeapfrog, at (with #2559) 2x more f1 evaluations per step, resulting in a 2.7x speedup.