Differences when solving in SI and dimensionless units

[henry2004y]

In principle, the Lorentz equation in SI units

$$ \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = \frac{q}{m}\left( \mathbf{v}\times\mathbf{B} + \mathbf{E} \right) $$

can be normalized to

$$ \frac{\mathrm{d}\mathbf{v}^\prime}{\mathrm{d}t^\prime} = \mathbf{v}^\prime\times\mathbf{B}^\prime + \mathbf{E}^\prime $$

with the following transformation

$$ \begin{aligned} \mathbf{v} &= \mathbf{v}^\prime V_0 \ t &= t^\prime t_0 = t^\prime \Omega^{-1} = t^\prime \frac{m}{qB_0} \ \mathbf{B} &= \mathbf{B}^\prime B_0 \ \mathbf{E} &= \mathbf{E}^\prime E_0 = \mathbf{E}^\prime V_0 B_0 \end{aligned} $$

All the coefficients here are written in SI units. Often we would choose $V_0 = c$ (speed of light).

However, consider the following example

using TestParticle, OrdinaryDiffEq, StaticArrays
using TestParticle: c, qᵢ, mᵢ

const B₀ = 1e-8 # [T]
const U₀ = c    # [m/s]
const E₀ = U₀*B₀ # [V/m]
const Ω = abs(qᵢ) * B₀ / mᵢ
const t₀ = 1 / Ω  # [s]
const l₀ = U₀ * t₀ # [m]

B(x) = SA[0, 0, B₀]
E(x) = SA[1e-8, 0.0, 0.0]

B_normalize(x) = SA[0, 0, 1.0]
E_normalize(x) = SA[1e-8/E₀, 0.0, 0.0]
#E_normalize(x) = SA[0.0, 0.0, 0.0]

x0 = [0.0, 0.0, 0.0] # [m]
v0 = [0.0, 0.01c, 0.0] # [m/s]
stateinit1 = [x0..., v0...]
tspan1 = (0, 2π/Ω) # [s]
# Solve in SI units
param1 = prepare(E, B, species=Proton)
prob1 = ODEProblem(trace!, stateinit1, tspan1, param1)
sol1 = solve(prob1, Vern9(); reltol=1e-4, abstol=1e-6)

# Solve in dimensionless units
param2 = prepare(E_normalize, B_normalize, 1.0; species=User)
x0 ./= l₀
v0 ./= U₀
tspan2 = (0, 2π/Ω/t₀)
stateinit2 = [x0..., v0...]

prob2 = ODEProblem(trace_normalized!, stateinit2, tspan2, param2)
sol2 = solve(prob2, Vern9(); reltol=1e-4, abstol=1e-6)

println("positions: ")
println(sol1.u[end][1:3])
println(sol2.u[end][1:3] .* l₀)
println("velocities: ")
println(sol1.u[end][4:6])
println(sol2.u[end][4:6] .* U₀)

The two methods do not give the same results. Where did we do wrong?

[henry2004y]

Actually when we plot both trajectories in SI units, they are very close to each other:

using GLMakie

f = Figure(fontsize=18)
ax = Axis3(f[1, 1],
    xlabel = "x",
    ylabel = "y",
    zlabel = "z",
    aspect = :data,
)

lines!(ax, sol1)
trange = range(tspan2..., length=40)
xplot = sol2.(trange, idxs=1) .* l₀
yplot = sol2.(trange, idxs=2) .* l₀
zplot = sol2.(trange, idxs=3) .* l₀
lines!(ax, xplot, yplot, zplot, linestyle=:dashdot, linewidth=5)

f

So the differences are probably solely numerical. I've also tried with other E field values.

[TCLiuu]

If it is a case given in this issue, I think the difference between the results obtained using the Vern9 algorithm should be numerical error.

[henry2004y]

One side finding from this experiment is that tracing in dimensionless units seems to be faster, because the timesteps are much larger.