Numerical Methods for Charged Particle Motion

This notebook demonstrates the implementation and analysis of numerical methods for simulating charged particle trajectories in magnetic fields using object-oriented programming.

%load_ext autoreload
%autoreload 2

# Add src to path for imports
import sys
sys.path.insert(0, '../src')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
from IPython.display import HTML

# Import numerical methods
from methods import BorisIntegrator, RK4Integrator, MultistepIntegrator, ABM4GuidingCenterIntegrator

# Import physical systems
from systems import System2D, TokamakSystem, GuidingCenter2D, GuidingCenterTokamak

# Import utility functions
from utils import (
    generate_reference_solution, compute_error, compute_convergence_study,
    compute_experimental_order, calculate_energy_errors, energy_error_guiding,
    plot_magnetic_fields, plot_trajectories, plot_errors, plot_energy_errors,
    plot_convergence_order, plot_work_precision,
    plot_error_guiding_center, plot_single_method_convergence, 
    plot_ABM_energy_error, plot_guiding_center_energy_errors,
)

Part A: Trajectory Visualization and Error Analysis

A.1 - Particle Trajectories

# Simulation parameters
h = 0.1
T = 1e5

# Initial conditions for 2D system
x0_2D = np.array([0.0, 1.0, 0.1])
v0_2D = np.array([0.09, 0.05, 0.20])

# Initial conditions for Tokamak system
x0_tokamak = np.array([1.20, 0.0, 0.0])
v0_tokamak = np.array([0, 4.816e-4, -2.059e-3])

# Create system instances
system_2D = System2D()
system_tokamak = TokamakSystem()

Magnetic Field Visualization

Before analyzing trajectories, we visualize the magnetic field structures for both systems.

import os
os.makedirs('../results/fields', exist_ok=True)

fig = plot_magnetic_fields(system_2D, system_tokamak)
fig.savefig('../results/fields/magnetic_fields.png', dpi=150, bbox_inches='tight')
plt.show()

# Boris method
boris = BorisIntegrator(h, T)
x_2D_boris, v_2D_boris = boris.integrate(x0_2D, v0_2D, system_2D)
x_tokamak_boris, v_tokamak_boris = boris.integrate(x0_tokamak, v0_tokamak, system_tokamak)

# Plot trajectories
trajectories_boris = {
    '2D': x_2D_boris,
    'tokamak': x_tokamak_boris,
    'method_name': 'Boris'
}
fig = plot_trajectories(trajectories_boris)
fig.savefig('../results/trajectories/boris_trajectories.png', dpi=300, bbox_inches='tight')
plt.show()

Boris Integration: 100%|██████████| 1000000/1000000 [01:05<00:00, 15292.43it/s]
Boris Integration: 100%|██████████| 1000000/1000000 [01:05<00:00, 15292.43it/s]
Boris Integration: 100%|██████████| 1000000/1000000 [01:08<00:00, 14681.56it/s]

# RK4 method
rk4 = RK4Integrator(h, T)
x_2D_rk4, v_2D_rk4 = rk4.integrate(x0_2D, v0_2D, system_2D)
x_tokamak_rk4, v_tokamak_rk4 = rk4.integrate(x0_tokamak, v0_tokamak, system_tokamak)

# Plot trajectories
trajectories_rk4 = {
    '2D': x_2D_rk4,
    'tokamak': x_tokamak_rk4,
    'method_name': 'RK4'
}
fig = plot_trajectories(trajectories_rk4)
fig.savefig('../results/trajectories/rk4_trajectories.png', dpi=300, bbox_inches='tight')
plt.show()

RK4 Integration: 100%|██████████| 1000000/1000000 [02:34<00:00, 6489.92it/s]
RK4 Integration: 100%|██████████| 1000000/1000000 [02:34<00:00, 6489.92it/s]
RK4 Integration: 100%|██████████| 1000000/1000000 [04:02<00:00, 4129.84it/s]

