# Non standard Hamiltonian dynamics - need high accuracy

I am solving dynamical ODEs that describe a Hamiltonian dynamics. However, I have that \dot{p}=f_1(q,p,t), \dot{q} = f_2(q,p,t), which does not apply to the requirements of the symplectic solvers available (or maybe am I missing something?)
(In general these equations I specified do not conserve energy, but my specific equation do conserve energy).

I need to solve these equations in very high accuracy (\sim 10^{-14}). I have tried to play with the algorithms such as Vern9(), Feagin12(), Feagin14() etc, but none of this give sufficient accuracy.
(I am not an expert, and I am not sure if my problem is stiff or not.)

How can I improve the accuracy of the solution?

One thing that might be necessary for this is using some form of higher precision arithmetic. eps(Float64)=2.2e-16, so that only gives you about 50 epsilon to play with, which isnβt a lot.

The thing is I am not reaching even 1e-7, even getting to 1e-12 would be much of an improvement.

Are you passing in arguments for atol and or rtol?

Yes. I am using 1e-14 in both.

Try using Rodas5 as your solver. If the issue is stiffness, that should fix it.

I would write it in dynamical ODE form and try DPRK12. Nystrom integrators utilize the partitioned structure for a bit more speed than standard methods while not assuming symplectic. Youβll likely want to change the types youβre using to not be Float64, since 1e-14 is just about at the edge of 64-bit accuracy range. Try changing your state and time types to BigFloat or Double64 (https://github.com/JuliaMath/DoubleFloats.jl) since that will likely help resolve the lower digits. Another thing to try is a newer solver in the ecosystem IRKGL16 (https://github.com/mikelehu/IRKGaussLegendre.jl) which is a 16th order Gauss-Legendre method: it could be useful here.

Not sure if stiffness has been identified as the issue here, but if it is, RadauIIA5() or ODEInterfaceDiffEq.jlβs radau() would be the methods to use for very high accuracy (the latter only supports 64-bit floats though, whereas RadauIIA5 can do arbitrary precision. But radau scales better to extreme accuracy cases)

2 Likes

I have a few questions regarding your advice:

1. How can I write in dynamical ODE form? I see that it does not apply to the first option ( Mathematical Specification of a Dynamical ODE Problem) since it is not of the form \dot{q} = p, not to the 2nd option ( Mathematical Specification of a 2nd Order ODE Problem) since the coupling is more involved, and not to the 3rd option ( Hamiltonian Problems) since the Hamiltonian is time dependent.

2. In order to use the Double64, I should just change the initial conditions and the time steps? i.e.

u0 = [Double64(q0), Double64(p0)]
tspan = (Double64(0.0), Double64(T))


Because I tried that, and it didnβt help the accuracy.

1. I had some error building IRKGL16
julia> Pkg.add(PackageSpec(url="https://github.com/mikelehu/IRKGaussLegendre.jl"))
Cloning git-repo https://github.com/mikelehu/IRKGaussLegendre.jl
Updating git-repo https://github.com/mikelehu/IRKGaussLegendre.jl
Updating git-repo https://github.com/mikelehu/IRKGaussLegendre.jl
Resolving package versions...
β Warning: julia version requirement for package ExamplePkg2 not satisfied
β @ Pkg.Operations /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/Pkg/src/Operations.jl:229
Installed LSODA ββββββββββββββββββ v0.6.1
Installed Formatting βββββββββββββ v0.4.1
Installed Espresso βββββββββββββββ v0.6.0
Installed ModelingToolkit ββββββββ v3.1.1
Installed Latexify βββββββββββββββ v0.13.1
Installed ODEInterfaceDiffEq βββββ v3.7.0
Installed ConstructionBase βββββββ v1.0.0
Installed JuliaVariables βββββββββ v0.2.0
Installed NameResolution βββββββββ v0.1.3
Installed ParameterizedFunctions β v5.3.0
Installed GeneralizedGenerated βββ v0.2.3
Installed PrettyPrint ββββββββββββ v0.1.0
Installed ODE ββββββββββββββββββββ v2.8.0
Installed SafeTestsets βββββββββββ v0.0.1
Installed Unitful ββββββββββββββββ v1.1.0
Installed MLStyle ββββββββββββββββ v0.3.1
Installed TaylorSeries βββββββββββ v0.10.3
Installed CanonicalTraits ββββββββ v0.2.1
Installed ODEInterface βββββββββββ v0.4.6
Installed TaylorIntegration ββββββ v0.8.2
Updating ~/.julia/environments/v1.3/Project.toml
[58bc7355] + IRKGaussLegendre v0.1.0 #master (https://github.com/mikelehu/IRKGaussLegendre.jl)
Updating ~/.julia/environments/v1.3/Manifest.toml
[a603d957] + CanonicalTraits v0.2.1
[187b0558] + ConstructionBase v1.0.0
[6912e4f1] + Espresso v0.6.0
[59287772] + Formatting v0.4.1
[6b9d7cbe] + GeneralizedGenerated v0.2.3
[58bc7355] + IRKGaussLegendre v0.1.0 #master (https://github.com/mikelehu/IRKGaussLegendre.jl)
[b14d175d] + JuliaVariables v0.2.0
[7f56f5a3] + LSODA v0.6.1
[23fbe1c1] + Latexify v0.13.1
[d8e11817] + MLStyle v0.3.1
[961ee093] + ModelingToolkit v3.1.1
[71a1bf82] + NameResolution v0.1.3
[c030b06c] + ODE v2.8.0
[54ca160b] + ODEInterface v0.4.6
[09606e27] + ODEInterfaceDiffEq v3.7.0
[65888b18] + ParameterizedFunctions v5.3.0
[8162dcfd] + PrettyPrint v0.1.0
[1bc83da4] + SafeTestsets v0.0.1
[92b13dbe] + TaylorIntegration v0.8.2
[6aa5eb33] + TaylorSeries v0.10.3
[1986cc42] + Unitful v1.1.0
Building LSODA ββββββββ ~/.julia/packages/LSODA/uaBH4/deps/build.log
Building ODEInterface β ~/.julia/packages/ODEInterface/xIa64/deps/build.log
β Error: Error building ODEInterface:
β ERROR: LoadError: Currently only gfortran is supported.
β Stacktrace:
β  [1] error(::String) at ./error.jl:33
β  [2] top-level scope at /Users/roi/.julia/packages/ODEInterface/xIa64/deps/build.jl:250
β  [3] include at ./boot.jl:328 [inlined]
β  [5] include(::Module, ::String) at ./Base.jl:31
β  [6] include(::String) at ./client.jl:424
β  [7] top-level scope at none:5
β in expression starting at /Users/roi/.julia/packages/ODEInterface/xIa64/deps/build.jl:249
β @ Pkg.Operations /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.3/Pkg/src/backwards_compatible_isolation.jl:649



So for now, I cannot really use this.

1. I tried RadauIIA5(), but it actually had worse accuracyβ¦

how low did you set the tolerances? 1e-24? There shouldnβt be a limit on accuracy if you do this with higher/arbitrary precision, so if that isnβt able to converge to what youβre saying is the solution, Iβd have to start questioning what youβre using as the truth values here.

itβs q' = f(p) so a little bit more expanded than that. Is that enough?

If not, it might be time to pull out https://docs.sciml.ai/latest/solvers/ode_solve/#Parallelized-Explicit-Extrapolation-Methods-1 , specifically ExtrapolationMidpointHairerWanner. You can set a pretty high max order and itβll multithread and hopefully do something good (your equation might be too small to multithread if itβs scalar, but I donβt know what youβre solving).

I set it as before to 1e-14.
I am comparing solutions to each other. I know that different initial conditions should maintain some qualities. For example, I know that the energy should be conserved and I see that the accuracy of that conservation is 1e-6.
I do not understand why the accuracy does not satisfy the limitation given by reltol and abstol.

I did not understand what you mean here.

Ahh yes, but the note that solver accuracy is step local while thatβs a global quantity. So you might have to do 1e-20 tolerances to get appropriate global energy error for example. All solvers are different, but preserving a quadratic quantity can be difficult if not built into the method. So try a lower tolerance.

One strategy that helps a ton is to combine with projections to the manifold. That can be done via: https://docs.sciml.ai/latest/features/callback_library/#Manifold-Conservation-and-Projection-1 . FWIW, when combined with a higher RK, thatβs the winning strategy in our energy conservation benchmarks:

https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html