Use Symplectic Integrator for ODE for time reversability constraint

Specific Domains Modelling & Simulations
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?

Thank you for your time,

2

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

http://docs.juliadiffeq.org/latest/types/dynamical_types.html#Mathematical-Specification-of-a-Dynamical-ODE-Problem-1

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

Thank you for your answer,

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

Thanks for the advice,

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

1 Like
7

Thank you @ChrisRackauckas for your help!

1 Like