# Multistep method
multistep = MultistepIntegrator(h, T)
x_2D_multi, v_2D_multi = multistep.integrate(x0_2D, v0_2D, system_2D)
x_tokamak_multi, v_tokamak_multi = multistep.integrate(x0_tokamak, v0_tokamak, system_tokamak)

# Plot trajectories
trajectories_multi = {
    '2D': x_2D_multi,
    'tokamak': x_tokamak_multi,
    'method_name': 'Multistep'
}
fig = plot_trajectories(trajectories_multi)
fig.savefig('../results/trajectories/multistep_trajectories.png', dpi=300, bbox_inches='tight')
plt.show()

Multistep Integration: 100%|██████████| 999993/999993 [00:54<00:00, 18241.65it/s]
Multistep Integration: 100%|██████████| 999993/999993 [00:54<00:00, 18241.65it/s]
Multistep Integration: 100%|██████████| 999993/999993 [00:58<00:00, 16962.53it/s]

The above plots show the trajectories of a charged particle in both a 2D magnetic field and a tokamak magnetic field using three different numerical methods: the Boris method, the RK4 method, and a 9-step method. The 2D system can be thought of as a horizontal cross-section of the tokamak system, the toroidal plane. The Tokamak system can be thought of as a 3D extension of the 2D system, where we view the magnetic field lines in the poloidal plane, the plane perpendicular to the toroidal direction of the tokamak. For both systems, the two first coordinates of the solutions are plotted against each other to visualize the particle's motion in the plane. For the tokamak system, the two first coordinates are flattened to visualize the three dimensional motion in a two dimensional plot. The particle follows helical trajectories in both cases, as expected for charged particles in magnetic fields. This helical path can be further simplified to into two components:

• A circular motion around the magnetic field lines in the poloidal plane.
• A circular motion along the magnetic field lines in the toroidal plane.

When viewing the tokamak system from the side (the poloidal plane), the trajectory appears as a banana shape, which is characteristic of particle motion in tokamaks due to the combined effects of the magnetic field configuration and the particle's velocity components. Tweaking the initial velocity components will change the shape and size of the banana orbits, and could even make the particle complete full loops around the tokamak.

All 3 methods, give similar results for both the 2D and tokamak systems. We see little to no difference in the trajectories produced by the different methods, apart from some drift in the 2D case. The error plots below give a better understanding of the differences between the methods.

A.2 - Position Error Analysis

# Generate reference solutions
print("Generating reference solutions...")
sol_2D = generate_reference_solution(x0_2D, v0_2D, T, system_2D)
sol_tokamak = generate_reference_solution(x0_tokamak, v0_tokamak, T, system_tokamak)

# Compute errors for all methods
print("Computing errors...")
err_2D_boris = compute_error(T, x_2D_boris, sol_2D, h)
err_tokamak_boris = compute_error(T, x_tokamak_boris, sol_tokamak, h)

err_2D_rk4 = compute_error(T, x_2D_rk4, sol_2D, h)
err_tokamak_rk4 = compute_error(T, x_tokamak_rk4, sol_tokamak, h)

err_2D_multi = compute_error(T, x_2D_multi, sol_2D, h)
err_tokamak_multi = compute_error(T, x_tokamak_multi, sol_tokamak, h)

# Create error dictionary for plotting
errors_dict = {
    '2D': {
        'Boris': err_2D_boris,
        'RK4': err_2D_rk4,
        'Multistep': err_2D_multi
    },
    'tokamak': {
        'Boris': err_tokamak_boris,
        'RK4': err_tokamak_rk4,
        'Multistep': err_tokamak_multi
    }
}

# Plot errors
fig = plot_errors(errors_dict, h, T)
fig.savefig('../results/errors/position_errors.png', dpi=300, bbox_inches='tight')
plt.show()

Generating reference solutions...
Computing errors...
Computing errors...

Position errors for the 2D and tokamak systems, for all 3 methods, are shown in the plots above. The position error is calculated as the Euclidean distance between the numerical solution and a reference solution obtained using a very small time step with the RK45 method.

From the plots, we can observe the following:

