\documentclass[11pt]{article}
\usepackage{amssymb,amsmath,amsthm,graphics,subfigure,graphicx}
\usepackage[latin1]{inputenc} 
\usepackage[english]{babel}     
\usepackage[T1]{fontenc}      
\usepackage{microtype}        
\usepackage{a4}
\usepackage{graphicx}
\usepackage{block}
\usepackage{paralist}
\usepackage{url}
\usepackage{hyperref}

%%% Punkter 
\newcounter{Punkt}
\renewcommand{\thePunkt}{\alph{Punkt}}
\newcommand{\Punktkolon}{)}
\newcommand{\PunktInn}{2em}
\newenvironment{punkt}%
 {\refstepcounter{Punkt}%
  \list{}{\setlength{\topsep}{0pt}}
  \item[{\bf \thePunkt\Punktkolon}]}
 {\endlist}

\newcommand{\dd}{\mathrm{d}}
\DeclareMathOperator{\Var}{Var}

\begin{document}
\begin{center}
    \textbf{\Large TMA4215 Numerical Mathematics} \\[5mm]
    Autumn 2009 \\[5mm]
    \textbf{Project 2} 
\end{center} 

Read the paper \textit{An Algorithmic Introduction to Numerical Simulation of Stochastic Differential equations} by Desmond J.\ Higham which can be found on the home page. Have focus on:
\begin{compactitem}
\item What is a stochastic differential equation? Where do they come from? 
\item How do we solve them numerically? Are there particular things to be aware of when solving SDEs?
\item How can you say something about the quality of the method (error, order and stability)?
\end{compactitem}
The \textsc{Matlab} files can be downloaded from \url{http://fox.maths.strath.ac.uk/~aas96106/algfiles.html}. 
When studying the \textsc{Matlab} code, you might find commands you are unfamiliar with. Use \texttt{doc} to get an explanation. It is strongly recommend that you run the files as you read the paper.

Rather than using the Milstein scheme presented in the paper, you will in the term project apply the following Runge--Kutta scheme:
\begin{align*}
   Y &= X_{j-1} + \sqrt{\Delta t}\, g(X_{j-1}) \\
   X_{j} &= X_{j-1} + \Delta t f(X_{j-1}) + \Delta W_j g(X_{j-1}) + 
   \frac{g(Y)-g(X_{j-1})}{2\sqrt{\Delta t}}\bigl((\Delta W_j)^2-\Delta t\bigr)
\end{align*} 
with $\Delta W_j=W(t_{j})-W(t_{j-1})$. 

\subsection*{Tasks}
\begin{compactitem}
\item Rewrite the \textsc{Matlab} file \texttt{emstrong.m} to justify that the Runge--Kutta method is of strong order~1.0.
\item Apply the method to the test problem assigned your group. 
\item Run 100 simulations. Plot 5 of the solution paths, as well as the mean of all the solution paths. Compare with the deterministic solution. What do you observe? 
\item \textbf{Extra:} If you have some spare time, see if you can do a linear mean square stability analysis of the Runge--Kutta method. (You will not lose any points by not doing this part, but you may earn some extra. The total number of points will not exceed 15.)
\end{compactitem}

\subsection*{Report}
Hand in a report consisting of maximum 6 pages (7 if you include a stability analysis).
Try to write the report as a small research paper, not as an assignment. 
The report should include a minimum of theoretical background in addition to the results of the tasks described above. You may include a few lines of \textsc{Matlab} code to emphasise a particular point. Otherwise, code listings should be avoided. If you use material from other sources, make sure to include references.  

\textit{Do not include your \textsc{Matlab} programs, not as an appendix, nor as separate files.} 
\newpage
\subsection*{Test equations}
Use the step-size $\Delta t=2^{-8}$ as default. You may want to experiment with the parameter~$\sigma$.
\begin{punkt}\textbf{Bonhoeffer--van der Pol oscillator}\\
This is simplification of the more famous Hodgkin and Huxley model for firing of a single neuron. The equation is given by
\begin{align*}
\dd X_1(t) &= c\bigl(X_1(t)+X_2(t)-\tfrac{1}{3}X_1(t)^3+z\bigr)\dd t+\sigma \dd W(t) & X_1(0) &= -1.9 \\
\dd X_2(t) &= -\tfrac{1}{c}\bigl(X_1(t)+bX_2(t)-a\bigr)\dd t, & X_2(0)&=1.2
\end{align*}  
Solve it  over the interval $0\leq t\leq 20$. 
Use $a=0.7$, $b=0.8$, $c=3.0$, and $z=-0.34$. Use $\sigma=0$ (deterministic) and $\sigma=0.5$.
\end{punkt}


\begin{punkt}\textbf{The Lotka--Volterra system}\\
This is the predatory-prey problem in which we have accepted a fluctuation in one of the growth parameters. The equation is given by
\begin{align*}
\dd X_1(t) &= \bigl(1-aX_2(t)\bigr)X_1(t)\dd t + \sigma X_1(t) \dd W(t) & X_1(0)& =0.8 \\
\dd X_2(t) &= \bigl(bX_1(t)-c\bigr)X_2(t)\dd t  & X_2(0)&= 1.0
\end{align*}
Solve the system over $0\leq t\leq 10$ and use the parameters $a=b=1$, $c=3$. Use $\sigma = 0$ and $\sigma = 0.2$. 
\end{punkt}

\begin{punkt}\textbf{Untitled}
\begin{align*}
\dd X_1(t)&=X_2(t) \dd t & X_1(0)&=1.0\\
a\, \dd X_2(t)&=\bigl(-X_1(t)-b\sin X_1(t) + \omega\bigr)\dd t+\sigma\dd W(t), &X_2(0)&=0.
\end{align*}
Let $0\leq t \leq 5$. Use $a=0.1$, $b = 1.0$, $\omega=\pi$. Try with $\sigma=0.0$ and $\sigma = 0.2$.  
\end{punkt}

\newpage
\begin{center}
\Large \bf Emergency Sheet
\end{center}
The \emph{expectation} or \emph{mean value} of a random variable $X$ with distribution $p(x)$ is given by
\[ \mu = \mathbb{E} X = \int_{-\infty}^{\infty}x p(x)\,\dd x \]
and the \emph{variance} is
\[ \sigma^2 = \Var(X) = \int_{-\infty}^{\infty} (x-\mu)^2p(x)\,\dd x. \]
A random variable is \emph{normally} or \emph{Gaussian} distributed if
\[  p(x) = \frac{1}{\sqrt{2\pi}\sigma} \mathrm{e}^{-(x-\mu)^2/(2\sigma^2)}. \]
A normal distributed random variable $X$ with mean $\mu$ and variance $\sigma^2$ is often denoted as $ X \sim N(\mu,\sigma^2)$. For normal distributed variables we get
\[ \Var(X)=\mathbb{E} X^2 - (\mathbb{E} X)^2. \]

\end{document}
