aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorRahiel Kasim <rahielkasim@gmail.com>2016-06-08 22:57:58 +0200
committerRahiel Kasim <rahielkasim@gmail.com>2016-06-08 22:57:58 +0200
commitb671ceb869a324b2572ba16b08fa4d0358904f33 (patch)
tree3c0edfd5f349443afee438b592852217ec820412
parent67505fd319a4d69853844ea6d271df8ff90d5451 (diff)
import reports from sharelatex
-rw-r--r--.gitignore1
-rw-r--r--Set 1/Report/Pictures/time_ind.pngbin0 -> 260259 bytes
-rw-r--r--Set 1/Report/Pictures/time_ind_conv.pngbin0 -> 363438 bytes
-rw-r--r--Set 1/Report/Pictures/wave_eq.pngbin0 -> 230119 bytes
-rw-r--r--Set 1/Report/diffusion1D.pngbin0 -> 68647 bytes
-rw-r--r--Set 1/Report/diffusion_1.pngbin0 -> 19198 bytes
-rw-r--r--Set 1/Report/diffusion_2.pngbin0 -> 19766 bytes
-rw-r--r--Set 1/Report/diffusion_3.pngbin0 -> 20524 bytes
-rw-r--r--Set 1/Report/diffusion_4.pngbin0 -> 22251 bytes
-rw-r--r--Set 1/Report/diffusion_5.pngbin0 -> 22312 bytes
-rw-r--r--Set 1/Report/main.tex228
-rw-r--r--Set 1/Report/references.bib63
-rw-r--r--Set 1/Scientific_Computing_1.pdfbin0 -> 1072470 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta00.pngbin0 -> 65325 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta05.pngbin0 -> 106214 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta10.pngbin0 -> 304530 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta15.pngbin0 -> 214343 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta20.pngbin0 -> 217986 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta200.pngbin0 -> 78975 bytes
-rw-r--r--Set 2/Report/Pictures/DLA_eta30.pngbin0 -> 187746 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_1000.pngbin0 -> 91158 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_1600.pngbin0 -> 143315 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_3000.pngbin0 -> 320669 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_350.pngbin0 -> 40890 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_5000.pngbin0 -> 465791 bytes
-rw-r--r--Set 2/Report/Pictures/gs2_50000.pngbin0 -> 435347 bytes
-rw-r--r--Set 2/Report/Pictures/gs_0.pngbin0 -> 695565 bytes
-rw-r--r--Set 2/Report/Pictures/gs_10000.pngbin0 -> 633056 bytes
-rw-r--r--Set 2/Report/Pictures/gs_1300.pngbin0 -> 89891 bytes
-rw-r--r--Set 2/Report/Pictures/gs_2500.pngbin0 -> 144293 bytes
-rw-r--r--Set 2/Report/Pictures/gs_4000.pngbin0 -> 251802 bytes
-rw-r--r--Set 2/Report/Pictures/gs_700.pngbin0 -> 61667 bytes
-rw-r--r--Set 2/Report/Pictures/mcDLA_01.pngbin0 -> 126189 bytes
-rw-r--r--Set 2/Report/Pictures/mcDLA_02.pngbin0 -> 69079 bytes
-rw-r--r--Set 2/Report/Pictures/mcDLA_05.pngbin0 -> 72822 bytes
-rw-r--r--Set 2/Report/Pictures/mcDLA_10.pngbin0 -> 120178 bytes
-rw-r--r--Set 2/Report/Pictures/mcDLA_10b.pngbin0 -> 139618 bytes
-rw-r--r--Set 2/Report/main.tex254
-rw-r--r--Set 2/Report/references.bib56
-rw-r--r--Set 2/Scientific_Computing_2.pdfbin0 -> 5351227 bytes
-rw-r--r--Set 3/Report/Pictures/L_dependence.pngbin0 -> 181660 bytes
-rw-r--r--Set 3/Report/Pictures/diffusion.pngbin0 -> 289799 bytes
-rw-r--r--Set 3/Report/Pictures/dstep_dependence.pngbin0 -> 130337 bytes
-rw-r--r--Set 3/Report/Pictures/eigenmodes_c.pngbin0 -> 1018289 bytes
-rw-r--r--Set 3/Report/Pictures/eigenmodes_r.pngbin0 -> 612337 bytes
-rw-r--r--Set 3/Report/Pictures/eigenmodes_s.pngbin0 -> 1067725 bytes
-rw-r--r--Set 3/Report/Pictures/energies.pngbin0 -> 161447 bytes
-rw-r--r--Set 3/Report/Pictures/finite_scaled_1.pngbin0 -> 47302 bytes
-rw-r--r--Set 3/Report/Pictures/finite_scaled_2.pngbin0 -> 53668 bytes
-rw-r--r--Set 3/Report/Pictures/finite_scaled_3.pngbin0 -> 57743 bytes
-rw-r--r--Set 3/Report/Pictures/finite_scaled_4.pngbin0 -> 60087 bytes
-rw-r--r--Set 3/Report/Pictures/infinite_1.pngbin0 -> 35184 bytes
-rw-r--r--Set 3/Report/Pictures/infinite_2.pngbin0 -> 42317 bytes
-rw-r--r--Set 3/Report/Pictures/infinite_3.pngbin0 -> 48916 bytes
-rw-r--r--Set 3/Report/Pictures/infinite_4.pngbin0 -> 64333 bytes
-rw-r--r--Set 3/Report/Pictures/parabolic_1.pngbin0 -> 57875 bytes
-rw-r--r--Set 3/Report/Pictures/parabolic_2.pngbin0 -> 60494 bytes
-rw-r--r--Set 3/Report/Pictures/parabolic_3.pngbin0 -> 63378 bytes
-rw-r--r--Set 3/Report/Pictures/parabolic_4.pngbin0 -> 66230 bytes
-rw-r--r--Set 3/Report/Pictures/stencil.pngbin0 -> 9492 bytes
-rw-r--r--Set 3/Report/Pictures/time_evolution_c.pngbin0 -> 941772 bytes
-rw-r--r--Set 3/Report/Pictures/time_evolution_r.pngbin0 -> 1060053 bytes
-rw-r--r--Set 3/Report/Pictures/time_evolution_s.pngbin0 -> 999724 bytes
-rw-r--r--Set 3/Report/main.tex441
-rw-r--r--Set 3/Report/references.bib70
-rw-r--r--Set 3/Scientific_Computing_3.pdfbin0 -> 9374317 bytes
66 files changed, 1112 insertions, 1 deletions
diff --git a/.gitignore b/.gitignore
deleted file mode 100644
index e33609d..0000000
--- a/.gitignore
+++ /dev/null
@@ -1 +0,0 @@
-*.png
diff --git a/Set 1/Report/Pictures/time_ind.png b/Set 1/Report/Pictures/time_ind.png
new file mode 100644
index 0000000..fa26423
--- /dev/null
+++ b/Set 1/Report/Pictures/time_ind.png
Binary files differ
diff --git a/Set 1/Report/Pictures/time_ind_conv.png b/Set 1/Report/Pictures/time_ind_conv.png
new file mode 100644
index 0000000..2685a8e
--- /dev/null
+++ b/Set 1/Report/Pictures/time_ind_conv.png
Binary files differ
diff --git a/Set 1/Report/Pictures/wave_eq.png b/Set 1/Report/Pictures/wave_eq.png
new file mode 100644
index 0000000..35b3d31
--- /dev/null
+++ b/Set 1/Report/Pictures/wave_eq.png
Binary files differ
diff --git a/Set 1/Report/diffusion1D.png b/Set 1/Report/diffusion1D.png
new file mode 100644
index 0000000..cacd34f
--- /dev/null
+++ b/Set 1/Report/diffusion1D.png
Binary files differ
diff --git a/Set 1/Report/diffusion_1.png b/Set 1/Report/diffusion_1.png
new file mode 100644
index 0000000..9f917d0
--- /dev/null
+++ b/Set 1/Report/diffusion_1.png
Binary files differ
diff --git a/Set 1/Report/diffusion_2.png b/Set 1/Report/diffusion_2.png
new file mode 100644
index 0000000..d7e5dc7
--- /dev/null
+++ b/Set 1/Report/diffusion_2.png
Binary files differ
diff --git a/Set 1/Report/diffusion_3.png b/Set 1/Report/diffusion_3.png
new file mode 100644
index 0000000..cc4eee7
--- /dev/null
+++ b/Set 1/Report/diffusion_3.png
Binary files differ
diff --git a/Set 1/Report/diffusion_4.png b/Set 1/Report/diffusion_4.png
new file mode 100644
index 0000000..97288dd
--- /dev/null
+++ b/Set 1/Report/diffusion_4.png
Binary files differ
diff --git a/Set 1/Report/diffusion_5.png b/Set 1/Report/diffusion_5.png
new file mode 100644
index 0000000..6171753
--- /dev/null
+++ b/Set 1/Report/diffusion_5.png
Binary files differ
diff --git a/Set 1/Report/main.tex b/Set 1/Report/main.tex
new file mode 100644
index 0000000..3a05cd5
--- /dev/null
+++ b/Set 1/Report/main.tex
@@ -0,0 +1,228 @@
+\documentclass[a4paper]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage{pdfpages, titling}
+\usepackage{array, float}
+\usepackage{cite, tablefootnote}
+\usepackage{graphicx, caption, subfigure,subfig, wrapfig}
+\usepackage{listings}
+\usepackage{color}
+\usepackage{mathtools, braket}
+\usepackage{amssymb}
+\usepackage[nottoc]{tocbibind} % references in the toc
+% \usepackage{subcaption}
+\usepackage{textcomp}
+\usepackage{minted}
+\usepackage[hidelinks]{hyperref}
+\graphicspath{{Pictures/}} % Specifies the directory where pictures are stored
+
+\newcommand{\vect}[1]{\boldsymbol{#1}}
+\newcommand{\subtitle}[1]{
+ \posttitle{
+ \par\end{center}
+ \begin{center}\large#1\end{center}
+ \vskip0.5em}
+}
+\lstset{frame=tb,
+ language=Java,
+ aboveskip=3mm,
+ belowskip=3mm,
+ showstringspaces=false,
+ columns=flexible,
+ basicstyle={\small\ttfamily},
+ numbers=none,
+ breaklines=true,
+ breakatwhitespace=true,
+ tabsize=3
+}
+
+\author{Jaro Camphuijsen (6042473) and Rahiel Kasim (10447539)}
+\title{Scientific Computing I}
+\subtitle{Exploring the SciCoPyth's mind}
+\begin{document}
+\maketitle
+
+\tableofcontents
+
+\section{Tools}
+\textbf{Rahiel}: I've been using Linux for many years now. My favourite programming tool (and favourite software all-round) is GNU Emacs. It is a very advanced text editor. You can better understand Emacs if you think of it as a Lisp virtual machine with many functions for text editing. This means that any editing operation is actually a lisp function, so you can use lisp itself to extend the editor in any way you want. You however don't have to do everything yourself because Emacs has a built-in package manager, with more than 3000 packages available. For newcomers I recommend using \url{www.spacemacs.org}, a community maintained Emacs configuration with a focus on Vim emulation.
+
+Another tool that is more easily adopted in someone's existing workflow is z (\url{www.github.com/rupa/z/}). It tracks which directories you use most often so you can more easily jump to them. For example, if I type \texttt{z bio} it cd's my shell to \texttt{\~{}/Dropbox/univ/Computational Biology}, and \texttt{z chrom} brings me to \texttt{\~{}/Code/archiveror/chromium}.
+
+I am going to become a better programmer simply by programming more and by studying more ways how people program to solve problems.
+\\ \\
+\noindent
+\textbf{Jaro}: I have used Linux in the past on some other machines however after two months in my Masters Computational Science of using Windows on my own laptop, Rahiel convinced me to switch to Linux. It has made my life easier from the very beginning, and with the text editor Emacs almost anything is possible. It has a steeper learning curve then other text editors, but once you know the basic key combination the more complex operations become very easy and less time consuming. \\
+For sharing our code Rahiel and I like to use GitHub~\cite{sunsistemoGH}, it keeps track of changes and you can find back older versions at any time. Another nice function of GitHub is adding so called "issues" to your repository, this you can keep track of bugs and it also gives you a to do list towards the end of your project. \\
+For reporting we like to use \LaTeX \ because of its easy and clean layout generation, it automatically places figures and tables at good locations and keeps track of the reference numbers. \LaTeX \ files can also be shared on GitHub, however better tools exist like \url{www.sharelatex.com} and \url{www.overleaf.com}, where you can edit your online \LaTeX \ files interactively.
+
+\section{Vibrating String}
+The time evolution of a vibrating string is given by the partial differential equation \ref{eq:wave1D}
+
+\begin{equation}
+\frac{\partial ^2 \Psi}{\partial t^2} = c^2 \frac{\partial ^2 \Psi}{\partial dx^2}
+\label{eq:wave1D}
+\end{equation}
+
+To get a unique solution for $\Psi (x,t)$ we have to declare boundary and initial conditions.
+
+% \begin{subequations}
+% \begin{align}
+% \Psi (x, t_1) = f_1(x)\\
+% \Psi '(x,t_2) = f_2(x)\\
+% \Psi (x_1, t) = g_1(t)\\
+% \Psi (x_2,t) = g_2(t)
+% \label{eq:wave1Dinit}
+% \end{align}
+% \end{subequations}
+
+\noindent\begin{minipage}{.5\linewidth}
+\begin{equation*}
+\Psi (x, 0) = f(x)
+\end{equation*}
+\end{minipage}
+\begin{minipage}{.5\linewidth}
+\begin{equation*}
+\Psi (0, t) = 0
+\end{equation*}
+\end{minipage}\\
+\begin{minipage}{.5\linewidth}
+\begin{equation*}
+\Psi '(x,0) = 0
+\end{equation*}
+\end{minipage}
+\begin{minipage}{.5\linewidth}
+\begin{equation*}
+\Psi (L,t) = 0
+\end{equation*}
+\end{minipage}
+\label{eq:wave1Dinit}
+
+For the initial position of the string we use three different cases:
+\begin{subequations}
+\begin{align}
+f(x) &= \sin (2 \pi x)\\
+\label{eq:init2} f(x) &= \sin (5 \pi x)\\
+f(x) &= \sin (5 \pi x) \qquad \text{if } \frac{1}{5} \le x \le \frac{2}{5} \text{, else } \Psi = 0
+\end{align}
+\end{subequations}
+
+Only in a few cases can this PDE be solved analytically, hence we would like to construct a numerical model. To do this we discretize the continuous problem. For time step $k$ and spatial step $i$ we can write for the partial second derivatives of $\Psi$:
+
+\begin{subequations}
+\begin{align}
+\frac{\partial ^2 \Psi}{\partial x^2} &= \frac{u^{k}_{i+1} - 2 u^{k}_{i} + u^{k}_{i-1}}{\Delta x^2}\\
+\frac{\partial ^2 \Psi}{\partial t^2} &= \frac{u^{k+1}_{i} - 2 u^{k}_{i} + u^{k-1}_{i}}{\Delta t^2}
+\end{align}
+\end{subequations}
+
+Now we can solve the discrete (1D) wave equation for $u^{k+1}_{i}$ so we can compute the new value of each spatial bin $i$ from the current and previous time steps and its surrounding bins \cite{sciComp}:
+
+\begin{equation}
+u^{k+1}_{i} = 2 u^{k}_{i} - u^{k-1}_{i} + c \left ( \frac{\Delta t}{\Delta x} \right ) ^2 \left ( u^{k}_{i+1} - 2 u^{k}_{i} + u^{k}_{i-1} \right ), \qquad i = 1, ... , n
+\end{equation}
+
+We've implemented this scheme and in figure \ref{fig:wave} we see solutions to the wave equation with the initial conditions of eq. \ref{eq:init2}. We also have animations that better elucidate the time evolution and the inevitable numerical instability.
+
+\begin{figure}
+\center{\includegraphics[height=8cm]{wave_eq}}
+\caption{Solutions to the wave equation calculated with the finite difference method described in the text.}
+\label{fig:wave}
+\end{figure}
+
+% Same in code:
+% \begin{minted}{python}
+% for i in range(1, N-1)
+% u2[i] = 2u1[i] - u0[i] + c (dt/dx)^2 (u1[i+1] - 2 u1[i] + u1[i-1])
+% \end{minted}
+% \label{code:wave1D}
+
+
+
+\section{The Time Dependent Diffusion Equation}
+
+\begin{equation}
+c^{k+1}_{i,j} = c^{k}_{i,j} + \frac{\Delta t}{\Delta x^2} \left ( u^{k}_{i+1,j} +u^{k}_{i-1,j} + u^{k}_{i,j+1} + u^{k}_{i,j+1}
+ \right ), \qquad i = 1, ... , n
+\end{equation}
+
+At the boundaries of the x domain from we use the following function values to make periodic boundary conditions:
+
+\begin{subequations}
+\begin{align}
+u^{k}_{N+1,j} = u^{k}_{0,j}\\
+u^{k}_{-1,j} = u^{k}_{N,j}
+\end{align}
+\end{subequations}
+
+for the boundaries of the y domain we lock the values of $u^{k}_{i,0}$ and $u^{k}_{i,N}$ to respectively 1 and 0. We then let the simulation only generate values in the y domain \{1, N-1 \}.
+
+The simulation was done with N=30 spatial bins in both the x and y direction, the value of the time step was determined by the stability constraint:
+
+\begin{equation}
+ \Delta t = 0.9 \frac{\Delta x^2}{4 D}
+\end{equation}
+
+For this specific configuration we can compute the analytic solution:
+
+\begin{equation}
+ c(x,t) = 1 - \sum\limits_{i=1}^{\inf} \left [ \text{erfc} \left (\frac{1 - x + 2i}{2 \sqrt{D t}} \right ) - \text{erfc} \left (\frac{1 + x + 2i}{2 \sqrt{D t}} \right ) \right]
+\end{equation}
+
+And so we can compare the analytic solution with our model. The result for at different points in time can be seen in figure \ref{fig:1Ddiffusion} where we dispose of the x dimension. This can be done without loss of information since $u$ is constant along the x dimension, due to the configuration of this simulation. The full two dimensional outcome can be found in figure \ref{fig:2Ddiffusion} where heatmaps are drawn for the same points in time as in figure \ref{fig:1Ddiffusion}. In addition an animated plot was made which can be viewed by running the code.
+
+\begin{figure}
+\center{\includegraphics[height=8cm]{diffusion1D.png}}
+\caption{Comparison of the analytic solution (dashed line) and the numerical model at several points in time for the time dependent diffusion equation }
+\label{fig:1Ddiffusion}
+\end{figure}
+
+\begin{figure}[h]
+\centering
+\subfloat[]{
+\includegraphics[width=0.4\textwidth]{diffusion_1.png}
+}
+\subfloat[]{
+\includegraphics[width=0.4\textwidth]{diffusion_2.png}
+}\\
+\subfloat[]{
+\includegraphics[width=0.4\textwidth]{diffusion_3.png}
+}
+\subfloat[]{
+\includegraphics[width=0.4\textwidth]{diffusion_4.png}
+}\\
+\subfloat[]{
+\includegraphics[width=0.4\textwidth]{diffusion_5.png}
+}
+\caption{Several stages in the time dependent diffusion simulation in a 2 dimensional heatmap, the plots correspond to respectively t = [0, 0.001, 0.01, 0.1, 1]}
+\label{fig:2Ddiffusion}
+\end{figure}
+
+
+\section{The Time Independent Diffusion Equation}
+We've implemented the Jacobi iteration, Gauss-Seidel method and SOR with $N=50$ and a tolerance of $\epsilon = 10^{-5}$. In figure \ref{fig:time_ind} we compare the results of the three methods with the analytic solution. We see that they all come very close to the analytic solution: the result of SOR is superimposed on the analytic solution ($\lim_{t \to \infty} c(y, t)=-y$), while the Gauss-Seidel is further down and the Jacobi iteration is still below that. The nice thing is that this is also the same order for efficiency of the methods, SOR is both the most efficient (least amount of iterations) and the most accurate (closest to the analytic solution).
+
+\begin{figure}
+\center{\includegraphics[height=8cm]{time_ind}}
+\caption{A cross section of the solution of the time independent diffusion equation for different methods..}
+\label{fig:time_ind}
+\end{figure}
+
+We can further investigate the efficiency of the methods we've implemented by looking at the convergence measure $\delta$ versus the number of iterations $k$. In figure \ref{fig:conv} we see this relation for the three methods, we again see that the SOR method is the most efficient as it converges at the lowest number of iterations.
+
+\begin{figure}
+\center{\includegraphics[height=8cm]{time_ind_conv}}
+\caption{The convergence measure vs. the number of iterations for the methods solving the time independent diffusion equation. We see that the SOR method converges with the least amount of iterations.}
+\label{fig:conv}
+\end{figure}
+
+We can find the optimal $\omega$ for the SOR method by doing the SOR computation with different values for $\omega$ between the specified range while looking for the $\omega$ that needs the least amount of iterations. One can find $\omega$ naively by trying all values, or by doing a (binary) search over the range. We found as optimal $\omega = 1.891$. We can also look at the optimal $\omega$ for different values of $N$, using the same method we found that smaller values of $N$ have a smaller optimal $\omega$ (for $N=40$ we found $\omega=1.87$) and for larger $N$ we also find larger $\omega$ (for $N=60$ we found $\omega=1.91$).
+
+We put a 10 x 10 square in the middle of the domain and find that the convergence is faster, it took 204 iterations and without an object it takes 250 iterations. We find a slower convergence of 224 iterations if we put two 5 x 10 rectangles in the domain. Using the same method to find the optimal $\omega$ we find with the 10 x 10 square the optimum $\omega=1.87$ with 159 iterations, so adding an object slightly lowered the optimal $\omega$. For the two 5 x 10 rectangles we find an optimum of $\omega=1.89$ with 222 iterations, so here the object very slightly lowered the optimal.
+
+For these two examples we found that adding objects lowered the iterations and may also lower the optimal $\omega$ value.
+
+\bibliography{references}
+\bibliographystyle{unsrt}
+
+\end{document}
diff --git a/Set 1/Report/references.bib b/Set 1/Report/references.bib
new file mode 100644
index 0000000..fc80fe0
--- /dev/null
+++ b/Set 1/Report/references.bib
@@ -0,0 +1,63 @@
+@misc{wiki:lj,
+ author = "Wikipedia",
+ title = "{Lennard-Jones potential --- Wikipedia{,} The Free Encyclopedia}",
+ year = "2015",
+ note = "[{\url{https://en.wikipedia.org/w/index.php?title=Lennard-Jones_potential&oldid=694525432}}; accessed 27-January-2016]"
+}
+
+
+@misc{overleaf,
+ author = "{Overleaf}",
+ title = "{Overleaf: Collaborative Writing and Publishing}",
+ note = "[{\url{https://www.overleaf.com/}}; accessed 12-April-2016]"
+}
+
+@article{multi,
+author = {Glenn J. Martyna and Mark E. Tuckerman and Douglas J. Tobias and Michael L. Klein},
+title = {Explicit reversible integrators for extended systems dynamics},
+journal = {Molecular Physics},
+volume = {87},
+number = {5},
+pages = {1117-1157},
+year = {1996},
+doi = {10.1080/00268979600100761},
+URL = {http://dx.doi.org/10.1080/00268979600100761},
+eprint = {http://dx.doi.org/10.1080/00268979600100761}
+}
+
+@book{numMeth,
+ author = "William F. Ames",
+ title = "{Numerical Methods for Partial Differential Equations}",
+ publisher = "Thomas Nelson and Sons LTD",
+ year = "1969"
+}
+
+@book{sciComp,
+ author = "Micheal T. Heath",
+ title = "{Scientific Computing: An Introductory Survey}",
+ publisher = "Elizabeth A. Jones",
+ edition = "2nd",
+ year = "2002",
+ ISBN = "0-07-112229-x"
+}
+
+
+@misc{sunsistemoGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/sunsistemo}}; accessed 12-April-2016]"
+}
+
+@misc{sunsistemo,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo - The N-Body Simulator}",
+ note = "[{\url{https://sunsistemo.js.org}}; accessed 12-April-2016]"
+}
+
+@misc{SciCoPythGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{SciCoPyth on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/SciCoPyth}}; accessed 12-April-2016]"
+}
+
+
diff --git a/Set 1/Scientific_Computing_1.pdf b/Set 1/Scientific_Computing_1.pdf
new file mode 100644
index 0000000..aba617e
--- /dev/null
+++ b/Set 1/Scientific_Computing_1.pdf
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta00.png b/Set 2/Report/Pictures/DLA_eta00.png
new file mode 100644
index 0000000..0479da4
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta00.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta05.png b/Set 2/Report/Pictures/DLA_eta05.png
new file mode 100644
index 0000000..334d471
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta05.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta10.png b/Set 2/Report/Pictures/DLA_eta10.png
new file mode 100644
index 0000000..6832e13
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta10.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta15.png b/Set 2/Report/Pictures/DLA_eta15.png
new file mode 100644
index 0000000..7063118
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta15.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta20.png b/Set 2/Report/Pictures/DLA_eta20.png
new file mode 100644
index 0000000..5f7543e
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta20.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta200.png b/Set 2/Report/Pictures/DLA_eta200.png
new file mode 100644
index 0000000..da15f54
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta200.png
Binary files differ
diff --git a/Set 2/Report/Pictures/DLA_eta30.png b/Set 2/Report/Pictures/DLA_eta30.png
new file mode 100644
index 0000000..25d727c
--- /dev/null
+++ b/Set 2/Report/Pictures/DLA_eta30.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_1000.png b/Set 2/Report/Pictures/gs2_1000.png
new file mode 100644
index 0000000..c7b47b2
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_1000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_1600.png b/Set 2/Report/Pictures/gs2_1600.png
new file mode 100644
index 0000000..dab648f
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_1600.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_3000.png b/Set 2/Report/Pictures/gs2_3000.png
new file mode 100644
index 0000000..b67fbc6
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_3000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_350.png b/Set 2/Report/Pictures/gs2_350.png
new file mode 100644
index 0000000..a268d50
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_350.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_5000.png b/Set 2/Report/Pictures/gs2_5000.png
new file mode 100644
index 0000000..68f5789
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_5000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs2_50000.png b/Set 2/Report/Pictures/gs2_50000.png
new file mode 100644
index 0000000..075151b
--- /dev/null
+++ b/Set 2/Report/Pictures/gs2_50000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_0.png b/Set 2/Report/Pictures/gs_0.png
new file mode 100644
index 0000000..07515c8
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_0.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_10000.png b/Set 2/Report/Pictures/gs_10000.png
new file mode 100644
index 0000000..1b8c47e
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_10000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_1300.png b/Set 2/Report/Pictures/gs_1300.png
new file mode 100644
index 0000000..655fdc3
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_1300.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_2500.png b/Set 2/Report/Pictures/gs_2500.png
new file mode 100644
index 0000000..acba7b0
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_2500.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_4000.png b/Set 2/Report/Pictures/gs_4000.png
new file mode 100644
index 0000000..4467f8f
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_4000.png
Binary files differ
diff --git a/Set 2/Report/Pictures/gs_700.png b/Set 2/Report/Pictures/gs_700.png
new file mode 100644
index 0000000..8aa85db
--- /dev/null
+++ b/Set 2/Report/Pictures/gs_700.png
Binary files differ
diff --git a/Set 2/Report/Pictures/mcDLA_01.png b/Set 2/Report/Pictures/mcDLA_01.png
new file mode 100644
index 0000000..a0b1c47
--- /dev/null
+++ b/Set 2/Report/Pictures/mcDLA_01.png
Binary files differ
diff --git a/Set 2/Report/Pictures/mcDLA_02.png b/Set 2/Report/Pictures/mcDLA_02.png
new file mode 100644
index 0000000..b2e36fc
--- /dev/null
+++ b/Set 2/Report/Pictures/mcDLA_02.png
Binary files differ
diff --git a/Set 2/Report/Pictures/mcDLA_05.png b/Set 2/Report/Pictures/mcDLA_05.png
new file mode 100644
index 0000000..2c7f310
--- /dev/null
+++ b/Set 2/Report/Pictures/mcDLA_05.png
Binary files differ
diff --git a/Set 2/Report/Pictures/mcDLA_10.png b/Set 2/Report/Pictures/mcDLA_10.png
new file mode 100644
index 0000000..2282035
--- /dev/null
+++ b/Set 2/Report/Pictures/mcDLA_10.png
Binary files differ
diff --git a/Set 2/Report/Pictures/mcDLA_10b.png b/Set 2/Report/Pictures/mcDLA_10b.png
new file mode 100644
index 0000000..2e1eea1
--- /dev/null
+++ b/Set 2/Report/Pictures/mcDLA_10b.png
Binary files differ
diff --git a/Set 2/Report/main.tex b/Set 2/Report/main.tex
new file mode 100644
index 0000000..d0bc766
--- /dev/null
+++ b/Set 2/Report/main.tex
@@ -0,0 +1,254 @@
+\documentclass[a4paper]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage{pdfpages, titling}
+\usepackage{array, float}
+\usepackage{cite, tablefootnote}
+\usepackage{graphicx, caption}
+\usepackage{listings}
+\usepackage{color}
+\usepackage{mathtools, braket}
+\usepackage{amssymb}
+\usepackage[nottoc]{tocbibind} % references in the toc
+\usepackage{subcaption}
+\usepackage{textcomp}
+\usepackage{minted}
+\usepackage[hidelinks]{hyperref}
+\graphicspath{{Pictures/}} % Specifies the directory where pictures are stored
+
+\newcommand{\vect}[1]{\boldsymbol{#1}}
+\newcommand{\subtitle}[1]{
+ \posttitle{
+ \par\end{center}
+ \begin{center}\large#1\end{center}
+ \vskip0.5em}
+}
+\lstset{frame=tb,
+ language=Java,
+ aboveskip=3mm,
+ belowskip=3mm,
+ showstringspaces=false,
+ columns=flexible,
+ basicstyle={\small\ttfamily},
+ numbers=none,
+ breaklines=true,
+ breakatwhitespace=true,
+ tabsize=3
+}
+
+\author{Jaro Camphuijsen (6042473) and Rahiel Kasim (10447539)}
+\title{Scientific Computing II}
+\subtitle{The Return of the SciCoPyth}
+
+\begin{document}
+\maketitle
+
+\tableofcontents
+
+\section{Diffusion Limited Aggregation}
+Diffusion Limited Aggregation (DLA) considers diffusing particles that that can stick to a certain structure and thereby expand the structure. Because we are dealing with diffusion we should solve the diffusion equation at each growth step to calculate the new concentration of particles.
+
+\begin{equation}
+\frac{\partial \phi (r,t)}{\partial t} = D \nabla ^2 \phi (r,t)
+\label{eq:diff}
+\end{equation}
+
+Because the growth of such a structure (e.g. a coral or mineral structure) is on a very long timescale and the structure growth per step is very small we can assume the diffusion reaches equilibrium before the next step and we can solve the time independent diffusion equation.
+
+The update scheme of a growth simulation of a certain number of \mintinline{python}{steps} is given in the following pseudo code snippet:
+
+\begin{minted}[]{python}
+def grow(eta, omega, steps):
+ for n in range(steps):
+ nutLattice = SOR(omega, objLattice)
+ objLattice = grow_step(eta, nutLattice)
+\end{minted}
+
+Here we have the nutrient lattice \mintinline{python}{nutLattice} which is updated each time step by the \mintinline{python}{SOR} function and contains the nutrient concentration at each cell, it uses the object lattice \mintinline{python}{objLattice} for its boundary conditions. The object lattice is updated using the \mintinline{python}{grow_step} function to calculate the probability for each cell to become part of the structure from the nutrient concentration in the cell.
+
+To simulate the growth we define two different cell states: a cell can be part of the structure or it can be empty. Each step we determine all grow candidates. For a cell to be a grow candidate it must be empty and has one of his neighbours be part of the structure. Each candidate has the following probability to switch it state from empty to be part of the structure:
+
+\begin{equation}\label{eq:grow_prob}
+ p_g (i,j) = \frac{c_{i,j}^{\eta}}{\sum_{\text{candidates}} c_{i,j}^{\eta}}
+\end{equation}
+
+For the calculation of the time independent diffusion equation we use the successive over-relaxation (SOR) function from the last assignment. In the previous assignment \cite{scico1} we also optimized the relaxation parameter $\omega$, however the optimal value of $\omega$ depends on the boundary values of the domain. These were fixed in the simple example of the previous exercise but since we have a growing object this time, the value of $\omega$ should in principle be optimized every time the structure has grown a bit. This to computationally demanding to be of any use.
+
+When we do a simulation with a growing parameter $0 < \eta < 1$ we get a problem using the SOR calculation. Due to the over relaxation it is possible to end up with a small negative concentration. This is no problem for other values of $\eta$, however since the root of a negative number is not defined in our program we end up with an error. To avoid this problem in calculation and because we can not truly optimize $\omega$ which would result in a minor speed improvement anyway, we should set $\omega = 1$ or use under relaxation for all simulations with $0 < \eta < 1$.
+
+\begin{figure}[ht]
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={7cm 1.8cm 7cm 1cm},clip]{DLA_eta00.png}
+ \caption{$\eta = 0$}
+ \label{fig:DLA_eta00}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={0cm 1cm 0cm 0cm},clip]{DLA_eta05.png}
+ \caption{$\eta = 0.5$}
+ \label{fig:DLA_eta05}
+ \end{subfigure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={7cm 1.8cm 7cm 1cm},clip]{DLA_eta10.png}
+ \caption{$\eta = 1$}
+ \label{fig:DLA_eta10}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={7cm 1.8cm 7cm 1cm},clip]{DLA_eta15.png}
+ \caption{$\eta = 1.5$}
+ \label{fig:DLA_eta15}
+ \end{subfigure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={7cm 1.8cm 7cm 1cm},clip]{DLA_eta20.png}
+ \caption{$\eta = 2$}
+ \label{fig:DLA_eta20}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={7cm 1.8cm 7cm 1cm},clip]{DLA_eta30.png}
+ \caption{$\eta = 3$}
+ \label{fig:DLA_eta30}
+ \end{subfigure}
+ \caption{Diffusion Limited Aggregation for different values of $\eta$}
+ \label{fig:DLA_eta}
+\end{figure}
+
+For varying growing parameter $\eta$ we get different structures as shown in figure \ref{fig:DLA_eta}. If we set $\eta$ to zero as in figure \ref{fig:DLA_eta00} we turn off the influence of the nutrient concentration. Each grow candidate has the same probability of growing which results in a circular bulk around the initial point. With increasing $\eta$ the bulk tends to grow towards more expanded structures with normal DLA growth for $\eta = 1$ in figure \ref{fig:DLA_eta10}. If we keep increasing $\eta$ the influence of the nutrient concentration increases and since we have a gradient directed to the top of the domain, the structure will grow more linearly with less pronounced branches. In the case of $\eta = 3$ as seen in figure \ref{fig:DLA_eta30}, the structure is almost completely linear this is because candidate cells at the top of the structure have a much higher probability to become part of the object due to the concentration gradient.
+
+
+
+\section{Monte Carlo DLA}
+A different method to simulate DLA uses random walkers and therefore is a Monte Carlo type simulation. This type of simulation can in principle be seen as a true imitation of diffusion since diffusing particles in the real world also follow a random walk. For this simulation we created a class \mintinline{python}{Walker} with methods to take a random step, check its neighbours and merge with the structure in the object lattice whenever it touches the object. Additionally a sticking parameter ($p_{stick}$) can be implemented which determines the probability that the walker merges if it touches the structure. If this probability is one we are left with the original Monte Carlo simulation. The result for different values of this sticking probability can be seen in figure \ref{fig:mcDLA}. We can clearly see the similarity between Monte-Carlo DLA with $p_{stick} = 1$ in figure \ref{fig:mcDLA_10} and the normal DLA simulation in figure \ref{fig:DLA_eta10}. Also we see that for lower values of $p_{stick}$ we get more linear structures with less and smaller branches. This is equivalent to raising the growing parameter $\eta$ in the normal DLA simulation.
+
+The advantage of the Monte Carlo method is that we do not have to solve the diffusion equations which is the bottle neck in the DLA simulation, however doing a pure random until the particle sticks takes a long time. We can improve on the runtime dramatically by adding a small bias for steps downward as is done in some simulations. \cite{davies} However this breaks the similarity with true diffusion since we introduce a net flow in the system. Of course we can use this flow as it is present in several natural systems (e.g. gravity pulling on the nutrients in the formation of corals).
+
+\begin{figure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={0cm 0cm 0cm 0cm},clip]{mcDLA_10.png}
+ \caption{$p_{stick} = 1$}
+ \label{fig:mcDLA_10}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={1cm 0.9cm 3.07cm 0cm},clip]{mcDLA_05.png}
+ \caption{$p_{stick} = 0.5$}
+ \label{fig:mcDLA_05}
+ \end{subfigure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=\linewidth, trim={1cm 0.9cm 3.07cm 0cm},clip]{mcDLA_02.png}
+ \caption{$p_{stick} = 0.2$}
+ \label{fig:mcDLA_02}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=5.6cm]{mcDLA_01.png}
+ \caption{$p_{stick} = 0.1$}
+ \label{fig:mcDLA_01}
+ \end{subfigure}
+ \caption{A Monte Carlo simulation of Diffusion Limited Aggregation for different values of the sticking probability $p_{stick}$. These structures can be compared to the structures in figure \ref{fig:DLA_eta10}-\ref{fig:DLA_eta20}}
+ \label{fig:mcDLA}
+\end{figure}
+
+\section{The Gray-Scott Model}
+The Gray-Scott Model is a reaction-diffusion system. It models the flow of two chemicals, $u$ and $v$, by showing their concentrations over time as they diffuse and react with each other. Reaction-diffusion systems look like the familiar diffusion equation but they have an extra term that is a function of the relevant concentrations (e.g. $f(u, v)$). The Gray-Scott model can be expressed as:
+\begin{equation} \label{eq:gray_scott}
+\begin{aligned}
+ \frac{\partial u}{\partial t} &= D_u \nabla^2 u - uv^2 + f(1 - u), \\
+ \frac{\partial v}{\partial t} &= D_v \nabla^2 v + vv^2 - (f + k) v,
+\end{aligned}
+\end{equation}
+where $D_u$ and $D_v$ are the diffusion constants of the chemicals and $f$ and $k$ parameters of the underlying chemical reaction. To solve these partial differential equations for $u$ and $v$ we discretize the equations using finite difference methods. We divide 2D space in $N$ chunks of $\delta x$, $\delta y$ and time in periods of $\delta t$. In particular, we substitute
+
+\begin{align*}
+ \frac{\partial u}{\partial t} &= \frac{u^{k+1}_i - u^k_i}{\delta t}, \\
+ \begin{split}
+ \nabla^2 u &= \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)^2 u, \\
+ &= \frac{1}{\delta x^2} (u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{i,j}),
+ \end{split}
+\end{align*}
+in system \ref{eq:gray_scott} to get the discretization:
+\begin{equation*} \label{eq:gray_scott_discrete}
+\begin{aligned}
+ u^{k+1}_{i,j} &= u^k_{i,j} + \delta t \left(\frac{D_u}{\delta x^2} (u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4 u_{i,j}) - uv^2 + f(1-u) \right) \\
+ v^{k+1}_{i,j} &= v^k_{i,j} + \delta t \left(\frac{D_v}{\delta x^2} (v_{i+1,j} + v_{i-1,j} + v_{i,j+1} + v_{i,j-1} - 4 v_{i,j}) + uv^2 - (f + k) v \right)
+\end{aligned}
+\end{equation*}
+
+We take periodic boundary conditions along the $x$ and $y$ directions: $u(0, y, t) = u(N \delta x, y, t)$ and $u(x, 0, t) = u(x, N \delta y, t)$, with the same for $v$. We've implemented the derived discretization with these boundary conditions.
+
+For the parameters we've chosen $\delta t = \delta x = \delta y = 1$, $N=300$, $D_u = 0.16$, $D_v 0.08$, $f = 0.035$ and $k = 0.060$. Our initial conditions were $u = 0.5$ everywhere and $v = 0.25$ in a square around the origin and $v = 0$ elsewhere. We added a small amount of noise to the initial conditions of $u$ and $v$ of magnitude 0.05.
+
+In figure \ref{fig:gs_time1} we see the concentrations of $u$ and $v$ over time. First at $t=0$ we see the square of the initial conditions with the added noise seen as speckles all over the domain. In the middle at $t=700$ we see that the noise has disappeared due to diffusion. The square itself has also disappeared, leaving behind four circles growing at its former corners. And at the bottom at $t=1300$ the circles have grown into four star-like shapes.
+
+The process continues in figure \ref{fig:gs_time2}. The four stars have merged to a nice symmetric emblem. In the middle we see some breaking of symmetry, this could indicate the effect of numerical error, as the initial conditions were isotropic. Finally at the bottom the pattern has grown to extend over the boundaries. Because we implemented periodic boundary conditions, this final pattern can be used on tiles with the edges fitting perfectly.
+
+\begin{figure}
+\includegraphics[width=12cm]{gs_0}
+\includegraphics[width=12cm]{gs_700}
+\includegraphics[width=12cm]{gs_1300}
+\caption{The patterns that arise in the Gray-Scott model with the parameters described in the text. At the top we see the concentrations at time $t=0$, in the middle at $t=700$ and at the bottom at $t=1300$.}
+\label{fig:gs_time1}
+\end{figure}
+
+\begin{figure}
+ \centering
+ \includegraphics[width=12cm]{gs_2500}
+ \includegraphics[width=12cm]{gs_4000}
+ \includegraphics[width=12cm]{gs_10000}
+ \caption{The continuation of figure \ref{fig:gs_time1}. Top: $t=2500$, middle: $t=4000$, bottom: $t=10000$.}
+ \label{fig:gs_time2}
+\end{figure}
+
+We now explore the patterns formed when the parameters are changed slightly: $f=0.030$ and $k=0.055$. Everything else is the same. This time we omit $u$ from the plots, because it does not add new information on the patterns given $v$. In figure \ref{fig:gs2} we see the formation of the patterns. The initial square grows into a waffle (b and c) and continues expanding until we end up with many small mountains (f), or spots. We conclude by noting that the small change in parameters lead to quite different patterns.
+
+\begin{figure}[ht]
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_350}
+ \caption{$t = 350$}
+ \label{fig:gs2_350}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_1000}
+ \caption{$t = 1000$}
+ \label{fig:gs2_1000}
+ \end{subfigure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_1600}
+ \caption{$t = 1600$}
+ \label{fig:gs2_1600}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_3000}
+ \caption{$t = 3000$}
+ \label{fig:gs2_3000}
+ \end{subfigure}
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_5000}
+ \caption{$t = 5000$}
+ \label{fig:gs2_5000}
+ \end{subfigure}%%
+ \begin{subfigure}[b]{0.5\linewidth}
+ \centering
+ \includegraphics[width=7cm]{gs2_50000}
+ \caption{$t = 50000$}
+ \label{fig:gs2_50000}
+ \end{subfigure}
+ \caption{The Gray-Scott model with $f=0.030$ and $k=0.055$.}
+ \label{fig:gs2}
+\end{figure}
+
+\bibliography{references}
+\bibliographystyle{unsrt}
+
+\end{document}
diff --git a/Set 2/Report/references.bib b/Set 2/Report/references.bib
new file mode 100644
index 0000000..1302e74
--- /dev/null
+++ b/Set 2/Report/references.bib
@@ -0,0 +1,56 @@
+@misc{wiki:lj,
+ author = "Wikipedia",
+ title = "{Lennard-Jones potential --- Wikipedia{,} The Free Encyclopedia}",
+ year = "2015",
+ note = "[{\url{https://en.wikipedia.org/w/index.php?title=Lennard-Jones_potential&oldid=694525432}}; accessed 27-January-2016]"
+}
+
+
+@misc{davies,
+ author = "{Jason Davies}",
+ title = "{Diffusion-limited aggregation}",
+ note = "[{\url{https://www.jasondavies.com/dla/}}; accessed 10-May-2016]"
+}
+
+@unpublished{scico1,
+author = "{J. Camphuijsen, R. Kasim}",
+title = "{Scientific Computing 1 - Exploring the SciCoPyth's mind}",
+year = {2016}
+}
+
+@book{numMeth,
+ author = "William F. Ames",
+ title = "{Numerical Methods for Partial Differential Equations}",
+ publisher = "Thomas Nelson and Sons LTD",
+ year = "1969"
+}
+
+@book{sciComp,
+ author = "Micheal T. Heath",
+ title = "{Scientific Computing: An Introductory Survey}",
+ publisher = "Elizabeth A. Jones",
+ edition = "2nd",
+ year = "2002",
+ ISBN = "0-07-112229-x"
+}
+
+
+@misc{sunsistemoGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/sunsistemo}}; accessed 12-April-2016]"
+}
+
+@misc{sunsistemo,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo - The N-Body Simulator}",
+ note = "[{\url{https://sunsistemo.js.org}}; accessed 12-April-2016]"
+}
+
+@misc{SciCoPythGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{SciCoPyth on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/SciCoPyth}}; accessed 12-April-2016]"
+}
+
+
diff --git a/Set 2/Scientific_Computing_2.pdf b/Set 2/Scientific_Computing_2.pdf
new file mode 100644
index 0000000..f3f165b
--- /dev/null
+++ b/Set 2/Scientific_Computing_2.pdf
Binary files differ
diff --git a/Set 3/Report/Pictures/L_dependence.png b/Set 3/Report/Pictures/L_dependence.png
new file mode 100644
index 0000000..b858722
--- /dev/null
+++ b/Set 3/Report/Pictures/L_dependence.png
Binary files differ
diff --git a/Set 3/Report/Pictures/diffusion.png b/Set 3/Report/Pictures/diffusion.png
new file mode 100644
index 0000000..572f85a
--- /dev/null
+++ b/Set 3/Report/Pictures/diffusion.png
Binary files differ
diff --git a/Set 3/Report/Pictures/dstep_dependence.png b/Set 3/Report/Pictures/dstep_dependence.png
new file mode 100644
index 0000000..797d937
--- /dev/null
+++ b/Set 3/Report/Pictures/dstep_dependence.png
Binary files differ
diff --git a/Set 3/Report/Pictures/eigenmodes_c.png b/Set 3/Report/Pictures/eigenmodes_c.png
new file mode 100644
index 0000000..29436f2
--- /dev/null
+++ b/Set 3/Report/Pictures/eigenmodes_c.png
Binary files differ
diff --git a/Set 3/Report/Pictures/eigenmodes_r.png b/Set 3/Report/Pictures/eigenmodes_r.png
new file mode 100644
index 0000000..8a96070
--- /dev/null
+++ b/Set 3/Report/Pictures/eigenmodes_r.png
Binary files differ
diff --git a/Set 3/Report/Pictures/eigenmodes_s.png b/Set 3/Report/Pictures/eigenmodes_s.png
new file mode 100644
index 0000000..ce0dcd0
--- /dev/null
+++ b/Set 3/Report/Pictures/eigenmodes_s.png
Binary files differ
diff --git a/Set 3/Report/Pictures/energies.png b/Set 3/Report/Pictures/energies.png
new file mode 100644
index 0000000..fb6e3ce
--- /dev/null
+++ b/Set 3/Report/Pictures/energies.png
Binary files differ
diff --git a/Set 3/Report/Pictures/finite_scaled_1.png b/Set 3/Report/Pictures/finite_scaled_1.png
new file mode 100644
index 0000000..fc05658
--- /dev/null
+++ b/Set 3/Report/Pictures/finite_scaled_1.png
Binary files differ
diff --git a/Set 3/Report/Pictures/finite_scaled_2.png b/Set 3/Report/Pictures/finite_scaled_2.png
new file mode 100644
index 0000000..af81031
--- /dev/null
+++ b/Set 3/Report/Pictures/finite_scaled_2.png
Binary files differ
diff --git a/Set 3/Report/Pictures/finite_scaled_3.png b/Set 3/Report/Pictures/finite_scaled_3.png
new file mode 100644
index 0000000..8398dc3
--- /dev/null
+++ b/Set 3/Report/Pictures/finite_scaled_3.png
Binary files differ
diff --git a/Set 3/Report/Pictures/finite_scaled_4.png b/Set 3/Report/Pictures/finite_scaled_4.png
new file mode 100644
index 0000000..2e094ca
--- /dev/null
+++ b/Set 3/Report/Pictures/finite_scaled_4.png
Binary files differ
diff --git a/Set 3/Report/Pictures/infinite_1.png b/Set 3/Report/Pictures/infinite_1.png
new file mode 100644
index 0000000..7887f36
--- /dev/null
+++ b/Set 3/Report/Pictures/infinite_1.png
Binary files differ
diff --git a/Set 3/Report/Pictures/infinite_2.png b/Set 3/Report/Pictures/infinite_2.png
new file mode 100644
index 0000000..bb8f5f9
--- /dev/null
+++ b/Set 3/Report/Pictures/infinite_2.png
Binary files differ
diff --git a/Set 3/Report/Pictures/infinite_3.png b/Set 3/Report/Pictures/infinite_3.png
new file mode 100644
index 0000000..35f7a76
--- /dev/null
+++ b/Set 3/Report/Pictures/infinite_3.png
Binary files differ
diff --git a/Set 3/Report/Pictures/infinite_4.png b/Set 3/Report/Pictures/infinite_4.png
new file mode 100644
index 0000000..2b66527
--- /dev/null
+++ b/Set 3/Report/Pictures/infinite_4.png
Binary files differ
diff --git a/Set 3/Report/Pictures/parabolic_1.png b/Set 3/Report/Pictures/parabolic_1.png
new file mode 100644
index 0000000..a2dbfb3
--- /dev/null
+++ b/Set 3/Report/Pictures/parabolic_1.png
Binary files differ
diff --git a/Set 3/Report/Pictures/parabolic_2.png b/Set 3/Report/Pictures/parabolic_2.png
new file mode 100644
index 0000000..097af22
--- /dev/null
+++ b/Set 3/Report/Pictures/parabolic_2.png
Binary files differ
diff --git a/Set 3/Report/Pictures/parabolic_3.png b/Set 3/Report/Pictures/parabolic_3.png
new file mode 100644
index 0000000..bc7dd3f
--- /dev/null
+++ b/Set 3/Report/Pictures/parabolic_3.png
Binary files differ
diff --git a/Set 3/Report/Pictures/parabolic_4.png b/Set 3/Report/Pictures/parabolic_4.png
new file mode 100644
index 0000000..de27c7e
--- /dev/null
+++ b/Set 3/Report/Pictures/parabolic_4.png
Binary files differ
diff --git a/Set 3/Report/Pictures/stencil.png b/Set 3/Report/Pictures/stencil.png
new file mode 100644
index 0000000..6e6dbb4
--- /dev/null
+++ b/Set 3/Report/Pictures/stencil.png
Binary files differ
diff --git a/Set 3/Report/Pictures/time_evolution_c.png b/Set 3/Report/Pictures/time_evolution_c.png
new file mode 100644
index 0000000..7982ce8
--- /dev/null
+++ b/Set 3/Report/Pictures/time_evolution_c.png
Binary files differ
diff --git a/Set 3/Report/Pictures/time_evolution_r.png b/Set 3/Report/Pictures/time_evolution_r.png
new file mode 100644
index 0000000..80296b4
--- /dev/null
+++ b/Set 3/Report/Pictures/time_evolution_r.png
Binary files differ
diff --git a/Set 3/Report/Pictures/time_evolution_s.png b/Set 3/Report/Pictures/time_evolution_s.png
new file mode 100644
index 0000000..606db51
--- /dev/null
+++ b/Set 3/Report/Pictures/time_evolution_s.png
Binary files differ
diff --git a/Set 3/Report/main.tex b/Set 3/Report/main.tex
new file mode 100644
index 0000000..400f92f
--- /dev/null
+++ b/Set 3/Report/main.tex
@@ -0,0 +1,441 @@
+\documentclass[a4paper]{article}
+\usepackage[utf8]{inputenc}
+\usepackage[english]{babel}
+\usepackage{pdfpages, titling}
+\usepackage{array, float, arydshln}
+\usepackage{cite, tablefootnote}
+\usepackage{graphicx, caption}
+\usepackage{listings}
+\usepackage{color}
+\usepackage{mathtools, braket}
+\usepackage{amssymb}
+\usepackage[nottoc]{tocbibind} % references in the toc
+\usepackage{subcaption}
+\usepackage{textcomp}
+\usepackage{minted}
+\usepackage[hidelinks]{hyperref}
+\graphicspath{{Pictures/}} % Specifies the directory where pictures are stored
+
+\newcommand{\vect}[1]{\boldsymbol{#1}}
+\newcommand{\subtitle}[1]{
+ \posttitle{
+ \par\end{center}
+ \begin{center}\large#1\end{center}
+ \vskip0.5em}
+}
+\lstset{frame=tb,
+ language=Java,
+ aboveskip=3mm,
+ belowskip=3mm,
+ showstringspaces=false,
+ columns=flexible,
+ basicstyle={\small\ttfamily},
+ numbers=none,
+ breaklines=true,
+ breakatwhitespace=true,
+ tabsize=3
+}
+
+\author{Jaro Camphuijsen (6042473) and Rahiel Kasim (10447539)}
+\title{Scientific Computing III}
+\subtitle{SciCoPyth Waving Goodbye!}
+
+\begin{document}
+\maketitle
+
+\tableofcontents
+
+\section{Eigenmodes of drums or membranes of different shapes}
+\subsection{Theory}
+When we want to find the eigenmodes of a vibrating surface, we want to solve the two dimensional wave equation and find a seperable solution:
+\begin{equation}
+\frac{1}{T(t)}\frac{\partial ^2 T(t)}{\partial t^2} = c^2 \frac{1}{v(x,y)} \nabla ^2 v(x,y)
+\label{eq:wave2d}
+\end{equation}
+
+Because we are interested in the time independent eigenmodes we want to solve only the right hand side of equation \ref{eq:wave2d} where we can replace the left hand side by a constant K.\footnote{Both the left and right hand side of equation \ref{eq:wave2d} are constant since one depends on only t and one on only x and y. This is the separable behaviour of the solution we are looking for.} We end up with the following equation:
+
+\begin{equation}
+\nabla ^2 v(x,y) = K v(x,y)
+\label{eq:wavexy}
+\end{equation}
+
+Which we can discretize:
+
+\begin{equation}
+\frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{\Delta x^2} + \frac{v_{i,j+1} - 2v_{i,j} + v_{i,j-1}}{\Delta y^2} =
+K v_{i,j}
+\label{eq:waveDiscr}
+\end{equation}
+
+By setting $\Delta x = \Delta y = \sqrt{\lambda / K}$ and we can simplify this to:
+\begin{equation}
+ v_{i+1,j} + v_{i,j+1} + v_{i-1,j} + v_{i,j-1} - 4v_{i,j} = \lambda v_{i,j}
+\label{eq:waveDiscrSimp}
+\end{equation}
+
+We can write this out explicitly for a $2\times2$ system and construct the equations in matrix form. All boundaries are set to zero so we can drop all terms with $i$ or $j = 0$ or $3$ in the matrix equation:
+
+\begin{equation}
+\begin{aligned}
+v_{2,1} + v_{1,2} + v_{0,1} + v_{1,0} - 4v_{1,1} = \lambda v_{1,1} \\
+v_{3,1} + v_{2,2} + v_{1,1} + v_{2,0} - 4v_{2,1} = \lambda v_{2,1} \\
+v_{2,2} + v_{1,3} + v_{0,2} + v_{1,1} - 4v_{1,2} = \lambda v_{1,2}\\
+v_{3,2} + v_{2,3} + v_{1,2} + v_{2,1} - 4v_{2,2} = \lambda v_{2,2}
+\end{aligned}
+\label{eq:waveDiscrExpl}
+\end{equation}
+
+
+
+\begin{equation}
+\begin{bmatrix*}[r]
+ -4 & 1 & 1 & 0 \\
+ 1 & -4 & 0 & 1 \\
+ 1 & 0 & -4 & 1 \\
+ 0 & 1 & 1 & -4
+\end{bmatrix*}
+\begin{bmatrix*}
+ v_{1,1}\\
+ v_{2,1}\\
+ v_{1,2}\\
+ v_{2,2}
+\end{bmatrix*}
+=
+\lambda
+\begin{bmatrix*}
+ v_{1,1}\\
+ v_{2,1}\\
+ v_{1,2}\\
+ v_{2,2}
+\end{bmatrix*}
+\label{eq:mat2x2}
+\end{equation}
+Equation \ref{eq:mat2x2} is a so called eigenvalue problem, it states that for certain vectors $v_{i,j}$ multiplied by the finite difference matrix, we obtain a scalar multiple of this same vector $\lambda v_{i,j}$. These vectors and scalar factors are called eigenvectors and eigenvalues of the matrix respectively. An eigenvalue problem can be solved directly using the appropriate methods. The eigenvectors represent the possible standing waves in the domain and their eigenvalues give us the frequency of such a wave.
+
+\begin{figure}
+\begin{centering}
+\includegraphics[width=8cm]{stencil.png}
+\caption{The stencil for a $4\times4$ system with all connection for the point $(1,2)$.}
+\label{fig:stencil4x4}
+\end{centering}
+\end{figure}
+
+In figure \ref{fig:stencil4x4} we see the stencil for a $4\times4$ system. The connections of point $(1,2)$ are the points $(1,1)$, $(2,2)$, $(1,3)$ and the boundary point $(0,2)$. However since all boundaries are zero we can drop this last point. The matrix representation of the eigenvalue problem of a square system of arbitrary size can be found in equation \ref{eq:mat4x4} in appendix \ref{app:A}. The connections of point $(1,2)$ in figure \ref{fig:stencil4x4} correspond to the second column of this matrix.
+
+The solution to the timedependent part in equation \ref{eq:wave2d} where we replace the right hand side by the same constant K, is given by equation \ref{eq:timedepsol}.
+
+\begin{equation}
+T(t) = A\sin{c\lambda t} + B\cos{c\lambda t}
+\label{eq:timedepsol}
+\end{equation}
+
+
+\subsection{Implementation}
+For a square domain with sides n cells, the finite difference matrix grows with $n^4$. As we can see the bigger the system the more zero valued entries we have in our finite difference matrix. When implementing this problem it is not very efficient to store the full matrix but we can use some other format instead. The SciPy library \cite{scipy} offers some useful matrix formats in its \mintinline{python}{sparse} module which store only the nonzero values. This type of matrix saves memory space and iteration loops since we now only have approximately $5L$ entries.
+
+We also want to find the eigenvalues for different domain shapes. For a circular domain we set all rows and columns corresponding with points outside a certain radius from the domain center to zero. The same is done for a rectangular domain, this time using a domain length of $2L$ so the short side of the rectangle is of length $L$.
+
+Once found, the eigenmodes can be sorted by their eigenvalue, we also want to ignore all modes with eigenvalue equal to zero since these are not interesting to us because there is no wave in that case.
+
+Once the eigenmodes are found we can evolve them in time using the simple solution of the time dependent part of the wave equation given by equation \ref{eq:timedepsol}. Here we set A to zero so the sine part is discarded and we are only left with the cosine part, we do not lose generality here since we only set the previously found steady state to be amplitude of the wave.
+
+
+\subsection{Results}
+In figure \ref{fig:eigenmodes_s}, \ref{fig:eigenmodes_r} and \ref{fig:eigenmodes_c} we can see some of the first eigenmodes for these systems.
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/eigenmodes_s.png}
+\caption{Some of the first eigenmodes calculated for a square domain sorted by their frequencies.}
+\label{fig:eigenmodes_s}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=9cm]{Pictures/eigenmodes_r.png}
+\caption{Some of the first eigenmodes calculated for a rectangular domain sorted by their frequencies.}
+\label{fig:eigenmodes_r}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=9cm]{Pictures/eigenmodes_c.png}
+\caption{Some of the first eigenmodes calculated for a circular domain sorted by their frequencies.}
+\label{fig:eigenmodes_c}
+\end{figure}
+
+In figure \ref{fig:L_dependence} we see how the eigenfrequencies depend on the size of the domain. We see it corresponds to a $\frac{1}{L^2}$ dependence, which makes sense since our domain area is of size $L^2$.
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/L_dependence.png}
+\caption{Dependence of the eigenfrequency on the size of the domain for three different shapes. A reference line is plotted to show the $1/L^2$ dependency.}
+\label{fig:L_dependence}
+\end{figure}
+
+In figure \ref{fig:dstep_dependence} we plotted the relation between
+the number of integration steps and the eigenfrequency. In theory this should not be of any influence since the system itself does not change, however due to limited resolution at low step numbers we get a numerical error in the eigenvalue. We see how the calculated eigenfrequency forms an assymptote which we can assume to be the true eigenfrequency of the system.
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/dstep_dependence.png}
+\caption{Dependence of the eigenfrequency on the number of discretization steps}
+\label{fig:dstep_dependence}
+\end{figure}
+
+The time evolution calculated by equation \ref{eq:timedepsol} is shown in figure \ref{fig:time_evolution_s}, \ref{fig:time_evolution_r} and \ref{fig:time_evolution_c} for different eigenmodes and several timesteps.
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/time_evolution_s.png}
+\caption{Some of the first time evolution calculated for a square domain sorted by their frequencies.}
+\label{fig:time_evolution_s}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/time_evolution_r.png}
+\caption{Some of the first time evolution calculated for a rectangular domain sorted by their frequencies.}
+\label{fig:time_evolution_r}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{Pictures/time_evolution_c.png}
+\caption{Some of the first time evolution calculated for a circular domain sorted by their frequencies.}
+\label{fig:time_evolution_c}
+\end{figure}
+
+
+\section{Schr\"odinger's equation}
+
+The time-independent Sch\"odinger equation for a single particle in a potential $V(x)$ is:
+
+\begin{equation} \label{eq:schrodinger}
+ -\frac{\hbar^2}{2m} \nabla^2 \psi(x) + V(x) \psi(x) = E \psi(x).
+\end{equation}
+
+This is an eigenvalue problem: there are infinitely many solutions $\psi_n(x)$, each with a different energy $E_n$. But the distinct $E_n$ are discrete values, the particle cannot have every energy. We will calculate the wave equations $\psi_n(x)$ and their energies $E_n$ by discretizing eq. \ref{eq:schrodinger}.
+
+Setting $\hbar = m = 1$ and using a centered finite difference method for the second spatial derivative, we find:
+
+\begin{equation} \label{eq:discrete_schrodinger}
+ (2 + 2 \Delta x^2 V(x)) \psi_i - \psi_{i-1} - \psi_{i+1} = 2 \Delta x^2 E \psi_i.
+\end{equation}
+
+We will first look at the infinite potential well, that is defined as:
+
+\begin{equation} \label{eq:inf_pot}
+ V(x) = 0 \text{ if} -a \leq x \leq x, \text{ infinite otherwise}
+\end{equation}
+
+We discretize space from $-a$ to $a$ with the points $i = 1, 2, \dots, N$, where $\Delta x = 2 a/ N$ and $x(i) = -a + i \Delta x$. We assume that $\psi(0) = \psi(N) = 0$, this models infinite potential barriers at the boundaries.
+
+We can write eq. \ref{eq:discrete_schrodinger} explicitly for all $i = 1, 2, \dots, N$ and derive a $N \times N$ matrix:
+
+\begin{equation} \label{eq:schrodinger_eigenvalue}
+\left[
+\begin{array}{ccccccc}
+\phi & -1 & 0 & 0 & 0 & \dots & 0\\
+-1 & \phi & -1 & 0 & 0 & \dots & 0\\
+0 & -1 & \phi & -1 & 0 & \dots & 0\\
+\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\
+0 & \dots & 0 & 0 & -1 & \phi & -1\\
+0 & \dots & 0 & 0 & 0 & -1 & \phi
+\end{array}
+\right]
+\begin{bmatrix*}
+\psi_1 \\
+\psi_2 \\
+\vdots \\
+\psi_N \\
+\end{bmatrix*}
+=
+2\Delta x^2 \vect{E}
+\begin{bmatrix*}
+\psi_1 \\
+\psi_2 \\
+\vdots \\
+\psi_N \\
+\end{bmatrix*},
+\end{equation}
+
+where $\phi = 2 + 2 \Delta x ^2 V(x)$. Note that in the first and last row we've used that $\psi_0 = \psi_N = 0$. In this form, the problem can be readily solved using standard routines to find the eigenvalues and eigenvectors.
+
+We've solved this for $N = 200$ and the infinite potential well as defined in eq. \ref{eq:inf_pot}, the computed probability distribution functions (the eigenvectors) can be seen in figure \ref{fig:infinite} in green. In blue are the real, analytic solutions. We see that they are very close. The potential is invisible because it is everywhere zero. The analytic solutions were computed using\footnote{adapted from equation 2.28 in \cite{griffiths}}:
+
+\begin{equation} \label{eq:inf_pot_analytic}
+ \psi_n(x) = \sqrt{\frac{2}{a}} \sin(n \pi (x - a) / (2 a)),
+\end{equation}
+where $n$ is the quantum number. Increasing numbers of $n$ belong to increasing energies $E_n$. For the infinite potential well we also have an analytic expression for the allowed energies\footnote{from equation 2.27 in \cite{griffiths}}:
+
+\begin{equation} \label{eq:inf_pot_en}
+ E_n = \frac{n^2 \pi^2 \hbar^2}{2m(2a)^2}.
+\end{equation}
+
+\begin{figure}
+\centering
+\includegraphics[width=6cm]{infinite_1}
+\includegraphics[width=6cm]{infinite_2}
+\includegraphics[width=6cm]{infinite_3}
+\includegraphics[width=6cm]{infinite_4}
+\caption{The solutions to the time-independent Schr\"odinger equation computed analytically and by solving the eigenvalue problem from eq. \ref{eq:schrodinger_eigenvalue}.}
+\label{fig:infinite}
+\end{figure}
+
+In figure \ref{fig:energies} we see the corresponding eigenvalues of the solutions in figure \ref{fig:infinite} with the analytic values calculated using eq. \ref{eq:inf_pot_en}. We see that for low eigenvalues till 50 the calculated and real values mostly overlap. From then on they start to diverge and the calculated energies reach a plateau. It seems that here only the first $N/4$ calculated energy values can be trusted.
+
+\begin{figure}
+\centering
+\includegraphics[width=10cm]{energies}
+\caption{The calculated and analytic energies for a particle in an infinite potential well.}
+\label{fig:energies}
+\end{figure}
+
+Let's now look at a finite potential well that is defined as:
+
+\begin{equation}
+ V(x) = 0 \text{ if} -a \leq x \leq a, \text{ $V_0$ otherwise},
+\end{equation}
+with $V_0 = 10$. The only difference in the implementation are the values for $\phi$ in the first and last rows of \ref{eq:schrodinger_eigenvalue}. The energies are a bit higher than for the infinite potential well and the wave equations are very similar to those in figure
+\ref{fig:infinite}. Because these results are rather uninteresting, we will magnify the effect of the finite potential by adjusting its domain:
+
+\begin{equation} \label{eq:finite_pot}
+ V(x) = 0 \text{ if} -4a/5 \leq x \leq 4a/5, \text{ $V_0$ otherwise},
+\end{equation}
+again with $V_0=10$. Now more rows in matrix \ref{eq:schrodinger_eigenvalue} model the effect of the finite potential well. The resulting wave functions can be seen in figure \ref{fig:finite}. We observe tunneling: there is a non-zero probability of the particle crossing the potential barrier.
+
+\begin{figure}
+\centering
+\includegraphics[width=6cm]{finite_scaled_1}
+\includegraphics[width=6cm]{finite_scaled_2}
+\includegraphics[width=6cm]{finite_scaled_3}
+\includegraphics[width=6cm]{finite_scaled_4}
+\caption{The solutions to the time-independent Schr\"odinger equation found by solving the eigenvalue problem from eq. \ref{eq:schrodinger_eigenvalue} with the finite potential well defined in eq. \ref{eq:finite_pot}. The potential is displayed in purple.}
+\label{fig:finite}
+\end{figure}
+
+Lastly we look at a parabolic potential well defined as:
+\begin{equation} \label{eq:finite_pot}
+ V(x) = bx^2,
+\end{equation}
+where $b$ specifies the magnitude of the potential. In figure \ref{fig:parabolic} we see the wave equations for a particle in such potential with $b=500$. When comparing with the wave equations in figure \ref{fig:infinite}, we see that the parabolic potential pushes the particle to the middle, where the potential and thus the energy of the particle would be lowest.
+
+The energies for both the finite and parabolic potentials are a bit higher than for the infinite potential. We don't explicitly show them because they would almost coincide with the infinite potential will as calculated in \ref{fig:energies}.
+
+\begin{figure}
+\centering
+\includegraphics[width=6cm]{parabolic_1}
+\includegraphics[width=6cm]{parabolic_2}
+\includegraphics[width=6cm]{parabolic_3}
+\includegraphics[width=6cm]{parabolic_4}
+\caption{Solutions to the time-independent Schr\"odinger equation for a particle in a parabolic potential well.}
+\label{fig:parabolic}
+\end{figure}
+
+\section{Direct methods for solving steady state problems}
+
+Here we look at solving the diffusion equation for the steady state using a direct method. The diffusion equation is:
+\begin{equation*} \label{eq:diffusion}
+ \frac{\partial c}{\partial t} = \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right) c,
+\end{equation*}
+for steady states we have $\partial c/ \partial t = 0$, so we have to solve $\nabla^2 c = 0$. We construct a matrix $M$ for the $\nabla^2$ as described in the appendix. We take a domain of 4 by 4, and discretize the $x$ and $y$ directions with $N=80$ equal units. The $M$ matrix then is a $N^2$ by $N^2$ matrix. For the boundary we construct a vector $b$ that has $N^2$ elements, where all values are zero, except the point (0.6, 1.2) which is 1. The solution for this system is shown in figure \ref{fig:diffusion}. We tried applying the circular domain by setting those pixels in $M$ that are beyond a radius of 2 to zero, but that failed because $M$ would then be a singular matrix with no solution.
+
+\begin{figure}
+\centering
+\includegraphics[width=14cm]{diffusion}
+\caption{The solution to the diffusion equation. The domain for both $x$ and $y$ directions are from -2 to 2.}
+\label{fig:diffusion}
+\end{figure}
+
+\newpage
+\appendix
+\section{Matrix representation of the wave equation}
+\label{app:A}
+
+The discretization of the two dimensional wave equation can be written in matrix form. The general form for any square systemsize is given as follows:
+\begin{equation}
+\left[
+\begin{array}{c:c:c:c:c}
+ A &I & 0 & \dots & 0 \\
+ \hdashline
+ I &A &I & \ddots & \vdots \\
+ \hdashline
+ 0 &I &\ddots & I & 0\\
+ \hdashline
+ \vdots &\ddots &I & A & I \\
+ \hdashline
+ 0 &\dots &0 & I & A \\
+\end{array}
+\right ]
+\begin{bmatrix*}
+ v_{1,1}\\
+ v_{2,1}\\
+ v_{1,2}\\
+ \vdots\\
+ v_{N,N}
+\end{bmatrix*}
+=
+\lambda
+\begin{bmatrix*}
+ v_{1,1}\\
+ v_{2,1}\\
+ v_{1,2}\\
+ \vdots\\
+ v_{N,N}
+\end{bmatrix*}
+\label{eq:genMatExp}
+\end{equation}
+
+Where I is the identity matrix of size $N \times N$, 0 is the zero matrix of the same size and A is an $N \times N$ matrix of the following form:
+
+\begin{equation}
+A =
+\left[
+\begin{array}{rrrrr}
+ -4 &1 & 0 & \dots & 0 \\
+ 1 &-4 &1 & \ddots & \vdots \\
+ 0 &1 &\ddots &\ddots& 0\\
+ \vdots &\ddots &\ddots & -4 & 1 \\
+ 0 &\dots &0 & 1 & -4 \\
+\end{array}
+\right]
+\label{eq:matA}
+\end{equation}
+
+For example the full finite difference matrix in a $4 \times 4$ domain looks like this:
+
+\begin{equation}
+\left[
+\tiny
+\begin{array}{rrrr|rrrr|rrrr|rrrr}
+ -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
+ 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
+ 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
+ 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
+\hline
+ 1 & 0 & 0 & 0 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
+ 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
+ 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
+ 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
+ \hline
+ 0 & 0 & 0 & 0 &1 & 0 & 0 & 0 & -4 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
+ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\
+ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\
+ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 0 & 0 & 0 & 1 \\
+ \hline
+ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &1 & 0 & 0 & 0 & -4 & 1 & 0 & 0 \\
+ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 \\
+ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 \\
+ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & -4 \\
+
+\end{array}
+\right]
+\label{eq:mat4x4}
+\end{equation}
+
+\bibliography{references}
+\bibliographystyle{unsrt}
+
+\end{document}
diff --git a/Set 3/Report/references.bib b/Set 3/Report/references.bib
new file mode 100644
index 0000000..5cdaf7f
--- /dev/null
+++ b/Set 3/Report/references.bib
@@ -0,0 +1,70 @@
+@misc{wiki:lj,
+ author = "Wikipedia",
+ title = "{Lennard-Jones potential --- Wikipedia{,} The Free Encyclopedia}",
+ year = "2015",
+ note = "[{\url{https://en.wikipedia.org/w/index.php?title=Lennard-Jones_potential&oldid=694525432}}; accessed 27-January-2016]"
+}
+
+@misc{scipy,
+ author = "Scipy developers",
+ title = "{scipy.org}",
+ year = "2016",
+ note = "[{\url{https://www.scipy.org/}}; accessed 22-May-2016]"
+}
+
+
+@misc{davies,
+ author = "{Jason Davies}",
+ title = "{Diffusion-limited aggregation}",
+ note = "[{\url{https://www.jasondavies.com/dla/}}; accessed 10-May-2016]"
+}
+
+@unpublished{scico1,
+author = "{J. Camphuijsen, R. Kasim}",
+title = "{Scientific Computing 1 - Exploring the SciCoPyth's mind}",
+year = {2016}
+}
+
+@book{griffiths,
+ author = "David J. Griffiths",
+ title = "{Introduction to Quantum Mechanics}",
+ publisher = "Pearson Education",
+ year = "2005"
+}
+
+@book{numMeth,
+ author = "William F. Ames",
+ title = "{Numerical Methods for Partial Differential Equations}",
+ publisher = "Thomas Nelson and Sons LTD",
+ year = "1969"
+}
+
+@book{sciComp,
+ author = "Micheal T. Heath",
+ title = "{Scientific Computing: An Introductory Survey}",
+ publisher = "Elizabeth A. Jones",
+ edition = "2nd",
+ year = "2002",
+ ISBN = "0-07-112229-x"
+}
+
+
+@misc{sunsistemoGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/sunsistemo}}; accessed 12-April-2016]"
+}
+
+@misc{sunsistemo,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{Sunsistemo - The N-Body Simulator}",
+ note = "[{\url{https://sunsistemo.js.org}}; accessed 12-April-2016]"
+}
+
+@misc{SciCoPythGH,
+ author = "{Jaro Camphuijsen, Rahiel Kasim}",
+ title = "{SciCoPyth on GitHub}",
+ note = "[{\url{https://github.com/sunsistemo/SciCoPyth}}; accessed 12-April-2016]"
+}
+
+
diff --git a/Set 3/Scientific_Computing_3.pdf b/Set 3/Scientific_Computing_3.pdf
new file mode 100644
index 0000000..5fb3e76
--- /dev/null
+++ b/Set 3/Scientific_Computing_3.pdf
Binary files differ