• The two symplectic methods, Boris and the 9-step method, are both bounded in their position errors over time for both the 2D and tokamak systems. This indicates that these methods are stable and conserve certain properties of the system, such as phase space volume, over long simulations.
• The RK4 method, on the other hand, shows a position error that grows over time for both systems. This growth in error suggests that the RK4 method may not be as stable as the symplectic methods for long-term simulations, although better in the short term.

Overall, we see periodically lower errors for all methods. This is likely due to the periodic nature of the particle trajectories, where the phase of the motion in the methods occasionally aligns closely with that of the reference solution, leading to reduced errors at those times.

A.3 - Energy Conservation

# Calculate energy errors
energy_errors_2D = calculate_energy_errors(
    {'Boris': (x_2D_boris, v_2D_boris),
     'RK4': (x_2D_rk4, v_2D_rk4),
     'Multistep': (x_2D_multi, v_2D_multi)},
    T, h, system_2D
)

energy_errors_tokamak = calculate_energy_errors(
    {'Boris': (x_tokamak_boris, v_tokamak_boris),
     'RK4': (x_tokamak_rk4, v_tokamak_rk4),
     'Multistep': (x_tokamak_multi, v_tokamak_multi)},
    T, h, system_tokamak
)

energy_errors_dict = {
    '2D': energy_errors_2D,
    'tokamak': energy_errors_tokamak
}

# Plot energy errors
fig = plot_energy_errors(energy_errors_dict, h, T)
fig.savefig('../results/errors/energy_errors.png', dpi=300, bbox_inches='tight')
plt.show()

The energy error for the symplectic methods remains bounded over time for both the 2D and tokamak systems, although with some oscillations. This indicates that these methods are effective at conserving energy, a key property of Hamiltonian systems, over long simulations. The oscillations in energy error are typical for symplectic integrators and reflect the trade-off between accuracy and conservation properties. The RK4 method, however, shows a clear drift in energy error over time for both systems. This drift indicates that the RK4 method does not conserve energy as effectively as the symplectic methods, leading to an accumulation of error in the energy calculation over long simulations.

Our multistep method behaves in accordance with theorem 3.2 from Hairer and Lubich (2017), which states that along numerical solutions of s-stable symmetric methods of order $p$, the total energy is conserved up to an order of $O(h^p)$ over times $O(h^{-p-2})$. This is evident in our results, where the energy error remains bounded, within a thick band of approximate thickness $O(1e-5)$, just with different constants for the 2D system and the Tokamak system. It does not exhiibit drift either, similar to Figure 5.1 from Hairer and Lubich (2017). This is further supported by our experiments with different step sizes, where smaller step sizes lead to smaller energy errors, consistent with the $O(h^p)$ behavior predicted by the theorem.

Part B: Convergence Analysis

B.1 - Convergence Order

# Convergence study parameters
T_conv = 1e3  # Shorter time for convergence study
h_base = 0.1
num_refinements = 3

# Generate reference solutions for convergence study
sol_2D_conv = generate_reference_solution(x0_2D, v0_2D, T_conv, system_2D)
sol_tokamak_conv = generate_reference_solution(x0_tokamak, v0_tokamak, T_conv, system_tokamak)

# Convergence studies for 2D system
print("Running convergence studies for 2D system...")

hs_2D_boris, errors_2D_boris, N_2D_boris, times_2D_boris = compute_convergence_study(BorisIntegrator, x0_2D, v0_2D, system_2D, sol_2D_conv, T_conv, h_base, num_refinements)
eocs_2D_boris = compute_experimental_order(errors_2D_boris, hs_2D_boris)

hs_2D_rk4, errors_2D_rk4, N_2D_rk4, times_2D_rk4 = compute_convergence_study(RK4Integrator, x0_2D, v0_2D, system_2D, sol_2D_conv, T_conv, h_base, num_refinements)
eocs_2D_rk4 = compute_experimental_order(errors_2D_rk4, hs_2D_rk4)

hs_2D_multi, errors_2D_multi, N_2D_multi, times_2D_multi = compute_convergence_study(MultistepIntegrator, x0_2D, v0_2D, system_2D, sol_2D_conv, T_conv, h_base, num_refinements)
eocs_2D_multi = compute_experimental_order(errors_2D_multi, hs_2D_multi)

