Project-3/latex/sections/methods.tex

\documentclass[../main.tex]{subfiles}
\resources % search for graphics in ../res/imgs/


\begin{document}

\section{Methods}

\subsection{Units and constants}

Before continuing, we need to define the units we'll be working with.
Since we are working with particles, we need small units to work with so the
numbers we are working with aren't so small that they could potentially lead
to large round-off errors in our simulation. The units that we will use are listed in Table~\ref{tab:units}.
\begin{table}[H]
    \begin{center}
        \begin{tabular}[c]{lll}
            Dimension & Unit & Symbol \\
            \hline
            Length & micrometer & $\mu m$ \\
            Time & microseconds & $\mu s$ \\
            Mass & atomic mass unit & $u$ \\
            Charge & the elementary charge & $e$ \\
            \hline
        \end{tabular}
    \end{center}
    \caption{The set of units we'll be working with.}\label{tab:units}
\end{table}
With these base units, we get
\begin{equation}
    k_e = 1.3893533 \cdot 10^5 \frac{u(\mu m)^3}{(\mu s)^2 e},
\end{equation}
and we get that the unit for magnetic field strength (Tesla, $T$) and electric potential (Volt, $V$) are
\begin{equation}
    \begin{split}
        T &= 9.64852558 \cdot 10^1 \frac{u}{(\mu s) e} \\
        V &= 9.64852558 \cdot 10^7 \frac{u (\mu m)^2}{(\mu s)^2 e}. \\
    \end{split}
    \label{eq:}
\end{equation}

\subsection{Dealing with a multi--particle system}

For a multi-particles system, we need to modify $\vb{F}$ to account for the 
force of other particles in the system acting upon each other. To do that, we
add another term to $\vb{F}$
\begin{equation}
    \vb{F}_i(t, \vb{v}_i, \vb{r}_i) = q_i \vb{E}(t, \vb{r}_i) + q_i \vb{v}_i \cross \vb{B} - \vb{E}_p(t, \vb{r}_i),
\end{equation}
where $i$ and $j$ are particle indices and
\begin{equation}
    \vb{E}_p(t, \vb{r}_i) = q_i k_e \sum_{j \neq i} 
    q_j \frac{\vb{r_i} - \vb{r_j}}{\left| \vb{r_i} - \vb{r_j} \right|^3}.
    \label{eq:}
\end{equation}

Newton's second law for a particle $i$ is then
\begin{equation}
    \frac{d^2\vb{r}_i}{dt^2} = \frac{\vb{F}_i\left(t, \frac{d\vb{r}_i}{dt}, \vb{r_i}\right)}{m_i},
    \label{eq:newtonlaw2}
\end{equation}
We can then rewrite the second order ODE from equation~\ref{eq:newtonlaw2} 
into a set of coupled first order ODEs. 
We now rewrite Newton's second law of motion as
\begin{equation}
    \begin{split}
        \frac{d\vb{r}_i}{dt} &= \vb{v}_i \\
        \frac{d\vb{v}_i}{dt} &= \frac{\vb{F}_i(t, \vb{v}_i, \vb{r}_i)}{m_i}.
    \end{split}
    \label{eq:coupled}
\end{equation}

\subsection{Forward Euler}

For a particle $i$, the forward Euler method for a coupled system is 
expressed as
\begin{equation}
    \begin{split}
        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + h \cdot \frac{d\vb{r}_{i,j}}{dt} = \vb{r}_{i,j} + h \cdot \vb{v}_{i,j} \\
        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + h \cdot \frac{\vb{v}_{i,j}}{dt} = \vb{v}_{i,j} + h \cdot \frac{\vb{F}\left( t_{j}, \vb{v}_{i,j}, \vb{r}_{i,j} \right)}{m},
    \end{split}
\end{equation}
for particle $i$ where $j$ is the current time step of the particle, 
$m$ is the mass of the particle, and $h$ is the step length.

When dealing with a multi-particle system, we need to ensure that we do not
update the position of any particles until every particle has calculated their
next step. An easy way of doing this is to create a copy of all the particles,
then update the copy, and when all the particles have calculated their next 
step, simply replace the particles with the copies. 
Algorithm~\ref{algo:forwardeuler} provides an overview on how that can be achieved. % Make this better

\begin{figure}[H]
    \begin{algorithm}[H]
    \caption{Forward Euler method}
    \label{algo:forwardeuler}
        \begin{algorithmic}
            \Procedure{Evolve forward Euler}{$particles, dt$}
            \State $N \leftarrow \text{Number of particles in } particles$
            \State $a \leftarrow \text{Calculate } \frac{\vb{F_i}}{m_i} \text{ for each particle in } particles$
            \For{ $i = 1, 2, \ldots , N$ }
            \State $particles_i.\vb{r} \leftarrow particles_i.\vb{r} + dt \cdot particles_i.\vb{v}$
            \State $particles_i.\vb{v} \leftarrow particles_i.\vb{v} + dt \cdot a_i$
            \EndFor
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
\end{figure}

\subsection{4th order Runge-Kutta}

For a particle $i$, we can express the 4th order Runge-Kutta (RK4) method as
\begin{equation}
    \begin{split}
        \vb{v}_{i,j+1} &= \vb{v}_{i,j} + \frac{h}{6} \left( \vb{k}_{\vb{v},1,i} 
            + 2\vb{k}_{\vb{v},2,i} + 2\vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} 
            \right) \\
        \vb{r}_{i,j+1} &= \vb{r}_{i,j} + \frac{h}{6} \left( \vb{k}_{\vb{r},1,i} 
            + 2\vb{k}_{\vb{r},2,i} + 2\vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} 
            \right),
    \end{split}
