# Use Symplectic Integrator for ODE for time reversability constraint

#1

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?

#2

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.

#3

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.

#4

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.

#5

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()?
Else I am thinking that I can probably provide the Jacobian of the model but I think we can only get Gradients with Interpolations.jl and not Jacobians.

#6

Yes, try Trapezoid(nlsolve = NLFunctional()) or Trapezoid(nlsolve = NLAnderson())

#7

Thank you @ChrisRackauckas for your help!

1 Like