# Convergence studies for Tokamak system
print("Running convergence studies for Tokamak system...")

hs_tokamak_boris, errors_tokamak_boris, N_tokamak_boris, times_tokamak_boris = compute_convergence_study(BorisIntegrator, x0_tokamak, v0_tokamak, system_tokamak, sol_tokamak_conv, T_conv, h_base, num_refinements)
eocs_tokamak_boris = compute_experimental_order(errors_tokamak_boris, hs_tokamak_boris)

hs_tokamak_rk4, errors_tokamak_rk4, N_tokamak_rk4, times_tokamak_rk4 = compute_convergence_study(RK4Integrator, x0_tokamak, v0_tokamak, system_tokamak, sol_tokamak_conv, T_conv, h_base, num_refinements)
eocs_tokamak_rk4 = compute_experimental_order(errors_tokamak_rk4, hs_tokamak_rk4)

hs_tokamak_multi, errors_tokamak_multi, N_tokamak_multi, times_tokamak_multi = compute_convergence_study(MultistepIntegrator, x0_tokamak, v0_tokamak, system_tokamak, sol_tokamak_conv, T_conv, h_base, num_refinements)
eocs_tokamak_multi = compute_experimental_order(errors_tokamak_multi, hs_tokamak_multi)

Running convergence studies for 2D system...

Convergence Study: 3it [00:06,  2.10s/it]
Convergence Study: 3it [00:06,  2.10s/it]
Convergence Study: 3it [00:14,  4.96s/it]
Convergence Study: 0it [00:00, ?it/s]
Convergence Study: 3it [00:04,  1.36s/it]
Convergence Study: 3it [00:04,  1.36s/it]

Running convergence studies for Tokamak system...

Convergence Study: 3it [00:06,  2.06s/it]
Convergence Study: 3it [00:06,  2.06s/it]
Convergence Study: 3it [00:14,  4.97s/it]
Convergence Study: 3it [00:14,  4.97s/it]
Convergence Study: 3it [00:04,  1.44s/it]
Convergence Study: 3it [00:04,  1.44s/it]

# Display results in tables
results_df_2D = pd.DataFrame({
    'h': hs_2D_boris,
    'Boris Error': errors_2D_boris,
    'Boris EOC': eocs_2D_boris,
    'RK4 Error': errors_2D_rk4,
    'RK4 EOC': eocs_2D_rk4,
    'Multistep Error': errors_2D_multi,
    'Multistep EOC': eocs_2D_multi,
    'Boris Time (s)': times_2D_boris,
    'RK4 Time (s)': times_2D_rk4,
    'Multistep Time (s)': times_2D_multi,
})

print("\nConvergence Study Results - 2D System:")
display(results_df_2D)

results_df_tokamak = pd.DataFrame({
    'h': hs_tokamak_boris,
    'Boris Error': errors_tokamak_boris,
    'Boris EOC': eocs_tokamak_boris,
    'RK4 Error': errors_tokamak_rk4,
    'RK4 EOC': eocs_tokamak_rk4,
    'Multistep Error': errors_tokamak_multi,
    'Multistep EOC': eocs_tokamak_multi,
    'Boris Time (s)': times_tokamak_boris,
    'RK4 Time (s)': times_tokamak_rk4,
    'Multistep Time (s)': times_tokamak_multi,
})

print("\nConvergence Study Results - Tokamak System:")
display(results_df_tokamak)

Convergence Study Results - 2D System:

	h	Boris Error	Boris EOC	RK4 Error	RK4 EOC	Multistep Error	Multistep EOC	Boris Time (s)	RK4 Time (s)	Multistep Time (s)
0	0.100	0.072403	NaN	2.836677e-04	NaN	0.011710	NaN	0.899319	2.083793	0.534446
1	0.050	0.018358	1.979620	1.080579e-05	4.714326	0.000699	4.066820	1.723635	4.146653	1.088193
2	0.025	0.004594	1.998652	4.502988e-07	4.584777	0.000043	4.030096	3.407846	8.296397	2.178048

