# Time-dependent Hamiltonian for symplectic integrators in DifferentialEquations.jl

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!

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.

2 Likes

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


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


?

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!

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?

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

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.

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.

3 Likes