lars/FYS4150
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Minor grammar corrections

Browse Source
This commit is contained in:
Lars Bogner
2025-10-22 14:19:55 +02:00
parent ef80ded820
commit 1dde2c2594
4 changed files with 23 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files
projects/project3/chapters/methods.tex
+14 -14
View File
@@ -1,5 +1,5 @@
\subsection*{Penning trap} \subsection*{Penning trap}
Penning traps use a combination of a magnetic and electric field to confine the trajectory of a charged particle to a finite volume over long periods of time\cite[p.~9-10]{vogel_particle_2024}. Finite volume in this context refers to volumes in the order of \unit{\mm \cubed} up to \unit{\cm \cubed}. Long time ranges imply a period that exceeds the period of motion of the particle by multiple orders. For the studied case we will use an ideal Penning trap, i.e. there are no inhomogenities in the fields. The magnetic field will be defined as Penning traps use a combination of a magnetic and electric field to confine the trajectory of a charged particle to a finite volume over long periods of time \cite[p.~9-10]{vogel_particle_2024}. Finite volume in this context refers to volumes in the order of \unit{\mm \cubed} up to \unit{\cm \cubed}. Long time ranges imply a period that exceeds the period of motion of the particle by multiple orders. For the studied case we will use an ideal Penning trap, i.e. there are no inhomogenities in the fields. The magnetic field will be defined as
\begin{equation} \label{eq:B_field} \begin{equation} \label{eq:B_field}
B(\vec r, t) = B_0 \Theta\left(d - |\vec r |\right) \cdot \vec e_z, B(\vec r, t) = B_0 \Theta\left(d - |\vec r |\right) \cdot \vec e_z,
\end{equation} \end{equation}
@@ -23,7 +23,7 @@ In the transversal plane, the equations of motion are coupled due to the Lorentz
\ddot x - \omega_0 \dot y - \frac{1}{2} \omega_z^2 x &\equiv \ddot x - \frac{q B_0}{m} \dot y - \frac{q V_0}{m d^2} x = 0, \label{eq:transv_eom_x}\\ \ddot x - \omega_0 \dot y - \frac{1}{2} \omega_z^2 x &\equiv \ddot x - \frac{q B_0}{m} \dot y - \frac{q V_0}{m d^2} x = 0, \label{eq:transv_eom_x}\\
\ddot y + \omega_0 \dot x - \frac{1}{2} \omega_z^2 y &\equiv \ddot y + \frac{q B_0}{m} \dot x - \frac{q V_0}{m d^2} y = 0 \label{eq:transv_eom_y}. \ddot y + \omega_0 \dot x - \frac{1}{2} \omega_z^2 y &\equiv \ddot y + \frac{q B_0}{m} \dot x - \frac{q V_0}{m d^2} y = 0 \label{eq:transv_eom_y}.
\end{align} \end{align}
Using the sum of equations: $\eqref{eq:transv_eom_x} + i \eqref{eq:transv_eom_y}$, we can use the general definition $f(t) \equiv x(t) + i y(t)$ to derive a single complex differential equation for the transversal motion Using the sum of equations: $\eqref{eq:transv_eom_x} + i \eqref{eq:transv_eom_y}$, we can utilize the general definition $f(t) \equiv x(t) + i y(t)$ to derive a single complex differential equation for the transversal motion
\begin{equation} \label{eq:complex_transv_eom} \begin{equation} \label{eq:complex_transv_eom}
\ddot f + i \omega_0 \dot f - \frac{1}{2} \omega_z^2 f = 0. \ddot f + i \omega_0 \dot f - \frac{1}{2} \omega_z^2 f = 0.
\end{equation} \end{equation}
@@ -40,7 +40,7 @@ This solution is valid for $\omega_0^2 > 2 \omega_z^2$. In the case of $\omega_0
\begin{equation} \begin{equation}
\frac{q}{m} B_0^2 > \frac{2 V_0}{d^2}. \frac{q}{m} B_0^2 > \frac{2 V_0}{d^2}.
\end{equation} \end{equation}
Given this condition, the solution for the transversal motion will be bound, i.e. $|f(t)| < \infty$ for all $t > 0$. Given the stability criterion, the two terms in \cref{eq:general_transv_solution} can be identified as sinusoidal motions with amplitudes $A_\pm$. If the two terms are in phase the upper limit of $|f| \equiv R_+$ will be $R_+ = A_+ + A_-$, while the lower limit for the case of antiparallel phases will be the absolute of the difference in the two amplitudes, $R_- = |A_+ - A_-|$. A analytical solution for a special case of initial conditions is derived in \cref{app:special_case_analytical_solution}. Given this condition, the solution for the transverse motion will be bounded, i.e. $|f(t)| < \infty$ for all $t > 0$. According to the stability criterion, the two terms in \cref{eq:general_transv_solution} can be identified as sinusoidal motions with amplitudes $A_\pm$. If the two terms are in phase the upper limit of $|f| \equiv R_+$ will be $R_+ = A_+ + A_-$, while the lower limit for the case of antiparallel phases will be the absolute of the difference in the two amplitudes, $R_- = |A_+ - A_-|$. A analytical solution for a special case of initial conditions is derived in \cref{app:special_case_analytical_solution}.
\subsubsection*{Multiple particles} \subsubsection*{Multiple particles}
For the case of multiple particles, it is important to include the particle-particle interactions. For this analysis we will only consider the Coulomb interaction between the particles, neglecting magnetic interactions. The force on particle $i$ due to all other particles $j$ is then given by For the case of multiple particles, it is important to include the particle-particle interactions. For this analysis we will only consider the Coulomb interaction between the particles, neglecting magnetic interactions. The force on particle $i$ due to all other particles $j$ is then given by
@@ -50,7 +50,7 @@ For the case of multiple particles, it is important to include the particle-part
\subsection*{Numerical Integration Methods} \subsection*{Numerical Integration Methods}
\subsubsection*{Forward Euler method} \subsubsection*{Forward Euler method}
The simplest numerical method to solve the ordinary differential equation system is the Forward Euler method. It is a first order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^2)$, with $h$ being the step size. The global error after $N$ steps is therefore on the order of $\mathcal{O}(h)$ \cite{fjordholm_numerical_nodate}. The method is explicit, meaning that the state of the system at the next time step can be calculated directly from the current state. Given a general ordinary differential equation of the form The simplest numerical method to solve the ordinary differential equation system is the Forward Euler method. It is a first order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^2)$, with $h$ being the step size. The global error after $N$ steps is therefore on the order of $\mathcal{O}(h)$ \cite{fjordholm_numerical_nodate}. This is an explicit method, meaning the state of the system at the next time step can be calculated directly from the current state. Given a general ordinary differential equation of the form
\begin{equation} \begin{equation}
\dot y(t) = f(t, y(t)), \dot y(t) = f(t, y(t)),
\end{equation} \end{equation}
@@ -63,7 +63,7 @@ For our second order differential equations, we first rewrite them as a system o
\vec r_{n+1} &= \vec r_n + h \vec v_n, \\ \vec r_{n+1} &= \vec r_n + h \vec v_n, \\
\vec v_{n+1} &= \vec v_n + h \frac{\vec F(\vec r_n, \vec v_n, t_n)}{m}. \vec v_{n+1} &= \vec v_n + h \frac{\vec F(\vec r_n, \vec v_n, t_n)}{m}.
\end{align} \end{align}
This method is computationally inexpensive, as we only need a single force evaluation per time step. However, it is known to be not very accurate. The computational expense of this method is minimal, as only a single force evaluation is required at each time step. However, it is acknowledged that this method is not particularly precise.
\subsubsection*{Runge-Kutta 4 method} \subsubsection*{Runge-Kutta 4 method}
The Runge-Kutta 4 (RK4) method is a popular and more accurate method for solving ordinary differential equations. It is a fourth-order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^5)$, and the global error after $N$ steps is on the order of $\mathcal{O}(h^4)$ \cite{fjordholm_numerical_nodate}. The RK4 method calculates the next state of the system using a weighted average of four different estimates of the slope (the derivative) at different points within the time step. Given a general ordinary differential equation of the form The Runge-Kutta 4 (RK4) method is a popular and more accurate method for solving ordinary differential equations. It is a fourth-order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^5)$, and the global error after $N$ steps is on the order of $\mathcal{O}(h^4)$ \cite{fjordholm_numerical_nodate}. The RK4 method calculates the next state of the system using a weighted average of four different estimates of the slope (the derivative) at different points within the time step. Given a general ordinary differential equation of the form
\begin{equation} \begin{equation}
@@ -77,7 +77,7 @@ the RK4 method updates the state of the system as follows \cite{cheever_fourth_2
k_4 &= h f(t_n + h, y_n + k_3), \\ k_4 &= h f(t_n + h, y_n + k_3), \\
y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4). y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4).
\end{align} \end{align}
As this method is again a generalized solver for ODE of the first order, we again rewrite our second order differential equations as a system of first order equations. As this method is, once more, a generalised solver for ODE of the first order, the second order differential equations are to be rewritten as a system of first order equations.
\subsubsection*{Velocity-Verlet method} \subsubsection*{Velocity-Verlet method}
In contrast to the previous two methods, the Velocity-Verlet method is a symplectic integrator which is specifically designed for second order differential equations of the form In contrast to the previous two methods, the Velocity-Verlet method is a symplectic integrator which is specifically designed for second order differential equations of the form
@@ -90,11 +90,11 @@ The method is time-reversible and conserves energy better over long time periods
\vec a_{n+1} &= \frac{\vec F(\vec r_{n+1}, t_{n+1})}{m}, \\ \vec a_{n+1} &= \frac{\vec F(\vec r_{n+1}, t_{n+1})}{m}, \\
\vec v_{n+1} &= \vec v_n + \frac{h}{2} (\vec a_n + \vec a_{n+1}), \vec v_{n+1} &= \vec v_n + \frac{h}{2} (\vec a_n + \vec a_{n+1}),
\end{align} \end{align}
It is important to note that the Velocity-Verlet integrator requires forces to be independent of velocity. In the case, that the algorithm is applicable, it is a very efficient method, as it only requires a single force evaluation per time step, while still being a second order method with a local truncation error per step on the order of $\mathcal{O}(h^3)$ and a global error after $N$ steps on the order of $\mathcal{O}(h^2)$\cite{kubar_lecture_2025}. It is important to note that the Velocity-Verlet integrator requires forces to be independent of velocity. In the case, that the algorithm is applicable, it is a very efficient method, as it only requires a single force evaluation per time step, while still being a second order method with a local truncation error per step on the order of $\mathcal{O}(h^3)$ and a global error after $N$ steps on the order of $\mathcal{O}(h^2)$ \cite{kubar_lecture_2025}.
Because of the symplectic nature, the low number of force evaluations and the good energy conservation properties, the Velocity-Verlet method is widely used and the quasi-standard in molecular dynamics simulations\cite{kubar_lecture_2025}. Because of the symplectic nature, the low number of force evaluations and the good energy conservation properties, the Velocity-Verlet method is widely used and the quasi-standard in molecular dynamics simulations \cite{kubar_lecture_2025}.
\subsubsection*{Boris algorithm} \subsubsection*{Boris algorithm}
The Boris algorithm is a widely used method for integrating the equations of motion of charged particles in electromagnetic fields \cite{webb_symplectic_2014,noauthor_pdf_2025}. It is particularly well-suited for problems where the Lorentz force plays a significant role, as it is designed to handle the velocity-dependent nature of the magnetic force. The Boris algorithm is a symplectic integrator\cite{webb_symplectic_2014}\footnote{The simplecticity of Boris algorithm is controversial with disagreement between papers \cite{noauthor_pdf_2025}. All publications agree on the fact, that Boris algorithm is phase-space preserving at the very least.}, which means it conserves the phase space volume and exhibits good long-term energy conservation properties. The algorithm updates the position and velocity of a charged particle in a magnetic field as follows: The Boris algorithm is a widely utilised technique for integrating the equations of motion of charged particles in electromagnetic fields \cite{webb_symplectic_2014,noauthor_pdf_2025}. The device is especially well-suited to problems in which the Lorentz force is significant, as it is designed to handle the velocity-dependent nature of the magnetic force. The Boris algorithm is a symplectic integrator \cite{webb_symplectic_2014}.\footnote{However, there is a lack of consensus on the symplecticity of the Boris algorithm, with different papers providing contradictory evidence \cite{noauthor_pdf_2025}. It is an irrefutable conclusion that all publications concur on the fact that the Boris algorithm is phase-space preserving at the very least.} This signifies that it conserves the phase space volume and exhibits favourable long-term energy conservation properties. The algorithm responsible for this calculation updates the position and velocity of a charged particle in a magnetic field as follows:
We define two auxiliary velocities $\vec v^-$ and $\vec v^+$, which represent the velocity before and after the magnetic field rotation, respectively. The algorithm proceeds in the following steps \cite{hoppock_iwhoppockboris-algorithm_2025}: We define two auxiliary velocities $\vec v^-$ and $\vec v^+$, which represent the velocity before and after the magnetic field rotation, respectively. The algorithm proceeds in the following steps \cite{hoppock_iwhoppockboris-algorithm_2025}:
\begin{align} \begin{align}
\vec v^- &= \vec v_n + \frac{q \vec E h}{2m}\\ \vec v^- &= \vec v_n + \frac{q \vec E h}{2m}\\
@@ -116,7 +116,7 @@ The final update of the velocity and position is then given by
The basic framework for the numerical analysis is based on a \texttt{PenningTrap} class, which contains all the particles present in the trap, as well as a parametrization of the electric and magnetic fields. The particles are represented by a \texttt{Particle} class, which contains the physical properties of the particle, as well as its current position and velocity. With this information, the \texttt{PenningTrap} class can calculate the forces acting on each particle, including the external fields and the particle-particle interactions. The particle-particle interactions can be toggled on and off, allowing for a simulation of both scenarios. The external fields can be modified by supplying a field-method of the form \texttt{external\_field(const arma::vec\& r, double t, const PenningTrap\& trap)}. The reference to the \texttt{PenningTrap} allows for the parameters of the field to be stored in the trap object. The implementation of the \texttt{Particle} and \texttt{PenningTrap} classes can be found in \texttt{/src/project3/include/classes.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/classes.cpp}. The basic framework for the numerical analysis is based on a \texttt{PenningTrap} class, which contains all the particles present in the trap, as well as a parametrization of the electric and magnetic fields. The particles are represented by a \texttt{Particle} class, which contains the physical properties of the particle, as well as its current position and velocity. With this information, the \texttt{PenningTrap} class can calculate the forces acting on each particle, including the external fields and the particle-particle interactions. The particle-particle interactions can be toggled on and off, allowing for a simulation of both scenarios. The external fields can be modified by supplying a field-method of the form \texttt{external\_field(const arma::vec\& r, double t, const PenningTrap\& trap)}. The reference to the \texttt{PenningTrap} allows for the parameters of the field to be stored in the trap object. The implementation of the \texttt{Particle} and \texttt{PenningTrap} classes can be found in \texttt{/src/project3/include/classes.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/classes.cpp}.
\subsection*{Numerical Methods Implementation} \subsection*{Numerical Methods Implementation}
All numerical methods are implemented as classes inheriting from a base class \texttt{Solver}. The base class contains a reference to the \texttt{PenningTrap} object, as well as the time step size. The recording of particle properties over time, like position and velocity is part of the general \texttt{Solver} class. Each derived class implements the \texttt{step()} method, which updates the state of the system by one time step using the respective numerical method. The implementation of the \texttt{Solver} class and its derived classes can be found in \texttt{/src/project3/include/solvers.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/solvers.cpp}. The following solvers are implemented:\texttt{Forward\-Euler\-Solver}, \texttt{RK4\-Solver}, \texttt{Velocity\-Verlet\-Solver} and \texttt{Analytical\-Solver}, which implements the analytical solution for a special case of initial conditions (see \cref{app:special_case_analytical_solution}). All numerical methods are implemented as classes inheriting from a base class \texttt{Solver}. The base class contains a reference to the \texttt{PenningTrap} object, as well as the time step size. The recording of particle properties over time, like position and velocity is part of the general \texttt{Solver} class. Each derived class implements the \texttt{step()} method, which updates the state of the system by one time step using the respective numerical method. The implementation of the \texttt{Solver} class and its derived classes can be found in \texttt{/src/project3/include/solvers.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/solvers.cpp}. The following solvers are implemented: \texttt{Forward\-Euler\-Solver}, \texttt{RK4\-Solver}, \texttt{Velocity\-Verlet\-Solver} and \texttt{Analytical\-Solver}, which implements the analytical solution for a special case of initial conditions (see \cref{app:special_case_analytical_solution}).
Unless otherwise specified, the simulations are run with the following parameters: Unless otherwise specified, the simulations are run with the following parameters:
\begin{itemize} \begin{itemize}
@@ -125,9 +125,9 @@ Unless otherwise specified, the simulations are run with the following parameter
\item Simulation parameters: $T = n_\mathrm{steps} \cdot h = \qty{50}{\us}$. \item Simulation parameters: $T = n_\mathrm{steps} \cdot h = \qty{50}{\us}$.
\end{itemize} \end{itemize}
\subsection*{Tools} \subsection*{Tools}
All simulations are written in C++ and compiled using \texttt{g++ (GCC) 15.2.1 20250808 (Red Hat 15.2.1-1)}. For optimization, the \texttt{-O3} flag is used. The code makes use of the Armadillo library\cite{sanderson_armadillo_2025,sanderson_practical_2019} for linear algebra operations and the argparse library\cite{pranav_p-ranavargparse_2025} for command line argument parsing. All simulations are written in C++ and compiled using \texttt{g++ (GCC) 15.2.1 20250808 (Red Hat 15.2.1-1)}. For optimization, the \texttt{-O3} flag is used. The code makes use of the Armadillo library \cite{sanderson_armadillo_2025,sanderson_practical_2019} for linear algebra operations and the argparse library \cite{pranav_p-ranavargparse_2025} for command line argument parsing.
To visualize the results and analyze the data, Python 3.13.7 with the libraries \texttt{numpy} \cite{harris_array_2020}, \texttt{matplotlib} \cite{hunter_matplotlib_2007} and \texttt{pandas} \cite{team_pandas-devpandas_2025} is used. To create easy to use command line interfaces for the Python scripts, the \texttt{typer} library\cite{ramirez_typer_nodate} is used. To visualize the results and analyze the data, Python 3.13.7 with the libraries \texttt{numpy} \cite{harris_array_2020}, \texttt{matplotlib} \cite{hunter_matplotlib_2007} and \texttt{pandas} \cite{team_pandas-devpandas_2025} is used. To create easy to use command line interfaces for the Python scripts, the \texttt{typer} library \cite{ramirez_typer_nodate} is used.
As a large number of simulations need to be run for different parameters, a bash script is used to automate the process. The script can be found in \texttt{/src/project3/generate\_results.sh} of the project repository. It uses the \texttt{pueue} task management tool\cite{beer_nukesorpueue_2025} to run multiple simulations in parallel, making use of all available CPU cores. To automate the plot generation, the script \texttt{generate\_plots.sh} is used, which can be found in the same directory. As a large number of simulations need to be run for different parameters, a bash script is used to automate the process. The script can be found in \texttt{/src/project3/generate\_results.sh} of the project repository. It uses the \texttt{pueue} task management tool \cite{beer_nukesorpueue_2025} to run multiple simulations in parallel, making use of all available CPU cores. To automate the plot generation, the script \texttt{generate\_plots.sh} is used, which can be found in the same directory.
To aid in the development process, the large language model GitHub Copilot\cite{noauthor_github_2025} was used to generate code snippets and provide suggestions for code completion. Other large language models were not used in the writing of this report. To aid in the development process, the large language model GitHub Copilot \cite{noauthor_github_2025} was used to generate code snippets and provide suggestions for code completion. Other large language models were not used in the writing of this report.
projects/project3/chapters/results.tex
+8 -8
View File
@@ -5,7 +5,7 @@ As there is no closed form solution to the trajectories in a Penning trap contai
\end{equation} \end{equation}
is used, where $\vec{r}_{\text{exact}}$ as derived in \cref{sec:methods}, and $\vec{r}_{\text{numerical}}$ is the numerical solution obtained from the implemented algorithms. The relative error is calculated at each time step, and the maximum value over the entire simulation is reported as the error for that particular simulation. is used, where $\vec{r}_{\text{exact}}$ as derived in \cref{sec:methods}, and $\vec{r}_{\text{numerical}}$ is the numerical solution obtained from the implemented algorithms. The relative error is calculated at each time step, and the maximum value over the entire simulation is reported as the error for that particular simulation.
Apart from the quantitative measure of the accuracy, the qualitative behavior of the trajectories is also examined. The numerical solutions are plotted and compared to the expected behavior. Especially properties like energy drift of the solution and the general shape of the trajectory can give insight into the performance of the algorithms. Apart from the quantitative measure of the accuracy, the qualitative behavior of the trajectories is also examined. The numerical solutions are plotted and compared to the expected behavior. Especially properties like the energy drift of the solution and the general shape of the trajectory can give insight into the performance of the algorithms.
\begin{figure} \begin{figure}
\centering \centering
@@ -13,13 +13,13 @@ Apart from the quantitative measure of the accuracy, the qualitative behavior of
\caption{Trajectory of two particles in a Penning trap, simulated using the RK4 method. The particles are initialized as: $\vec r_1 = (\qty{20}{\um}, \qty{0}{um}, \qty{20}{\um})$, $\vec v_1 = (\qty{0}{\um/\us}, \qty{25}{\um/\us}, \qty{0}{\um/\us})$, $\vec r_2 = (\qty{25}{\um}, \qty{25}{\um}, \qty{0}{\um})$, $\vec v_2 = (\qty{0}{\um/\us}, \qty{40}{\um/\us}, \qty{5}{\um/\us})$.} \caption{Trajectory of two particles in a Penning trap, simulated using the RK4 method. The particles are initialized as: $\vec r_1 = (\qty{20}{\um}, \qty{0}{um}, \qty{20}{\um})$, $\vec v_1 = (\qty{0}{\um/\us}, \qty{25}{\um/\us}, \qty{0}{\um/\us})$, $\vec r_2 = (\qty{25}{\um}, \qty{25}{\um}, \qty{0}{\um})$, $\vec v_2 = (\qty{0}{\um/\us}, \qty{40}{\um/\us}, \qty{5}{\um/\us})$.}
\label{fig:two_particle_trajectory_rk4} \label{fig:two_particle_trajectory_rk4}
\end{figure} \end{figure}
The Runge-Kutta 4 method shows a very good qualitative behavior, with no apparent inaccuracies in the trajectory, as shown in \cref{fig:two_particle_trajectory_rk4}. Comparing against the same trajectories calculated using the Velocity-Verlet method, we can immediately see shortcomings of the latter. Trajectories for the Velocity-Verlet algorithm, the Boris algorithm and Euler algorithm are included in \cref{fig:two_particle_trajectory_verlet,fig:two_particle_trajectory_boris,fig:two_particle_trajectory_euler} in \cref{app:further_results}. The trajectory shows a significant drift over time, which is not present in the RK4 solution. The Verlet solution rapidly increases its radius, which is not physical in a system with conserved energy. While one would expect the symplectic Verlet method to perform better in this regard, the velocity-dependant Lorentz force is not suited for computations using the Velocity-Verlet algorithm, as discussed in \cref{sec:methods}. The Runge-Kutta 4 method shows a very good qualitative behavior, with no apparent inaccuracies in the trajectory, as shown in \cref{fig:two_particle_trajectory_rk4}. A comparison of the same trajectories calculated using the Velocity-Verlet method reveals the latter's shortcomings. Trajectories for the Velocity-Verlet algorithm, the Boris algorithm and Euler algorithm are included in \cref{fig:two_particle_trajectory_verlet,fig:two_particle_trajectory_boris,fig:two_particle_trajectory_euler} in \cref{app:further_results}. The trajectory shows a significant drift over time, which is not present in the RK4 solution. The Verlet solution rapidly increases its radius, which is not physical in a system with conserved energy. While one would expect the symplectic Verlet method to perform better in this regard, the velocity-dependent Lorentz force is not suited for computations using the Velocity-Verlet algorithm, as discussed in \cref{sec:methods}.
A similar behavior of non-physical qualitative behavior is observed for the Euler method, which also shows a significant energy drift. The Boris algorithm, which is specifically designed to handle electromagnetic forces, shows a qualitatively correct behavior, with no apparent energy drift and a stable orbit over time. A similar behavior of non-physical qualitative behavior is observed for the Euler method, which also shows a significant energy drift. The Boris algorithm, which is specifically designed to handle electromagnetic forces, shows a qualitatively correct behavior, with no apparent energy drift and a stable orbit over time.
As well as the trajectory in the transversal plane, the osciallatory behavior in the $z$ direction and the phase space usage in all three dimensions is examined. As there is no influence of the magnetic field in the $z$ direction, the oscillations should be perfectly harmonic. Furthermore should the Velocity-Verlet method exhibit correct behavior, as there is no velocity-dependant force in this direction. The very good performance of RK4, Boris and Velocity-Verlet in this regard is confirmed, while the Euler method shows energy drift and a non-harmonic behavior with increasing amplitude over time. As well as the trajectory in the transversal plane, the oscillatory behavior in the $z$ direction and the phase space usage in all three dimensions is examined. As there is no influence of the magnetic field in the $z$ direction, the oscillations should be perfectly harmonic. Furthermore, should the Velocity-Verlet method exhibit correct behavior, as there is no velocity-dependent force in this direction. The very good performance of RK4, Boris and Velocity-Verlet in this regard is confirmed, while the Euler method shows energy drift and a non-harmonic behavior with increasing amplitude over time.
The numerical analysis of the trajectories shows similar behavior for the two algorithms designed to handle electromagnetic forces, RK4 and Boris. The evolution of the relative error $\epsilon$ as defined in \cref{eq:rel_error} over time is shown in \cref{fig:relative_error}. Both algorithms show a steady osciallation in the relative error around \numrange{1e-3}{2e-3}, with RK4 showing a slightly lower error overall. The Velocity-Verlet and Euler methods show a significantly higher error, with the Euler method reaching a relative error of up to \num{0.1} at certain times. The Velocity-Verlet method shows a relative error of up to \num{0.2}. The numerical analysis of the trajectories shows similar behavior for the two algorithms designed to handle electromagnetic forces, RK4 and Boris. The evolution of the relative error $\epsilon$ as defined in \cref{eq:rel_error} over time is shown in \cref{fig:relative_error}. Both algorithms show a steady oscillation in the relative error around \numrange{1e-3}{2e-3}, with RK4 showing a slightly lower error overall. The Velocity-Verlet and Euler methods show a significantly higher error, with the Euler method reaching a relative error of up to \num{0.1} at certain times. The Velocity-Verlet method shows a relative error of up to \num{0.2}.
\begin{figure} \begin{figure}
\centering \centering
@@ -29,7 +29,7 @@ The numerical analysis of the trajectories shows similar behavior for the two al
\end{figure} \end{figure}
\subsection*{Many-Body Simulations} \subsection*{Many-Body Simulations}
In addition to the case of perfect trajectories without Coulomb interactions, the behavior of a cloud of up to \num{100} particles is examined. This use case is of particular interest, as it is closer to real-world applications of Penning traps, such as in mass spectrometry. The limit of increasing particle number is also of interest, as the computational complexity increases significantly with the number of particles, due to the expected $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. Additionally, the collective behavior of particles, i.e. the statistical properties of the ensemble, and the accuracy thereof is studied. In addition to the case of perfect trajectories in the absence of Coulomb interactions, the behavior of a cloud of up to \num{100} particles is examined. This use case is of particular interest, as it is closer to real-world applications of Penning traps, such as in mass spectrometry. The limit of the increasing particle number is also a subject of interest, as the computational complexity increases significantly with the number of particles. This is due to the expected $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. Furthermore, the collective behavior of particles, i.e. the statistical properties of the ensemble, and the accuracy thereof, is studied.
All particles are initialized with random positions and velocities, sampled from a normal distribution. The standard deviation of the position distribution is set to \qty{50}{\um} in all three dimensions, while the standard deviation of the velocity distribution is set to \qty{50}{\um/\us}. To compare the influence of Coulomb interactions on the simulation performance and statistical properties, simulations with and without interactions are performed. All particles are initialized with random positions and velocities, sampled from a normal distribution. The standard deviation of the position distribution is set to \qty{50}{\um} in all three dimensions, while the standard deviation of the velocity distribution is set to \qty{50}{\um/\us}. To compare the influence of Coulomb interactions on the simulation performance and statistical properties, simulations with and without interactions are performed.
@@ -65,7 +65,7 @@ A further behavior that can be studied in the many-body simulations is the reson
\end{tabular} \end{tabular}
\end{table} \end{table}
For different step sizes, i.e. a different number of total time steps, the total run time of a simulation with different numbers of particles is measured. For measuring the run time, the \texttt{ctime} module with the \texttt{std::clock()} function is used. All runs were performed on the same machine, with an AMD Ryzen 7 4750U CPU and \qty{16}{\giga\byte} of RAM, running Fedora 42. For parallelization of the runs the scheduler \texttt{pueue} was used with \num{14} simultaneous jobs. To reduce the influence of other processes on the machine, the runs were performed at night when the machine was otherwise idle. Furthermore, the runs were interweaved in a way, that no solving algorithm could take dvantage from time dependant properties of the machine, such as thermal throttling. The output during the simulation was enabled, as the output is reused for the analysis of the statistical properties of the particle ensemble. The I/O operations were performed outside of the timed section of the code, to get a more accurate measure of the run time of the algorithms themselves. The times were averaged over \num{4} runs with different time steps, and the relative time per particle and time step was calculated. The results are shown in \cref{tab:solver_scaling}. For different step sizes, i.e. a different number of total time steps, the total run time of a simulation with different numbers of particles is measured. For measuring the run time, the \texttt{ctime} module with the \texttt{std::clock()} function is used. All runs were performed on the same machine, with an AMD Ryzen 7 4750U CPU and \qty{16}{\giga\byte} of RAM, running Fedora 42. For parallelization of the runs the scheduler \texttt{pueue} was used with \num{14} simultaneous jobs. To reduce the influence of other processes on the machine, the runs were performed at night when the machine was otherwise idle. Furthermore, the runs were interwoven in a way, that no solving algorithm could take advantage from time dependent properties of the machine, such as thermal throttling. The output during the simulation was enabled, as the output is reused for the analysis of the statistical properties of the particle ensemble. The I/O operations were performed outside the timed section of the code, to get a more accurate measure of the run time of the algorithms themselves. The times were averaged over \num{4} runs with different time steps, and the relative time per particle and time step was calculated. The results are shown in \cref{tab:solver_scaling}.
The results show a clear linear scaling of the relative time per particle and time step with increasing number of particles, as expected from the $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. The Euler, Velocity-Verlet and Boris methods show a very similar performance, with the RK4 method being slower by a factor of approximately \num{2.8}. The similar performance of the Euler, Velocity-Verlet and Boris methods is expected, as all three methods require a single force calculation per time step. The RK4 method requires four force calculations per time step, which explains the significantly higher run time. From the scaling behavior for simulations with Coulomb interactions, we can derive that the major contributing factor to the run time is the calculation of the interactions. The results show a clear linear scaling of the relative time per particle and time step with increasing number of particles, as expected from the $\mathcal{O}(N^2)$ scaling of the Coulomb interaction calculations. The Euler, Velocity-Verlet and Boris methods show a very similar performance, with the RK4 method being slower by a factor of approximately \num{2.8}. The similar performance of the Euler, Velocity-Verlet and Boris methods is expected, as all three methods require a single force calculation per time step. The RK4 method requires four force calculations per time step, which explains the significantly higher run time. From the scaling behavior for simulations with Coulomb interactions, we can derive that the major contributing factor to the run time is the calculation of the interactions.
@@ -100,7 +100,7 @@ A qualitative analysis of the symplectic properties is performed by examining th
\end{equation} \end{equation}
over the simulation duration. As the system is isolated, the total energy should be conserved. A non-conservation of the total energy indicates a non-symplectic behavior of the algorithm. over the simulation duration. As the system is isolated, the total energy should be conserved. A non-conservation of the total energy indicates a non-symplectic behavior of the algorithm.
\Cref{fig:total_energy} shows the total energy development of the system for the different algorithms. A more detailed analysis on the dependance on step size is shown in \cref{app:further_results} with \cref{fig:boris_many_particles_energy,fig:euler_many_particles_energy,fig:verlet_many_particles_energy,fig:rk4_many_particles_energy}. The RK4 and Boris methods show a very stable total energy over the entire simulation duration, with slight oscillations around a mean value. The Velocity-Verlet method shows a significant drift to increasing total energy over time, indicating a non-symplectic behavior in this system. The Euler method shows a similar behavior, with a slightly less pronounced increase in total energy over time. \Cref{fig:total_energy} shows the total energy development of the system for the different algorithms. A more detailed analysis on the dependence on step size is shown in \cref{app:further_results} with \cref{fig:boris_many_particles_energy,fig:euler_many_particles_energy,fig:verlet_many_particles_energy,fig:rk4_many_particles_energy}. The RK4 and Boris methods show a very stable total energy over the entire simulation duration, with slight oscillations around a mean value. The Velocity-Verlet method shows a significant drift to increasing total energy over time, indicating a non-symplectic behavior in this system. The Euler method shows a similar behavior, with a slightly less pronounced increase in total energy over time.
\begin{figure} \begin{figure}
\centering \centering
\includegraphics[width=\columnwidth]{../../src/project3/python/plots/many_particles_100_particles_32000_steps_enabled_interactions_total_energies.pdf} \includegraphics[width=\columnwidth]{../../src/project3/python/plots/many_particles_100_particles_32000_steps_enabled_interactions_total_energies.pdf}
@@ -108,6 +108,6 @@ over the simulation duration. As the system is isolated, the total energy should
\label{fig:total_energy} \label{fig:total_energy}
\end{figure} \end{figure}
The non-symplectic behavior of the Velocity-Verlet method in this system is expected due to the velocity-dependant Lorentz force, as discussed in \cref{sec:methods}. The non-symplectic behavior of the Euler method is also expected, as it is not a symplectic algorithm. The good energy conservation of the RK4 method is noteworthy, as it is not a symplectic algorithm either. However, a conservation of total energy in a single simulation for a specific time span does not imply symplectic properties of the algorithm. In general RK4 is not considered a symplectic algorithm, and may show non-symplectic behavior in other systems or over longer time spans. The non-symplectic behavior of the Velocity-Verlet method in this system is expected due to the velocity-dependent Lorentz force, as discussed in \cref{sec:methods}. The non-symplectic behavior of the Euler method is also expected, as it is not a symplectic algorithm. The good energy conservation of the RK4 method is noteworthy, as it is not a symplectic algorithm either. However, a conservation of total energy in a single simulation for a specific time span does not imply symplectic properties of the algorithm. In general RK4 is not considered a symplectic algorithm, and may show non-symplectic behavior in other systems or over longer time spans.
Boris method also shows a very good conservation of total energy, which is expected due to its design for electromagnetic systems and symplectic properties therein. Boris method also shows a very good conservation of total energy, which is expected due to its design for electromagnetic systems and symplectic properties therein.
projects/project3/main.pdf
BIN
View File
Binary file not shown.
projects/project3/main.tex
+1 -1
View File
@@ -97,7 +97,7 @@
%This is how we create an abstract section. %This is how we create an abstract section.
\begin{abstract} \begin{abstract}
Using integration algorithms for ordinary differential equations, we predict the time evolution of charged particles in the electromagnetic fields of a Penning trap. Comparing the computational cost and numerical accuracy of the Forward Euler, Velocity-Verlet, Boris and Runge-Kutta 4th order methods, we find that the Boris algorithm is the most efficient choice for this problem. It is phase-space preserving, energy conserving and reduces the computational cost by up to \qty{700}{\percent} compared to general purpose approaches. To study ensemble properties of the trapped particles, we simulate the dynamics of 100 particles to find resonances and observe energy conservation of the different algorithms. Comparing to analytical solutions, at \num{32000} time steps of $h = \qty{1.56}{\ns}$, the Boris algorithm achieves a relative numerical error of the trajectory of less than \num{2e-3}. Using integration algorithms for ordinary differential equations, we predict the time evolution of charged particles in the electromagnetic fields of a Penning trap. Comparing the computational cost and the numerical accuracy of the Forward Euler, Velocity-Verlet, Boris and Runge-Kutta 4th order methods, we find that the Boris algorithm is the most efficient choice for this problem. It is phase-space preserving, energy conserving and reduces the computational cost by up to \qty{700}{\percent} compared to general purpose approaches. To study ensemble properties of the trapped particles, we simulate the dynamics of 100 particles to find resonances and observe energy conservation of the different algorithms. Comparing to analytical solutions, at \num{32000} time steps of $h = \qty{1.56}{\ns}$, the Boris algorithm achieves a relative numerical error of the trajectory of less than \num{2e-3}.
\end{abstract} \end{abstract}
\maketitle \maketitle