Convergence Study Results - Tokamak System:

	h	Boris Error	Boris EOC	RK4 Error	RK4 EOC	Multistep Error	Multistep EOC	Boris Time (s)	RK4 Time (s)	Multistep Time (s)
0	0.100	0.001214	NaN	8.173172e-07	NaN	1.609712e-04	NaN	0.859745	2.215161	0.606354
1	0.050	0.000306	1.986393	5.195859e-08	3.975462	9.957382e-06	4.014893	1.734650	4.230139	1.175879
2	0.025	0.000077	1.999139	3.144178e-09	4.046608	6.205354e-07	4.004181	3.418559	8.301303	2.360825

# Prepare convergence data
results_dict_2D = {
    'Boris': {'hs': hs_2D_boris, 'errors': errors_2D_boris},
    'RK4': {'hs': hs_2D_rk4, 'errors': errors_2D_rk4},
    'Multistep': {'hs': hs_2D_multi, 'errors': errors_2D_multi}
}

results_dict_tokamak = {
    'Boris': {'hs': hs_tokamak_boris, 'errors': errors_tokamak_boris},
    'RK4': {'hs': hs_tokamak_rk4, 'errors': errors_tokamak_rk4},
    'Multistep': {'hs': hs_tokamak_multi, 'errors': errors_tokamak_multi}
}

# Plot convergence order
fig = plot_convergence_order(results_dict_2D, results_dict_tokamak)
fig.savefig('../results/convergence/convergence_order.png', dpi=300, bbox_inches='tight')
plt.show()

Plotting the global errors against the time step size h on a log-log scale allows us to visualize the convergence rates of the different numerical methods used for simulating the charged particle motion in both the 2D and tokamak systems. From the plots, we can observe the following:

• The Boris method shows a convergence rate of approximately 2, indicating that it is a second-order method. This means that as the time step size h is halved, the global error decreases by a factor of about 4. The RK4- and 9-step methods both exhibit a convergence rate of approximately 4, indicating that they are fourth-order methods, as expected.
• The difference in global errors between the methods is likely a result of their inherent accuracy and stability properties. Higher-order methods like RK4 and the 9-step method generally provide more accurate results for a given time step size compared to lower-order methods like the Boris method. This is because higher-order methods can better capture the dynamics of the system with fewer discretization errors.

B.2 - Work-Precision Diagram

# Prepare work-precision data
work_precision_data_2D = {
    'Boris': {'errors': errors_2D_boris, 'times': times_2D_boris},
    'RK4': {'errors': errors_2D_rk4, 'times': times_2D_rk4},
    'Multistep': {'errors': errors_2D_multi, 'times': times_2D_multi}
}

work_precision_data_tokamak = {
    'Boris': {'errors': errors_tokamak_boris, 'times': times_tokamak_boris},
    'RK4': {'errors': errors_tokamak_rk4, 'times': times_tokamak_rk4},
    'Multistep': {'errors': errors_tokamak_multi, 'times': times_tokamak_multi}
}

# Plot work-precision diagram
fig = plot_work_precision(work_precision_data_2D, work_precision_data_tokamak)
fig.savefig('../results/convergence/work_precision.png', dpi=300, bbox_inches='tight')
plt.show()

The work-precision diagrams illustrate the relationship between computational work (measured in function evaluations) and the achieved global error for the different numerical methods applied to both the 2D and tokamak systems. From the diagrams, we can observe the following:

• The Boris method, being a second-order method, generally requires fewer function evaluations to achieve a given level of accuracy compared to the higher-order methods. This makes it computationally efficient for moderate accuracy requirements.
• The RK4 method, while more accurate for a given time step size, requires significantly more function evaluations to reach the same error levels as the Boris method. This indicates that while RK4 is more accurate, it is also more computationally intensive.
• The 9-step method strikes a balance between accuracy and computational work, being the most efficient among the three methods, whilst still achieving good accuracy.

