# Non standard Hamiltonian dynamics - need high accuracy

**URL:** <https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811>\
**Category:** General Usage\
**Tags:** question, diffeq\
**Created:** [May 5, 2020, 2:43pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811 "2020-05-05T14:43:04Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [May 5, 2020, 2:43pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/1 "2020-05-05T14:43:04Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2020, 2:50pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/2 "2020-05-05T14:50:51Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [May 5, 2020, 2:55pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/3 "2020-05-05T14:55:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2020, 2:57pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/4 "2020-05-05T14:57:05Z")

</div>

Are you passing in arguments for `atol` and or `rtol`?

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [May 5, 2020, 2:59pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/5 "2020-05-05T14:59:45Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2020, 3:09pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/6 "2020-05-05T15:09:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 5, 2020, 3:35pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/7 "2020-05-05T15:35:21Z")

</div>

I would write it in [dynamical ODE form](https://docs.sciml.ai/latest/types/dynamical_types/) and try [DPRK12](https://docs.sciml.ai/latest/solvers/dynamical_solve/#Runge-Kutta-Nystr%C3%B6m-Integrators-1). 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` ([GitHub - JuliaMath/DoubleFloats.jl: math with more good bits](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` ([GitHub - mikelehu/IRKGaussLegendre.jl: Implicit Runge-Kutta Gauss-Legendre 16th order (Julia)](https://github.com/mikelehu/IRKGaussLegendre.jl)) which is a 16th order Gauss-Legendre method: it could be useful here.

> [@Oscar\_Smith](#):
>
> Try using `Rodas5` as your solver. If the issue is stiffness, that should fix it.

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)

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [May 5, 2020, 4:41pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/8 "2020-05-05T16:41:10Z")

</div>

I have a few questions regarding your advice:

1. How can I write in [dynamical ODE form](https://docs.sciml.ai/latest/types/dynamical_types/)? I see that it does not apply to the first option ( [Mathematical Specification of a Dynamical ODE Problem](https://docs.sciml.ai/latest/types/dynamical_types/#Mathematical-Specification-of-a-Dynamical-ODE-Problem-1)) since it is not of the form \dot{q} = p, not to the 2nd option ( [Mathematical Specification of a 2nd Order ODE Problem](https://docs.sciml.ai/latest/types/dynamical_types/#Mathematical-Specification-of-a-2nd-Order-ODE-Problem-1)) since the coupling is more involved, and not to the 3rd option ( [Hamiltonian Problems](https://docs.sciml.ai/latest/types/dynamical_types/#Hamiltonian-Problems-1)) 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.

```julia
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
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]
│ [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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 5, 2020, 5:01pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/9 "2020-05-05T17:01:13Z")

</div>

> [@roi.holtzman](#):
>
> Because I tried that, and it didn’t help the 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.

> [@roi.holtzman](#):
>
> since it is not of the form ˙q=p\dot{q} = p

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](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).

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [May 5, 2020, 5:51pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/10 "2020-05-05T17:51:26Z")

</div>

> [@ChrisRackauckas](#):
>
> how low did you set the tolerances? `1e-24` ?

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`.

> [@ChrisRackauckas](#):
>
> it’s `q' = f(p)` so a little bit more expanded than that

I did not understand what you mean here.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 5, 2020, 8:41pm UTC](https://discourse.julialang.org/t/non-standard-hamiltonian-dynamics-need-high-accuracy/38811/11 "2020-05-05T20:41:22Z")

</div>

> [@roi.holtzman](#):
>
> For example, I know that the energy should be conserved and I see that the accuracy of that conservation is `1e-6` .

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](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](https://benchmarks.sciml.ai/html/DynamicalODE/Henon-Heiles_energy_conservation_benchmark.html)

[https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole\_boson\_Hamiltonian\_energy\_conservation\_benchmark.html](https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.html)