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
N=2
#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)
matzdot=block_matrix(SR,[[dHdp.T],[-dHdq.T]])
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
#integration
def main():
z0 = [1.,-0.5,-2.,1.]
print("z0: ")
print(z0)
solverz = ode(dynz).set_integrator('dopri5')
solverz.set_initial_value(z0,0.)
solverz.integrate(Tf)
zf = solverz.y
print("zmid: ")
print(zf)
solverz.set_initial_value(zf,Tf)
solverz.integrate(0.)
print("z0: ")
print(solverz.y)
return solverz.y
main()
Here is the output:
z0:
[1.00000000000000, -0.500000000000000, -2.00000000000000, 1.00000000000000]
zmid:
[ 1.35249010e+00 -2.07489135e+00 -3.88286701e+12 8.51614447e+11]
z0:
[ 0.99956526 -0.50018124 -2.00172766 1.00161313]