Questions on the DynamicalODE integrators

Specific Domains Modelling & Simulations
,
1

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. While testing for the integrators, especially those symplectic ones, I found some puzzling results regarding energy conservation:

# 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> 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?

2

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

3

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

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

4

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

5

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?

6

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.