These observations rely highly on the implementation of the methods, as for example calculating the jacobian matrix for the 9-step method analytically yields a significant speedup compared to using autograd. With our current implementation, the 9-step method, with a good implementation, should be the best choice for simulating charged particle motion in magnetic fields when considering both accuracy and computation time. However, optimizing the implementations could change the results of these diagrams.

Part C: Guiding Center Approximation (Optional)

Guiding Center Method

We got inspired by the work of Hoppe, Iantchenko, and Stranberg (2015) on guiding center dynamics in tokamaks and decided to implement a simple guiding center model for the tokamak system. The guiding center approximation simplifies the motion of charged particles in a magnetic field by averaging out the fast gyromotion around magnetic field lines, allowing us to focus on the slower drift motion of the particle's guiding center. In other words, instead of tracking the particle's full 6D-phase-space position ($x$, $v$), we now only need to track its 4D guiding center position ($x$, $v_\parallel$). The reduced model is:

$$\frac{d\mathbf{X}}{dt} = \frac{v_\parallel \mathbf{B}}{|\mathbf{B}|} + \frac{\mu}{q|\mathbf{B}|^2} \mathbf{B} \times \nabla |\mathbf{B}|$$

$$ \frac{dv_\parallel}{dt} = -\frac{\mu}{m} \mathbf{B} \cdot \nabla |\mathbf{B}| $$

$$ \mu = \frac{m v_\perp^2}{2|\mathbf{B}|} $$

where $\mathbf{X}$ is the guiding center position, $v_\parallel$ is the velocity component parallel to the magnetic field, $v_\perp$ is the initial velocity component perpendicular to the magnetic field, $\mu$ is the magnetic moment, $q$ is the particle charge, $m$ is the particle mass, and $\mathbf{B}$ is the magnetic field vector.

To solve this system, we cannot use the Boris method, as that relies on splitting the Lorentz force, which is not present in the guiding center equations. We also wanted to use a new method, with a slightly different approach. Therefore we implemented the Adams-Bashforth-Moulton (ABM) predictor-corrector method, as talked about briefly in the lectures. This method is a 4th order multi-step approach that uses previous points explicitly to predict the next point and then corrects it implicitly, providing a good balance between accuracy and computational efficiency.

# Guiding center simulation parameters
h_gc = 10.0
T_gc = 1e5

# Create guiding center systems
gc_system_2D = GuidingCenter2D()
gc_system_tokamak = GuidingCenterTokamak()

# Create ABM4 integrator
abm4 = ABM4GuidingCenterIntegrator(h_gc, T_gc)

# Run simulations
print("Running ABM4 for 2D system...")
Xs_2D, vpars_2D, mu_2D = abm4.integrate(x0_2D, v0_2D, gc_system_2D)

print("Running ABM4 for Tokamak system...")
Xs_tokamak, vpars_tokamak, mu_tokamak = abm4.integrate(x0_tokamak, v0_tokamak, gc_system_tokamak)

# Plot trajectories
trajectories_abm4 = {
    '2D': Xs_2D,
    'tokamak': Xs_tokamak,
    'method_name': 'ABM4'
}
fig = plot_trajectories(trajectories_abm4)
fig.savefig('../results/guiding_center/gc_trajectories.png', dpi=300, bbox_inches='tight')
plt.show()

Running ABM4 for 2D system...

Guiding, h=10.0: 100%|██████████| 9993/9993 [00:04<00:00, 2307.87it/s]
Guiding, h=10.0: 100%|██████████| 9993/9993 [00:04<00:00, 2307.87it/s]

Running ABM4 for Tokamak system...

Guiding, h=10.0: 100%|██████████| 9993/9993 [00:12<00:00, 832.31it/s]

The plots above show the guiding center trajectories in the tokamak magnetic field, i.e. the center of the banana orbits. The guiding center motion captures the overall drift of the particle in the tokamak, without resolving the fast gyromotion around the magnetic field lines. This results in a smoother trajectory that still reflects the influence of the magnetic field in the toroidal plane.

