diff --git a/latex/main.pdf b/latex/main.pdf index d1fe6a2..c6f27db 100644 Binary files a/latex/main.pdf and b/latex/main.pdf differ diff --git a/latex/references/references.bib b/latex/references/references.bib index b8d6413..6a49b46 100644 --- a/latex/references/references.bib +++ b/latex/references/references.bib @@ -74,4 +74,14 @@ publisher = {Farleia Forlag}, year = {2018}, pages = {162--163} -} \ No newline at end of file +} + + +@article{rk4_method, + author = "Morten Hjorth-Jensen", + title = "Computational Physics, Lecture Notes Fall 2015", + journal = "Department of Physics, University of Oslo", + year = "2015", + chapter = "8.4", + pages = "250--252", +} diff --git a/latex/sections/methods.tex b/latex/sections/methods.tex index 409b5bc..42348f6 100644 --- a/latex/sections/methods.tex +++ b/latex/sections/methods.tex @@ -49,11 +49,6 @@ Since \eqref{eq:motion_x} and \eqref{eq:motion_y} are coupled, we want to rewrit \end{align*} - - - - - Physical properties given by newtons second law \eqref{eq:newton_second} \begin{equation}\label{eq:general_solution} f(t) = A_{+}e^{-i(\omega_{+} t + \phi_{+})} + A_{-}e^{-i(\omega_{-} t + \phi_{-})} @@ -66,4 +61,230 @@ The particle moves and its position can be determined using newton. where the el \subsection*{Tools} We used matplotlib -\end{document} \ No newline at end of file +\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}