Time-dependent Hamiltonian for symplectic integrators in DifferentialEquations.jl

[EmadMasroor]

Hello,

I am using the DifferentialEquations.jl scaffolding to solve a Hamiltonian dynamical system. It is important to conserve the Hamiltonian to a reasonable accuracy, and so I am using symplectic integrators (SymplecticEuler or VerletLeapfrog work well for me) together with the constructor HamiltonianProblem{T}(H,p0,q0,tspan;kwargs...). This works well as long as there is no time-dependence.

However, I am at a loss as to how I can use this approach for a time-dependent Hamiltonian. The time-dependence is not very complicated (periodically switches), and I would like for my Hamiltonian function to look something like this:

function H(r,θ,p)
    switchtime = p
    if sin((π/switchtime)*t)>=0
        H = # some function of r,θ which does not depend on time
    else
        H = # another function of r,θ which does not depend on time
    end
end

But of course, since t is not a recognized variable, and H only takes the arguments p,q,params, this does not work.

I could easily accomplish this if I use the simple ODE constructor, e.g. ODEProblem(f::ODEFunction,u0,tspan,p=NullParameters();kwargs...) by including the explicit time dependence inside f(u,p,t) , but this would not allow me to use symplectic methods and my Hamiltonian diverges too quickly.

So my question is, is it possible to introduce simple time-dependence in the Hamiltonian for a HamiltonianProblem ? If the answer is no (which would be a bummer – I think this is the most elegant way to implement my ODEs), then can someone explain how I could do this with one of the other symplectic-friendly constructors, SecondOrderODEProblem or DynamicalODEProblem ?

My equations are ODEs of the form \dot{r} = f_1(r,\theta;t), \dot{\theta} = f_2(r,\theta;t), which I can re-cast as Hamiltonian with H = H(r,\theta;t). In the Hamiltonian formulation, r is my ‘momentum’ and \theta is my coordinate.

I’ll be happy to share more code and a MWE if it will help someone understand the question better!

[ChrisRackauckas]

@EmadMasroor

So my question is, is it possible to introduce simple time-dependence in the Hamiltonian for a HamiltonianProblem ? If the answer is no (which would be a bummer – I think this is the most elegant way to implement my ODEs), then can someone explain how I could do this with one of the other symplectic-friendly constructors, SecondOrderODEProblem or DynamicalODEProblem ?

I’m not sure the symplectic algorithms will work on a time-dependent Hamiltonian without order loss?

If you want to give it a try, you can define a DynamicalODEProblem where the two parts are defined by the gradients of the Hamiltonians. https://github.com/SciML/DiffEqPhysics.jl/blob/master/src/hamiltonian.jl should be very helpful.

[EmadMasroor]

Can I calculate the gradients of the Hamiltonian symbolically (I use SymEngine), and still get the benefits of a symplectic method? (another approach would be, for example, using ForwardDiff on the Hamiltonian like here, which seems a little daunting to me and I’d prefer explicitly defining the r.h.s. symbolically if I can)

Assuming I can do that, here’s what I tried. First, I defined the two functions f1 and f2 as in-place functions for the r.h.s.:

function f1!(dr,r,θ,p,t)
    dr = # some complicated function of r and θ.
end
function f2!(dθ,r,θ,p,t)
    dθ = # some complicated function of r and θ.
end

Then, I set up the problem like so:

r0 = 0.7
θ0 = 3π/4
p = 0
tspan = (0,1.0)
prob1=DynamicalODEProblem(f1!,f2!,r0,θ0,tspan,p)

And then I solve it using

sol1 = solve(prob1, SymplecticEuler(),dt=0.01,reltol=1e-6)

However, this gives an error:

MethodError: no method matching similar(::Float64, ::Type{Float64})

[ChrisRackauckas]

@EmadMasroor

function f1!(dr,r,θ,p,t)
    dr = # some complicated function of r and θ.
end
function f2!(dθ,r,θ,p,t)
    dθ = # some complicated function of r and θ.
end

if you write it like that it should be a mutating function. Did you mean:

function f1!(r,θ,p,t)
    # some complicated function of r and θ.
end
function f2!(r,θ,p,t)
    # some complicated function of r and θ.
end

?

[EmadMasroor]

Oh, right … that was a silly mistake. For one thing, the documentation made me think that it must be a mutating function. But it’s really unnecessary to have them be mutating functions here, so I’ll go ahead and change that.

Thanks! It’s working now, and I’ll go ahead and put in my time-dependence now. I appreciate you!

[EmadMasroor]

HI Chris,

It appears that using DynamicalODEProblem doesn’t preserve symplecticness when I introduce time-dependence. Here is a MWE:

The two functions on the R.H.S.:

function f1(r,θ,p,t)
    switchtime = p
    if sin((π/switchtime)*t)>=0
        dr = (-0.8*r*sin(θ)/sqrt(0.64 - 1.6*r*cos(θ) + r^2) -  0.8*(-r*sin(θ)/sqrt(0.64 - 1.6*r*cos(θ) + r^2) + 0.8*r*sin(θ)*(0.8 - r*cos(θ))/(0.64 - 1.6*r*cos(θ) + r^2)^(3/2)) - 1.0*(-r*sin(θ)/sqrt(1.5625 - 2.5*r*cos(θ) + r^2) + 1.25*r*sin(θ)*(1.25 - r*cos(θ))/(1.5625 - 2.5*r*cos(θ) + r^2)^(3/2)))/(r^2*sin(θ))
    else 
        dr = -(0.8*r*sin(θ)/sqrt(0.64 + 1.6*r*cos(θ) + r^2) - 0.8*(r*sin(θ)/sqrt(0.64 + 1.6*r*cos(θ) + r^2) - 0.8*r*sin(θ)*(0.8 + r*cos(θ))/(0.64 + 1.6*r*cos(θ) + r^2)^(3/2)) - 1.0*(r*sin(θ)/sqrt(1.5625 + 2.5*r*cos(θ) + r^2) - 1.25*r*sin(θ)*(1.25 + r*cos(θ))/(1.5625 + 2.5*r*cos(θ) + r^2)^(3/2)))/(r^2*sin(θ))
    end