The benefit of using the guiding center appoximation is that it allows for larger time steps in the simulation, as we are not constrained by the need to resolve the fast gyromotion. This can lead to significant computational savings, especially for long-term simulations of particle dynamics in tokamaks. In our simulations we used 100x larger time steps for the 2D system and 200x larger time steps for the tokamak system compared to the full particle motion simulations, while still capturing the essential dynamics of the guiding center motion.

Guiding Center Error Analysis

We can compare the guiding center approximation results with full particle orbit simulations to assess accuracy.

sol_2D_guiding_center = generate_reference_solution(x0_2D, v0_2D, T_gc, gc_system_2D)
sol_tokamak_guiding_center = generate_reference_solution(x0_tokamak, v0_tokamak, T_gc, gc_system_tokamak)

err_gc_2D = compute_error(T_gc, Xs_2D, sol_2D_guiding_center, h_gc)
err_gc_tokamak = compute_error(T_gc, Xs_tokamak, sol_tokamak_guiding_center, h_gc)

# Plot errors
fig = plot_error_guiding_center(err_gc_2D, err_gc_tokamak, h_gc, h_gc, T_gc)
fig.savefig('../results/guiding_center/gc_errors.png', dpi=300, bbox_inches='tight')
plt.show()

The position error plot shows the global error of the guiding center ABM method as a function of the time step size h. The error increases as time increaes for both the 2D and the tokamak systems, similar to the RK4 method for the full particle motion, although at a slower rate for the tokamak system here. This hints that the Adams-Bashforth-Moulton method is not symplectic (which it is not). The error for the 2D system is significantly higher than for the tokmak system, which likely is a twofold effect of the high step size used and the guiding center moving way faster in the 2D system compared to the tokamak system.

Guiding Center Energy Conservation

energy_err_gc_2D = energy_error_guiding(Xs_2D, vpars_2D, gc_system_2D.B, gc_system_2D.U, mu_2D)
energy_err_gc_tokamak = energy_error_guiding(Xs_tokamak, vpars_tokamak, gc_system_tokamak.B, gc_system_tokamak.U, mu_tokamak)

fig = plot_ABM_energy_error(energy_err_gc_2D, energy_err_gc_tokamak, h_gc, h_gc, T_gc)
fig.savefig('../results/guiding_center/gc_energy_errors.png', dpi=300, bbox_inches='tight')
plt.show()

The energy error plot shows the energy conservation error for the guiding center ABM method in both the 2D and tokamak systems. As expected for a non-symplectic method, the energy increases over time for the 2D system. However, for the tokamak system, the energy error remains bounded over time, similar to the symplectic methods used for the full particle motion. This, in addition to the very small errors observed, seems to be the result of something else. The energy in the guiding center model is given by

$$ E = \tfrac12 v_\parallel^2 + \mu |\mathbf{B}| + U(\mathbf{X}). $$

Differentiating along the trajectory gives

$$ \frac{dE}{dt} = v_\parallel \dot v_\parallel + \mu\, \nabla |\mathbf{B}| \cdot \dot{\mathbf{X}} + \nabla U \cdot \dot{\mathbf{X}}. $$

Using the guiding–center equations

$$ \dot{\mathbf{X}} = v_\parallel \mathbf{b} + \frac{\mu}{q|\mathbf{B}|}\,\mathbf{b}\times\nabla |\mathbf{B}|, \qquad \dot v_\parallel = -\mu\,\mathbf{b}\cdot\nabla |\mathbf{B}| - \mathbf{b}\cdot \nabla U, $$

we get

$$ \frac{dE}{dt} = \frac{\mu}{q|\mathbf{B}|}\, \nabla U \cdot \left( \mathbf{b}\times\nabla|\mathbf{B}| \right). $$ In the tokamak setup, $U = 0$ and $\nabla U = 0$, giving

$$ \frac{dE}{dt} = 0. $$

We can therefore conclude that the energy should be exactly conserved in the tokamak guiding center model, explaining the bounded energy error observed in the plot. Any small errors seen are likely due to numerical inaccuracies in the implementation of the ABM method. The fact that we got relatively normal position errors, but very small energy errors, further supports this conclusion.