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.
1 Like
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 (GitHub - JuliaMath/DoubleFloats.jl: math with more good bits) 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)) 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)
3 Likes
I have a few questions regarding your advice:
-
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.
-
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.
- 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]
β [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.
- 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
https://benchmarks.sciml.ai/html/DynamicalODE/Quadrupole_boson_Hamiltonian_energy_conservation_benchmark.html
1 Like