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 ( since that will likely help resolve the lower digits. Another thing to try is a newer solver in the ecosystem IRKGL16 ( 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)


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=""))
   Cloning git-repo ``
  Updating git-repo ``
  Updating git-repo ``
 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 (
  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 (
  [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]
β”‚  [4] include_relative(::Module, ::String) at ./loading.jl:1105
β”‚  [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 , 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: . FWIW, when combined with a higher RK, that’s the winning strategy in our energy conservation benchmarks:

1 Like