I write a test program to integrate foward on $[0,T_f]$ and then backward on $[T_f,0]$ from the endpoint of the forward integration an Hamiltonian system: $$ \dot q(t) = \frac{\partial H}{\partial p}(q(t),p(t)), \qquad \dot p(t) = -\frac{\partial H}{\partial q}(q(t),p(t)) $$ and the expression of $H$ (given in the program below doesn't really matter).

I compute symbolically the derivatives of the Hamiltonian and the jacobian of the dynamics using SageMath. I use the "dopri5" routine of the scipy python library for the numerical integration.

The program does what we want but there are some numerical errors since we don't get exactly the initial from which we started and I think theses errors come from the integrator ? We see that the intermediate point (zmid) is of order +1e12 so it causes the error, but what could be the other cause ?

To be more precise and accurate, I modified the following settings:

  • reduce atol, rtol and increase nsteps of the "dopri5" routine

Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale in the dynamics ...

from scipy.integrate import ode
from sage.all import *

#Dimension of the state variables
#Final time
Tf = 10.
#Symbolic variables of the Hamiltonian
var('q1 q2 p1 p2 v')
qs = [q1,q2]
ps = [p1,p2]
zs = qs + ps

#Hamiltonian function and jacobian of the dynamics
def hfun():
    H = p1*(cos(q1^2)+1/3*sin(q2)^2) + p2*(cos(q1)-2*sin(q2*q1)^2)
    dHdp = jacobian(H,ps)
    dHdq = jacobian(H,qs)
    jaczdot = jacobian(matzdot,zs)
    return H, jaczdot

#Hamiltonian dynamics
def dynz(t,z):
    H,aux = hfun()
    jacHp = jacobian(H,ps)
    jacHp = list(jacHp[0].subs({zs[i]:z[i] for i in range(0,2*N)}))
    dzdt = jacHp
    jacHq = jacobian(H,qs)
    dzdt.extend(list(-jacHq[0].subs({zs[i]:z[i] for i in range(0,2*N)})))
    return dzdt

def main():
    z0 = [1.,-0.5,-2.,1.]
    print("z0: ")

    solverz = ode(dynz).set_integrator('dopri5')

    zf = solverz.y
    print("zmid: ")

    print("z0: ")

    return solverz.y


Here is the output:

[1.00000000000000, -0.500000000000000, -2.00000000000000, 1.00000000000000]
[ 1.35249010e+00 -2.07489135e+00 -3.88286701e+12  8.51614447e+11]
[ 0.99956526 -0.50018124 -2.00172766  1.00161313]
  • 2
    $\begingroup$ Have you considered using a symplectic integrator? $\endgroup$ – Brian Borchers May 26 at 14:50

Do you have more suggestions to improve the precision and accuracy of the method ? What are the cause of such errors from the mathematical point of view ? For instance if we have different time scale in the dynamics ...

You need to use a reversible ODE solver method if you want to do this. I actually recently showed in a blog post that there are many cases where you can expect this to fail without a reversible integrator, such as when solving the Lorenz equation.

There are quite a few resources on time reversible methods:

and the best source is probably Hairer's Geometric Numerical Integration book which goes into detail in pretty much all ways, I highly recommend giving it a read. Basically, just because a numerical method is consistent (the error in a step goes to zero) doesn't mean that the method is stable (that the errors don't grow over time).

A symplectic method is reversible, but there are non-symplectic methods which are reversible as well (Trapezoid). Reversibility also requires that the set of time steps is chosen in a way that is symmetric, which essentially excludes any adaptive integration. Symplectic can have some extra advantages on Hamiltonian systems, so you're probably looking for a symplectic integrator.

SciPy doesn't include any reversible integrators, and does have the ability to turn off adaptive integration. Harier does have some Fortran codes for symplectic integrators but I'm not sure of any Python bindings. GSL has a few, but again not sure about Python bindings. DifferentialEquations.jl has a load of symplectic integrators, and a way to use it all from Python via diffeqpy. You can directly define a problem for use by symplectic integrators just by using HamiltonianProblem from DiffEqPhysics. A nice resource on using these can be found in the tutorials, and to use it from diffeqpy you essentially just append de. on package calls.

Footnote: this idea to do reverse integration seems to have recently gotten popular in deep learning and probabilistic programming circles for calculation of adjoints due to an award winning paper at NeurIPS, but you should think twice because it's not stable! Instead, one should use a system in which the calculation of the adjoints only requires forward passes, like is done in DifferentialEquations.jl or SUNDIALS where it's done by checkpointing or interpolation. So if you're calculating adjoints for Hamltionian Monte Carlo (HMC), keep this in mind.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.