\end{equation}
where
\begin{equation}
    \begin{split}
        \vb{k}_{\vb{v},1,i} &= \frac{\vb{F}_i(t_j, \vb{v}_{i,j}, 
            \vb{r}_{i,j})}{m} \\
        \vb{k}_{\vb{r},1,i} &= \vb{v}_{i,j} \\
        \vb{k}_{\vb{v},2,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}, \vb{r}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{r},1,i}}{2})}{m} \\
        \vb{k}_{\vb{r},2,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2} \\
        \vb{k}_{\vb{v},3,i} &= \frac{\vb{F}_i(t_j+\frac{h}{2}, \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2}, \vb{r}_{i,j} 
            + h \frac{\cdot \vb{k}_{\vb{r},2,i}}{2})}{m} \\
        \vb{k}_{\vb{r},3,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},2,i}}{2} \\
        \vb{k}_{\vb{v},4,i} &= \frac{\vb{F}_i(t_j+h, \vb{v}_{i,j} 
            + h \cdot \vb{k}_{\vb{v},3,i}, \vb{r}_{i,j} 
            + h \cdot \vb{k}_{\vb{r},3,i})}{m} \\
        \vb{k}_{\vb{r},4,i} &= \vb{v}_{i,j} 
            + h \cdot \frac{\vb{k}_{\vb{v},1,i}}{2}.
    \end{split}
\end{equation}

In order to find each $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$, 
we need to first compute all $\vb{k}_{\vb{r},i}$ and $\vb{k}_{\vb{v},i}$ 
for all particles, then we can update the particles in order to compute 
$\vb{k}_{\vb{r},i+1}$ and $\vb{k}_{\vb{v},i+1}$. In order for the algorithm
to work, we need to save a copy of each particle before starting so that we
can update the particles correctly for each step.

This approach would require 8 loops to be able to complete the calculation since
we cannot update the particles until after all $\vb{k}$ values have been 
computed, however if we create a temporary array that holds particles, we can
put the updated particles in there, and then use that array in the next loop,
and would reduce the required amount of loops down to 4.

\begin{figure}
    \begin{algorithm}[H]
    \caption{RK4 method}
    \label{algo:rk4}
        \begin{algorithmic}
            \Procedure{Evolve RK4}{$particles, dt$}
            \State $N \leftarrow \text{Number of particles inside the Penning trap}$
            \State $orig\_p \leftarrow \text{Copy of particles}$ 
            \State $tmp\_p \leftarrow \text{Array of particles of size }N$
            \State $\vb{k}_{\vb{r}} \leftarrow \text{2D array of vectors of size } 4 \cross N$
            \State $\vb{k}_{\vb{v}} \leftarrow \text{2D array of vectors of size } 4 \cross N$

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},1,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},1,i} \leftarrow \frac{\vb{F}_i}{m_i}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{r},1,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{v},1,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},2,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},2,i} \leftarrow \frac{\vb{F}_i}{m_i}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{r},2,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} 
                + \frac{dt}{2} \cdot \vb{k}_{\vb{v},2,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},3,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},3,i} \leftarrow \frac{\vb{F}_i}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + dt \cdot \vb{k}_{\vb{r},3,i}$
            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + dt \cdot \vb{k}_{\vb{v},3,i}$
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Update particles}

            \For{ $i = 1, 2, \ldots, N$ }
            \State $\vb{k}_{\vb{r},4,i} \leftarrow particles_i.\vb{v}$
            \State $\vb{k}_{\vb{v},4,i} \leftarrow \frac{\vb{F}}{m}$

            \State $tmp\_p_i.\vb{r} \leftarrow orig\_p_i.\vb{r} + \frac{dt}{6} 
                \cdot \left( \vb{k}_{\vb{r},1,i} + \vb{k}_{\vb{r},2,i} 
                + \vb{k}_{\vb{r},3,i} + \vb{k}_{\vb{r},4,i} \right)$

            \State $tmp\_p_i.\vb{v} \leftarrow orig\_p_i.\vb{v} + \frac{dt}{6} 
                \cdot \left( \vb{k}_{\vb{v},1,i} + \vb{k}_{\vb{v},2,i} 
                + \vb{k}_{\vb{v},3,i} + \vb{k}_{\vb{v},4,i} \right)$
            
            \EndFor

            \State $particles \leftarrow tmp\_p$ \Comment{Final update}
            \EndProcedure
        \end{algorithmic}
    \end{algorithm}
    $\vb{F}$ in the algorithm does not take any arguments as it uses the 
    velocities and positions of the particles inside the array $particles$ to
    calculate the total force acting on particle $i$.
\end{figure}


\subsection{Testing the simulation}

\subsection{Relative error and error convergance rate}

\biblio
\end{document}