# Time-dependent Hamiltonian for symplectic integrators in DifferentialEquations.jl
**URL:** https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987
**Category:** General Usage
**Tags:** package
**Created:** [May 22, 2020, 7:23pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987 "2020-05-22T19:23:44Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)
#### Post date: [May 22, 2020, 7:23pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/1 "2020-05-22T19:23:44Z")
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:
```julia
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!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 22, 2020, 7:38pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/2 "2020-05-22T19:38:41Z")
> [@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](https://github.com/SciML/DiffEqPhysics.jl/blob/master/src/hamiltonian.jl) should be very helpful.
---
### Author:  [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)
#### Post date: [May 22, 2020, 9:35pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/3 "2020-05-22T21:35:43Z")
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](https://nextjournal.com/sosiris-diffeq/diffeq-models#kepler-problem-and-sympletic-solvers), 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.:
```julia
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:
```julia
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
```julia
sol1 = solve(prob1, SymplecticEuler(),dt=0.01,reltol=1e-6)
```
However, this gives an error:
```julia
MethodError: no method matching similar(::Float64, ::Type{Float64})
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 22, 2020, 9:44pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/4 "2020-05-22T21:44:36Z")
> [@EmadMasroor](#):
>
> ```julia
> 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:
```julia
function f1!(r,θ,p,t)
# some complicated function of r and θ.
end
function f2!(r,θ,p,t)
# some complicated function of r and θ.
end
```
?
---
### Author:  [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)
#### Post date: [May 23, 2020, 5:37am UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/5 "2020-05-23T05:37:54Z")
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! 🙂
---
### Author:  [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)
#### Post date: [May 26, 2020, 6:51pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/6 "2020-05-26T18:51:37Z")
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.:**
```julia
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:
```julia
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
```julia
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:
```julia
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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2020, 9:22pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/7 "2020-05-26T21:22:41Z")
> [@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?
---
### Author:  [@EmadMasroor](https://discourse.julialang.org/u/EmadMasroor)
#### Post date: [May 26, 2020, 9:55pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/8 "2020-05-26T21:55:31Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2020, 9:58pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/9 "2020-05-26T21:58:14Z")
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.
---
### Author:  [@Jithin\_George](https://discourse.julialang.org/u/Jithin_George)
#### Post date: [February 13, 2026, 8:06pm UTC](https://discourse.julialang.org/t/time-dependent-hamiltonian-for-symplectic-integrators-in-differentialequations-jl/39987/10 "2026-02-13T20:06:55Z")
Actually for time-dependent Hamiltonians which are separable, the Candy Rozmus integrator is really good. See Candy, J., and W. Rozmus. “A symplectic integration algorithm for separable Hamiltonian functions.” _Journal of Computational Physics_ 92.1 (1991): 230-256. and McLachlan, Robert I., and Pau Atela. “The accuracy of symplectic integrators.” _Nonlinearity_ 5.2 (1992): 541.