end
function f2(r,θ,p,t)
    switchtime = p
    if sin((π/switchtime)*t)>=0
        dθ = -(1.0*(1 + (-1/2)*(2*r - 1.6*cos(θ))/sqrt(0.64 - 1.6*r*cos(θ) + r^2)) - 0.8*(cos(θ)/sqrt(0.64 - 1.6*r*cos(θ) + r^2) + (1/2)*(0.8 - r*cos(θ))*(2*r - 1.6*cos(θ))/(0.64 - 1.6*r*cos(θ) + r^2)^(3/2)) - 1.0*(cos(θ)/sqrt(1.5625 - 2.5*r*cos(θ) + r^2) + (1/2)*(1.25 - r*cos(θ))*(2*r - 2.5*cos(θ))/(1.5625 - 2.5*r*cos(θ) + r^2)^(3/2)))/(r*sin(θ))
    else
        dθ = (1.0*(1 + (-1/2)*(2*r + 1.6*cos(θ))/sqrt(0.64 + 1.6*r*cos(θ) + r^2)) - 0.8*(-cos(θ)/sqrt(0.64 + 1.6*r*cos(θ) + r^2) + (1/2)*(0.8 + r*cos(θ))*(2*r + 1.6*cos(θ))/(0.64 + 1.6*r*cos(θ) + r^2)^(3/2)) - 1.0*(-cos(θ)/sqrt(1.5625 + 2.5*r*cos(θ) + r^2) + (1/2)*(1.25 + r*cos(θ))*(2*r + 2.5*cos(θ))/(1.5625 + 2.5*r*cos(θ) + r^2)^(3/2)))/(r*sin(θ))
    end
end

which I solve for a short time interval using the commands:

prob1=DynamicalODEProblem(f1,f2,0.7,3π/4,(0,0.4),0.3)
sol1 = solve(prob1, SymplecticEuler(),dt=0.001,reltol=1e-6,abstol=1e-6)

The two functions physically represent a radial velocity and an azimuthal velocity, which I derived using SymEngine from an axisymmetric streamfunction. Both the functions ‘switch over’ between “right-handed” and “left-handed” versions, corresponding to a “right-handed” and “left-handed” Hamiltonian. While I do not use this Hamiltonian (or streamfunction) directly in the ODE, I do use it in order to check that my trajectories remain on streamlines like they are supposed to. That Hamiltonian is defined usingL

function H1(r,θ,p,t)
    switchtime = p
    if sin((π/switchtime)*t)>=0
        H = -0.8*(1 - (0.8 - r*cos(θ))/sqrt(0.64 - 1.6*r*cos(θ) + r^2)) -
    1.0*(1 - (1.25 - r*cos(θ))/sqrt(1.5625 - 2.5*r*cos(θ) + r^2)) + 1.0*(0.8 +
        r - sqrt(0.64 - 1.6*r*cos(θ) + r^2))
    else
        H = -(-0.8*(1 - (0.8 + r*cos(θ))/sqrt(0.64 + 1.6*r*cos(θ) + r^2)) - 1.0*(1 - (1.25 + r*cos(θ))/sqrt(1.5625 + 2.5*r*cos(θ) + r^2)) + 1.0*(0.8 + r - sqrt(0.64 + 1.6*r*cos(θ) + r^2)))
    end
    return H
end

I then wrote a function showTrajectory which plots the information in a convenient way, using the function H1 to calculate the value of the streamfunction:

function showTrajectory(soln)
    r = soln[1,:]
    θ = soln[2,:]
    p1 = plot(θ, r, proj = :polar, color=:black, seriestype=:scatter,
        markersize=0.05,legend=:false,size=(700,700))
    p2 = plot(soln.t,H1.(r,θ,0.3,soln.t),ylim=(-1,1), title="Streamfunction",
        seriestype=:scatter,markersize=0.1,legend=:none,size=(500,250))
    plot(p1,p2)
end

Yielding plots like this one (switchtime is 0.3, total integration time is 0.4)

You can see that the trajectory starts doing crazy things as soon as the ‘switching’ kicks in. The streamfunction should be preserved at a constant (but different) value after the switch has kicked in.

My definition of \psi has a distributed singularity and I eventually need this to behave relatively well in the neighborhood of a singularity, too…

By the way, this places me at basically the same position as when I was using the regular ODEProblem constructor.

Would it help to do manifold projection here? What else am I not doing right?

[ChrisRackauckas]

@EmadMasroor

Would it help to do manifold projection here? What else am I not doing right?

Probably. Are you sure the system is symplectic when accounting for events?

[EmadMasroor]

I’m not sure if I can say it is. I am following fluid particle trajectories in an axisymmetric domain according to a pre-defined streamfunction \psi_1, and every switchtime seconds I want the particles to flow according to \psi_2. I think this means that it should still be symplectic when accounting for this event…

I should probably mention here that once the integrator is working, I intend to add a callback which detects when a particle has been ‘sucked in’ to a sink, and then ejects it from a source, at the same streamline at which it entered the sink.

[ChrisRackauckas]

Symplectic is a fairly easy property to break. Even adaptive time stepping is something that’s not necessarily compatible with symplecticness (without adding time as a part of the Hamiltonian in special ways), so I would be surprised if these kinds of changes all keep that property intact. You might be better off with a very high order DPRKN method and projections applied every once in awhile.