# Questions on the DynamicalODE integrators

**URL:** <https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721>\
**Category:** Modelling & Simulations\
**Tags:** question, diffeq\
**Created:** [March 28, 2023, 3:54pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721 "2023-03-28T15:54:03Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 28, 2023, 3:54pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/1 "2023-03-28T15:54:03Z")

</div>

Hi!

I have a question about the integrators in the differential equation community. The integrators for `DynamicalODEProblem` are listed in [Dynamical, Hamiltonian, and 2nd Order ODE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/latest/solvers/dynamical_solve/). While testing for the integrators, especially those symplectic ones, I found some puzzling results regarding energy conservation:

```julia
# Tracing charged particle in a static EM field.

using TestParticle
using OrdinaryDiffEq
using StaticArrays
using LinearAlgebra: ×
#using Plots

const B₀ = 1e-8 # [T]

function location!(dx, x, v, p::TestParticle.TPTuple, t)
   dx .= v
end

function lorentz!(dv, x, v, p::TestParticle.TPTuple, t)
   q, m, E, B = p
   dv .= q/m*(E(x, t) + v × (B(x, t)))
end

## Initialize field

function uniform_B(x)
   return SA[0, 0, B₀]
end

function uniform_E(x)
   return SA[0.0, 0.0, 0.0]
end

"Check energy conservation."
function getenergy(sol)
   vx, vy, vz = sol.u[1][4:6]
   KE1 = √(vx^2 + vy^2 + vz^2)
   vx, vy, vz = sol.u[end][4:6]
   KE2 = √(vx^2 + vy^2 + vz^2)
   KE1, KE2
end

## Initialize particles

x0 = [0.0, 0, 0]
v0 = [0.0, 1e2, 0.0]
stateinit = [x0..., v0...]
tspan_proton = (0.0, 2000.0)

param_proton = prepare(uniform_E, uniform_B, species=Proton)

## Solve for the trajectories

prob_p = DynamicalODEProblem(location!, lorentz!, x0, v0, tspan_proton, param_proton)

Ωᵢ = TestParticle.qᵢ * B₀ / TestParticle.mᵢ
Tᵢ = 2π / Ωᵢ
println("Number of gyrations: ", tspan_proton[2] / Tᵢ)

function main(prob_p)
   @info "Tsit5"
   @time sol_p = solve(prob_p, Tsit5(); save_idxs=[1,2,3,4,5,6]);

   KE1, KE2 = getenergy(sol_p)
   println("Energy: ", KE1, " to ", round(KE2, digits=4))
   @info "KahanLi6"
   @time sol_p = solve(prob_p, KahanLi6(); dt=0.01, save_idxs=[1,2,3,4,5,6]);

   KE1, KE2 = getenergy(sol_p)
   println("Energy: ", KE1, " to ", round(KE2, digits=4))
   @info "Trapezoid"
   @time sol_p = solve(prob_p, Trapezoid(); save_idxs=[1,2,3,4,5,6]);

   KE1, KE2 = getenergy(sol_p)
   println("Energy: ", KE1, " to ", round(KE2, digits=4))
   @info "Vern6"
   @time sol_p = solve(prob_p, Vern6(); save_idxs=[1,2,3,4,5,6]); # higher accuracy than Tsit5
   KE1, KE2 = getenergy(sol_p)
   println("Energy: ", KE1, " to ", round(KE2, digits=4))
end

main(prob_p)
#=
# Visualization
xyz = plot(sol_p, vars=(1,2,3), xlabel="x", ylabel="y", zlabel="z",
   aspect_ratio=:equal, label="proton", lw=2)
=#

```

```julia
julia> main(prob_p);
[ Info: Tsit5
  0.010167 seconds (61.31 k allocations: 3.842 MiB)
Energy: 100.0 to 109.8673
[ Info: KahanLi6
  0.458471 seconds (2.80 M allocations: 237.040 MiB, 50.05% gc time)
Energy: 100.0 to 887502.8753
[ Info: Trapezoid
  0.082830 seconds (1.05 M allocations: 25.859 MiB)
Energy: 100.0 to 100.5663
[ Info: Vern6
  0.004061 seconds (51.38 k allocations: 3.406 MiB)
Energy: 100.0 to 97.2829

```

This shows that the symplectic integrator `KahanLi6` does not conserve energy in my case, while others perform reasonably well (\< 10% energy change) for a single proton that gyrates for 300 periods in a static magnetic field. Did I set something wrong when using the symplectic integrators?

---

<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:** [March 28, 2023, 4:50pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/2 "2023-03-28T16:50:04Z")

</div>

> [@henry2004y](#):
>
> This shows that the symplectic integrator `KahanLi6` does not conserve energy in my case, while others perform reasonably well (\< 10% energy change) for a single proton that gyrates for 300 periods in a static magnetic field. Did I set something wrong when using the symplectic integrators?

Symplectic integrators do not conserve energy. Though that looks like an instability: does decreasing dt help?

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 28, 2023, 5:22pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/3 "2023-03-28T17:22:49Z")

</div>

Yes decreasing dt helps (`dt=0.01` in the original post):

```julia
[ Info: KahanLi6
# dt = 0.001
Energy: 100.0 to 248.2101
# dt = 0.0001
Energy: 100.0 to 109.5171

```

With `Tsit5` the time step quickly grows from 7.2233651437407934e-6 in the first step to about 0.5 in the 10th step.

---

<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:** [March 28, 2023, 5:42pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/4 "2023-03-28T17:42:12Z")

</div>

Okay, then that might’ve just been above the stability limit on KahanLi6.

---

<div class="post-metadata">

**Author:** ![henry2004y](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henry2004y/32/9284_2.png) [@henry2004y](https://discourse.julialang.org/u/henry2004y)\
**Post date:** [March 28, 2023, 5:44pm UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/5 "2023-03-28T17:44:35Z")

</div>

However it seems that the time step limit for all listed dynamical ODE integrators is pretty small compared with `Tsit5()` or other general methods? Is that expected?

---

<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:** [March 29, 2023, 2:44am UTC](https://discourse.julialang.org/t/questions-on-the-dynamicalode-integrators/96721/6 "2023-03-29T02:44:06Z")

</div>

Somewhat. It really depends on the kind of problem being solved but indeed some explicit methods have been designed with a larger or more robust stability region in mind.