# Use Symplectic Integrator for ODE for time reversability constraint

Hi everyone,

I would like to use Symplectic Integrators such as VelocityVerlet
or VerletLeapfrog from DifferentialEquations.jl to ensure the time-reversibility of the equation of motion.
Nonetheless, my problem can not be set as a DynamicalODEProblem.

My equation of motion is just:
\frac{du}{dt} = f2(u) where v = f2(u) is a model for the velocity which depends only of the position u )

Since my function f2 is not a function of the velocity v but only of the position u, I can not use these symplectic integrators.
Is there a way to get around this issue?

where did you get this idea? I think you flipped the equations. The docs say

This is a Partitioned ODE partitioned into two groups, so the functions should be specified as f1(dv,v,u,p,t) and f2(du,v,u,p,t) (in the inplace form), where f1 is independent of v (unless specified by the solver), and f2 is independent of u and t .

which matches your problem. Position depends on velocity, velocity depends on position and possibly time.

Sorry I should have been clearer, I mean that I don’t have 2 ODEs to solve (velocity v and position u ) but only one ODE to solve for the position u. Indeed, the velocity is given by a kinematic model and not by an ODE, i.e. I know in advance the velocity v at each position u : v=f2(u),

Therefore I don’t know what equation to set for \frac{dv}{dt}
I want to do something like this:

function f1(dv,v,u, p ,t)
dv[1:2] = [0.0, 0.0]
end

function f2(du,v,u,p,t)
#Velocity(u) known model for velocity
du[1:2] = Velocity(u)
end


Maybe I don’t need to use DynamicalODEProblem, I just want to ensure the time-reversibility of my solution.

Oh, then it’s not a good problem for these integrators. You’ll want to use a reversible method for first order ODEs, like Trapezoid(). You can set that to fixed point iteration if the problem isn’t stiff.

My model is an interpolated function using CubicSplineInterpolation of Interpolations.jl, therefore it doesn’t support ForwardDiff used by default in the implicit time scheme Trapezoid.
Do you know a derivative-free method I can use in Trapezoid()?
Yes, try Trapezoid(nlsolve = NLFunctional()) or Trapezoid(nlsolve = NLAnderson())