% $Id: numModellingTheory.tex 1339 2007-05-09 05:07:30Z joehope $

% Copyright (C) 2000-2007
%
% Code contributed by Greg Collecutt, Joseph Hope and the xmds-devel team
%
% This file is part of xmds.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software 
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

\begin{latexonly}
\begin{savequote}[8cm]
%\href{http://nearingzero.net/}{http://nearingzero.net/}
\includegraphics[width=80mm]{figures/defferential}
\end{savequote}
\end{latexonly}

\chapter{Numerical Modelling Theory} 
\label{chap:numericalModelling}

\section{Differential Equations} 
\label{sec:differentialEquations}

The majority of numerical models for real world problems are based on
one or more differential equations involving the various parameters
concerned. The large majority of which may be written in the very
general form:
\begin{equation}
\left[ A\left(x,y\right) \frac{\partial^2}{\partial x^2} +
       2B\left(x,y\right) \frac{\partial^2}{\partial x \partial y} +
       C\left(x,y\right) \frac{\partial^2}{\partial y^2} \right]
       a^i\left(x,y\right) = f^i\left(x,y,\vect{a},\frac{\partial
       a^i}{\partial x},\frac{\partial a^i}{\partial y}\right).
\end{equation}
These come under three sub-classifications~\cite{Kreyszig:1999:1}
according the the parameters $A$, $B$, and $C$:
\begin{enumerate}
\item The PDE is said to be {\em elliptic} in regions of the $x$-$y$
plane where $AC-B^2 > 0$.
\item The PDE is said to be {\em parabolic} in regions of the $x$-$y$
plane where $AC-B^2 = 0$.
\item The PDE is said to be {\em hyperbolic} in regions of the $x$-$y$
plane where $AC-B^2 < 0$.
\end{enumerate}

This model for classifying PDEs is only two dimensional, but many
problems in physics exhibit up to four dimensions. However, the
majority of PDEs in physics are able to be expressed in a form where
one of the dimensions is special; special in that there are only ever
first order derivatives involving that dimension. A two dimensional
form of such a PDE looks like:
\begin{equation}
  C\left(x,y\right) \frac{\partial^2}{\partial y^2} a^i\left(x,y\right)
  = f^i\left(x,y,\vect{a},\frac{\partial a^i}{\partial
    x},\frac{\partial a^i}{\partial y}\right),
\end{equation}
and such an equation falls into the category of ``parabolic'' by the
above classification. For this reason many three and four dimensional
PDEs in physics are referred to as parabolic, although technically
this classification exists only for PDEs in two dimensions.  The more
general form for such ``parabolic-like'' differential equations is:
\begin{equation}
  \frac{\partial}{\partial x^0} a^i(\vect{x}) =
  f^i\left[\vect{x},\vect{a}(\vect{x})\right],
\end{equation}
where the functionals $f^i$ may include first or second order partial
derivatives of the {\em transverse} dimensions, $x^{i\neq0}$.  This
set of equations can be generalised to allow the functionals $f^i$ to
include partial derivatives of any order. Given that the equations are
in this general form, they may be propagated in the first order
dimension, $x^0$, using one of the algorithms detailed further on.

\subsection{Boundary Conditions} 
\label{sec:boundaryConditions}

In the case of partial differential equations, the equations must
apply either over all of the space spanned by the transverse
dimensions, or else over a confined region of this space. Either way,
computational resources cannot solve over an infinite region of space,
and so boundaries of some sort must be imposed. Therefore, a
description of what happens at these boundaries is essential for
solving the problem. There are three common types of boundary
conditions:
\begin{enumerate}
\item {\bf Dirichlet boundary conditions} This is where the value of a
component is fixed at a particular boundary:
\begin{equation}
  a^i\left(\vect{x} | x^k=x^k_{\pm}\right) = a^i_{k\pm}.
\end{equation}
\item {\bf Neumann boundary conditions} This is where the
cross-boundary gradient of a component is fixed at a particular
boundary:
\begin{equation}
  \frac{\partial}{\partial x^k} a^i\left(\vect{x} | x^k=x^k_{\pm}\right)
  = b^i_{k\pm}.
\end{equation}
\item {\bf Periodic boundary conditions} This is where the value of a
component (and all of its cross-boundary derivatives) is the same at
each end of the dimension in question.  It is as if this dimension was
originally circular, except that it has been cut for the purpose of
discretisation.
\begin{equation}
  \left(\frac{\partial}{\partial x^k}\right)^n a^i\left(\vect{x} |
  x^k=x^k_{+}\right) = \left(\frac{\partial}{\partial
    x^k}\right)^n a^i\left(\vect{x} | x^k=x^k_{-}\right), \quad
  n=0,1,2...
\end{equation}
\end{enumerate}
If an infinite region of space is absolutely necessary, then it can be
mapped to a finite region of space using a transformation such as
$x=\tan(\epsilon)$, which maps an infinite range in $x$ to a finite
range, $[-\frac{\pi}{2},+\frac{\pi}{2}]$, in $\epsilon$. There are,
of course, a few caveats when performing such a mapping.

Finally, the last boundary condition is the initial state of the
components. The values (or form) of all the components must be
specified for some point in the dimension of propagation, $x^0$.

\section{Stochastic Equations} 
\label{sec:stochasticEquations}

Stochastic differential equations also include noise sources,
$\vect{\xi}(\vect{x})$, which may be real or complex:
\begin{equation}
  \frac{\partial}{\partial x^0} a^i(\vect{x}) =
  f^i\left[\vect{x},\vect{a}(\vect{x})\right] +
  g^{ij}\left[\vect{a}(\vect{x})\right]\xi^j(\vect{x}).
\end{equation}
These can make propagation a little more technically challenging, but
not impossible, and enables a vast range of complex phenomena to be
modelled. As a result an individual propagation path becomes dependent
on the underlying noise sources used, and usually many such paths
using independent noises are combined and averaged to obtain the
desired result. The noise sources are normally delta-correlated in the
propagation dimension, though not always in the transverse dimensions,
and not always with each other, so that generally:
\begin{align}
  \left<\xi^i(\vect{x})\xi^j(\vect{x}')\right> & =
  \delta(x^{0'}-x^0)C_{ij}(\vect{x}'_\bot,\vect{x}_\bot);\\
  \left<\xi^i(\vect{x})\xi^{j*}(\vect{x}')\right> & =
  \delta(x^{0'}-x^0)D_{ij}(\vect{x}'_\bot,\vect{x}_\bot),
\end{align}
where as usual the subscript $\bot$ denotes the transverse dimensions.

\subsection{Ito vs Stratonovich Calculus} 
\label{sec:itoStratonovichCalculus}

The integration of a stochastic equation as expressed above is not
well defined in normal calculus. In fact, there are multiple ways in
which such an integration can be defined as the limit of some discrete
process, and they are not equivalent. The two main forms of stochastic
calculus are Ito calculus and Stratonovich calculus. Ito calculus is
better suited to mathematical proofs, but Stratonovich calculus has
the major advantage that derivatives follow the normal chain
rule. This typically simplifies the numerical technique for solving a
Stratonovich equation with a high-order method. Both types of
equations arise naturally from modelling physical systems.
Fortunately, it is possible to transform a set of equations in Ito
calculus to Stratonovich calculus (and vice-versa).

If a set of Ito equations is given by:
\begin{equation}
(I) \qquad \frac{\partial}{\partial x^0} a^i =
f^i\left[\vect{a}\right] + g^{ij}\left[\vect{a}\right]\xi^j,
\end{equation}
then the corresponding set of Stratonovich equations is
\begin{equation}
(S) \qquad \frac{\partial}{\partial x^0} a^i =
f^i\left[\vect{a}\right]-\frac{1}{2}\sum_{jk}
g^{kj}\left[\vect{a}\right] \frac{\partial}{\partial a^k}
g^{ij}\left[\vect{a}\right] +
g^{ij}\left[\vect{a}(\vect{x})\right]\xi^j(\vect{x}).
\end{equation}
Numerical integration of stochastic equations is more computationally
intensive than deterministic equations of equivalent size and
complexity.

\section{Numerical Methods for Differential Equations} 
\label{sec:numericalMethods}

We now briefly describe several numerical methods for integrating
deterministic and stochastic differential equations.

\subsection{The Euler and Inverse Euler Methods} 
\label{sec:eulerMethods}

Consider the ordinary differential equation:
\begin{equation}
\frac{\partial}{\partial t} a(t) = f[t,a],
\end{equation}
where f is some functional of $t$ and/or $a$. Given that a ``current''
value of $a$ (say $a_n = a(t_n)$) is known then the simplest method to
propagate $a$ forward in time is to use the current slope
($f[t_n,a_n]$) to calculate the new value $a_{n+1}$:
\begin{equation}
a_{n+1} = a_n + h \: f[t_n,a_n],
\end{equation}
where $h$ is the time step so that $t_{n+1} = t_n + h$. This is known
at the explicit, or {\em forward time}, Euler method. It is very
simple, but not very accurate or stable. Stability can be improved to
an extent by thinking in reverse. Rather than using the slope at the
current point in time to calculate the next value of $a$, use the
slope and the {\em next} point in time. But this, of course, is not
known until the next value of $a$ is known, so the process becomes an
iterative, or {\em implicit} one. Initially the current value of $a$
is used as next value, and then the step:
\begin{equation}
a_{n+1} = a_n + h \: f[t_{n+1},a_{n+1}],
\end{equation}
is iterated until convergence is obtained. This method is known as the
{\em implicit}, or {\em backward time}, or {\em inverse}, Euler
method. For certain classes of differential equation, implicit methods
are more stable than their explicit equivalents, but overall the
implicit Euler method is no more accurate than the explicit Euler
method.

The explicit and implicit Euler methods do not account for the second
and higher order derivatives of $a$. Thus the error per step is of
order $h^2$, which means the error over a given integration interval
is of order $h$. Consequently these methods are known as first order
methods.

Both of these methods can be used to integrate stochastic equations by
treating the noise terms as part of a larger function $f$ which
depends on a pseudo-random number which changes each time step.  This
technique converges to the Ito integral with a convergence of order
$h^{1/2}$.  In general, the order of the stochastic convergence of a
numerical method will be less than or equal to half the order of the
deterministic convergence.

\subsection{The Improved Euler Method} 
\label{sec:improvedEulerMethod}

Here two evaluations of the derivative $f[t,a]$ are made. The first
one is based on the current value $a$, and is used to estimate the
next value of $a$ as in the explicit method. This estimated next value
is then used to calculate the next value of the derivative, which in
turn is used to make a second estimate of the next value of $a$ as in
the implicit method. These two estimates are then averaged:
\begin{align}
a^{(1)}_{n+1} & = a_n + h \: f[t_n,a_n];\\
a^{(2)}_{n+1} & = a_n + h \: f[t_{n+1},a^{(1)}_{n+1}];\\
a_{n+1} & = \frac{1}{2}\left(a^{(1)}_{n+1} + a^{(2)}_{n+1}\right).
\end{align}
The two estimates of $f[t,a]$ account for the second order derivative
of $a$, so the error per step is reduced to order $h^3$, and the
integral is correct to order $h^2$.  This method is also known as the
improved Euler-Cauchy method, or Huen's method. It is also known as
the second-order Runge-Kutta method.

Unfortunately, this method is not guaranteed to converge for
stochastic integration when the noise terms are naively added to the
function on the RHS.

\subsection{The Semi-Implicit Method} 
\label{sec:siMethod}

This method is very similar to the improved Euler method. Essentially
the method lies half way between the explicit and implicit methods in
that the new value for $a$ is calculated based on the derivative at
the mid-point of the step. Similarly to the implicit method, the
current value of $a$ is used initially as the mid-point value,
$a_m=a(t_m); t_m = t + h/2$, and then the step:
\begin{equation}
a_m = a_n + \frac{1}{2} \: h \: f[t_m,a_m],
\end{equation}
is iterated until convergence is obtained. The next value for $a$ is
then:
\begin{equation}
a_{n+1} = 2a_m - a_n.
\end{equation}

Similar to the improved Euler method, the semi-implicit method is also
a second order method. It has the added advantage that a na\"{\i}ve
inclusion of the noise terms in the function $f$ produces an algorithm
that converges to the Stratonovich integral with global error of the
order of $h$.

\subsection{The Fourth-Order Runge-Kutta Method} 
\label{sec:rk4Method}

The fourth-order Runge-Kutta method has been the workhorse for
numerical modelling for many years. The method uses four evaluations
of the derivative to account for the first, second, third, and fourth
derivatives---making it a fourth order method. The method is simply:
\begin{align}
\vect{k}_1 &= h \: \vect{f}[t_n,a_n]; \\ 
\vect{k}_2 &= h \:
\vect{f}[\frac{1}{2}(t_n + t_{n+1}),\vect{a}_n+\frac{1}{2}\vect{k}_1]; \\ 
\vect{k}_3 &= h \: \vect{f}[\frac{1}{2}(t_n +
t_{n+1}),\vect{a}_n+\frac{1}{2}\vect{k}_2]; \\ 
\vect{k}_4 &= h \: \vect{f}[t_{n+1},\vect{a}_n+\vect{k}_3]; \\ 
\vect{a}_{n+1} &= \vect{a}_n + \frac{1}{6}\left(\vect{k}_1 +
2\vect{k}_2 + 2\vect{k}_3 + \vect{k}_4\right). \\
\end{align}

Higher order methods are possible, but usually the extra accuracy does
not warrant the extra computation. In other words, this method usually
achieves the fastest evaluation for a given propagation accuracy.
Note that even an optimised implementation of this algorithm will
require about three times as much memory as a lower order method.

Although empirically there is a subclass of Stratonovich stochastic
problems that converge with this algorithm to second order, there is
no guarantee of convergence for stochastic equations.

\subsection{The Adaptive Fourth-Order Runge-Kutta Method}
\label{sec:ark45Method}

For problems that exhibit very diverse behavior of the solution, say with slow change for large parts of the propagation and rapid variations in others, it is very useful to employ an algorithm that can adjust its stepsize. This allows the simulation to speed through smooth uninteresting countryside in a few great strides but to tiptoe with many small steps through treacherous terrain.

The adaptive stepsize method implemented in \xmds in an embedded fourth-fifth order Runge-Kutta method (ARK45). It takes advantage of the fact that it is possible to combine certain six function evaluations that can result in a fifth order Runge-Kutta method in another way to give a result that is accurate to fourth order. The ARK45 algorithm calculates the forth and fifth order solution and uses the difference between them as an estimate of the current discretisation error.

The structure is: 
\begin{align}
\vect{k}_1 &= h \: \vect{f}[t_n,a_n]; \nonumber \\ 
\vect{k}_2 &= h \:
\vect{f}[t_n+a_{2} h,\vect{a}_n+b_{21}\vect{k}_1];\nonumber \\ 
\vect{k}_3 &= h \: \vect{f}[t_n + a_{3}h,\vect{a}_n+b_{31}\vect{k}_1+b_{32}\vect{k}_2]; \nonumber\\ 
\vect{k}_4 &= h \: \vect{f}[t_n + a_{4}h,\vect{a}_n+b_{41}\vect{k}_1+b_{42}\vect{k}_2+b_{43}\vect{k}_3];\nonumber\\ 
\vect{k}_5 &= h \: \vect{f}[t_n + a_{5}h,\vect{a}_n+b_{51}\vect{k}_1+b_{52}\vect{k}_2+b_{53}\vect{k}_3+b_{54}\vect{k}_4]; \nonumber\\ 
\vect{k}_6 &= h \: \vect{f}[t_n + a_{6}h,\vect{a}_n+b_{61}\vect{k}_1+b_{62}\vect{k}_2+b_{63}\vect{k}_3 +b_{64}\vect{k}_4+b_{65}\vect{k}_5];\nonumber \\ 
\vect{a}_{n+1} &= \vect{a}_n +c_{1}\vect{k}_1+c_{2}\vect{k}_2+c_{3}\vect{k}_3+c_{4}\vect{k}_4+c_{5}\vect{k}_5 +c_{6}\vect{k}_6+ {\cal O}(h^{6});\nonumber\\
\vect{a}^{*}_{n+1} &= \vect{a}_n +c_{1}^{*}\vect{k}_1+c_{2}^{*}\vect{k}_2+c_{3}^{*}\vect{k}_3+c_{4}^{*}\vect{k}_4+c_{5}^{*}\vect{k}_5 +c_{6}^{*}\vect{k}_6+ {\cal O}(h^{5}).
\end{align}
%
Here $\vect{a}_{n+1}$ is the fifth order- and $\vect{a}^{*}_{n+1}$ the fourth order solution. The coefficients that appear are listed in table \ref{tab:ck}. To adjust the time step \xmds calculates the relative difference between the two solutions $\Delta_{m}$. If partial differential equations are being solved it will determine the maximum of the relative errors of all grid points, omitting points where the function is less than $c_{thresh} \times M[a]$, where $c_{thresh}$ is a small threshold value (see chapter \ref{chap:languageRef}) and $M[a]$ is the peak value of the function across the grid. If multiple functions are present, each peak value is determined separately. 

\begin{table}\centering
\begin{tabular}{|c|c|ccccc|c|c|}
\hline
$i\rule[-2.3mm]{0cm}{7mm}$ &$\;a_i\;$ & & & $b_{ij}$ & & &$\;c_i\;$ & $\;c_i^*\;$ \\
\hline
 $\;1\rule[-2.3mm]{0mm}{8mm}\;$&  &  & & & & & $\frac{37}{378}$ & $\frac{2825}{27648}$ \\
2\rule[-2.3mm]{0cm}{7mm}& $\frac{1}{5}$ & $\frac{1}{5}$ & & & & & 0 & 0 \\
3\rule[-2.3mm]{0cm}{7mm}& $\frac{3}{10}$ & $\frac{3}{40}$ & $\frac{9}{40}$& & & &
 $\frac{250}{621}$ & $ \frac{18575}{48384}$\\
4\rule[-2.3mm]{0cm}{7mm}& $\frac{3}{5}$ & $\frac{3}{10}$ & $-\frac{9}{10}$&$\frac{6}{5}$ & & &
 $\frac{125}{594}$ & $\; \frac{13525}{55296}\;$\\
5\rule[-2.3mm]{0cm}{7mm}& $1$ & $-\frac{11}{54}$ & $\frac{5}{2}$&$-\frac{70}{27}$ &$\frac{35}{27}$ & &
 $0$ & $ \frac{277}{14336}$\\
6\rule[-3.3mm]{0cm}{8mm}& $\frac{7}{8}$ & $\frac{1631}{55296}$ & $\frac{175}{512}$&
$\frac{575}{13824}$ &$\frac{44275}{110592}$ & $\frac{253}{4096}$ &
 $\;\frac{512}{1771}\;$ & $ \frac{1}{4}$\\
\hline
\multicolumn{2}{|c|}{\rule[-2.3mm]{0cm}{7mm}j\;=} & 1& 2& 3 & 4& 5& \multicolumn{2}{c|}{}\\
\hline
\end{tabular}
\caption{Cash-Karp parameters for the embedded Runga-Kutta method.}
\label{tab:ck}\end{table}

If the prescribed accuracy is $ \Delta_{\rm tol}$, a time step of size $\delta \tau_{n}$ is accepted if $\Delta_m \leq \Delta_{\rm tol}$. It is rejected if  $\Delta_m > \Delta_{\rm tol}$, and the step must be calculated again. The size of the next time step $\delta \tau_{n+1}$ is determined by:
\begin{eqnarray}
\delta \tau_{n+1} = \left\{
\begin{array}{ll}
\displaystyle 0.92 \delta \tau_n \left|\frac{\Delta_{\rm
tol}}{\Delta_{m}}\right|^{1/5} & \mbox{if } \Delta_{m} \leq \Delta_{\rm tol},\\
\\
\displaystyle 0.92 \delta \tau_n 
\left|\frac{\Delta_{\rm tol}}{\Delta_{m}}\right|^{1/4}
& \mbox{if } \Delta_{m} >\Delta_{\rm tol}.\\
\end{array}
\right.\end{eqnarray}
%
This algorithm will consume even more memory than the fourth-order Runge-Kutta, but this easily pays off in situations where the rate of change of the solution varies a lot, or when the appropriate step size is unknown. 

As with the fourth-order Runge-Kutta  there is no guarantee of convergence for stochastic equations.

\section{Numerical Methods for Partial Differential Equations}
\label{sec:numericalMethodsPde}

The methods in the previous section all deal with ordinary
differential equations quite well, but partial differential equations
are considerably more complex, since the derivative depends also upon
the values for $a$ at other points in the transverse space.  In
principle, once the transverse spaces have been divided into a finite
lattice, sets of partial differential equations are just larger sets
of ordinary differential equations, but there are many subtleties
involved with achieving this discretisation.

The critical issue is how the evolution due to the transverse
derivatives is performed. This issue divides algorithms into two
classes: explicit picture and interaction picture. We begin here by
explaining just how the transverse derivatives may be calculated, and
then outline how the Semi-Implicit and the Runge-Kutta algorithms can
be applied in each of these cases.

\subsection{Evaluating Transverse Derivatives} 
\label{sec:transverseDerivatives}

The transverse dimensions must be discretised to form a
lattice. Lattices with more points enable more detail, but also take
up more memory, and more computational resource. The transverse
derivatives must be evaluated at each lattice point, and there must be
some method for doing so that uses the values of the field from the
other lattice points. There are two main methods for doing this.

Firstly, there is the {\bf matrix} method. If we let $a_k$ represent
the discretised values of $a(x)$ at the lattice points $x_k$, then we
may evaluate the first order derivative like so:
\begin{align}
\left(\frac{\partial a(x)}{\partial x}\right)_{x_k}
& = \frac{a_{k+1}-a_{k-1}}{x_{k+1}-x_{k-1}}\\
& = M_{jk}a_j, \quad k \neq 1,n.
\end{align}
Note the symmetry of the weightings about the point where the
derivative is being evaluated, causing the method to be {\em space
centered}. Note also that, unless periodic boundary conditions are
specified, it is not not possible to implement a space centered method
at the ends of the lattice. The best that can be done is to implement
a non-space centered method at the ends, and insist that the
derivatives vanish at the lattice boundaries for the technique to be
valid. {\bf M} may be expressed as a tri-diagonal square matrix, which
for periodic boundary conditions is:
\begin{align}
{\bf M} & = \frac{1}{2 \Delta x} \left[
\begin{array}{ccccccc}
 0 &  1 & 0 & 0 & ... & ... & -1\\
-1 &  0 & 1 & 0 & ... & ... & 0\\
 0 & -1 & 0 & 1 & ... & ... & 0\\
... & ... & ... & ... & ... & ... & ...\\
0 & ... & ... & -1 &  0 &  1 & 0\\
0 & ... & ... & 0  & -1 &  0 & 1\\
1 & ... & ... & 0  &  0 & -1 & 0\\
\end{array}
\right].
\end{align}
And similarly for evaluating the second derivative we might use:
\begin{align}
{\bf M} & = \frac{1}{(\Delta x)^2} \left[
\begin{array}{ccccccc}
-2 &  1 & 0  & 0 & ... & ... & 1\\
 1 & -2 & 1  & 0 & ... & ... & 0\\
0  &  1 & -2 & 1 & ... & ... & 0\\
... & ... & ... & ... & ... & ... & ...\\
0 & ... & ... &  1 & -2 &  1 & 0\\
0 & ... & ... &  0 &  1 & -2 & 1\\
1 & ... & ... &  0 &  0 &  1 & -2\\
\end{array}
\right].
\end{align}

These methods evaluate the derivatives using only the adjacent points
in space, but in principle it is possible to use progressively more
distant neighbours to produce an estimate of each derivative which
takes into account contributions from higher order derivatives.  In
practice, these methods are used rarely due to the availability of a
method based on the {\bf Fourier Transform}. The partial derivate with
respect to dimension $x^j$ is calculated quite simply using
\eqn{eq:l_in_k}:
\begin{align}
\frac{\partial}{\partial x^j}a(\vect{x}) & = \frac{\partial}{\partial
x^j} \frac{1}{\left(2\pi\right)^\frac{N}{2}}\int d\vect{k} e^{i
\vect{k}\cdot\vect{x}} \tilde{a}(\vect{k}) \nonumber \\ 
& =
\frac{1}{\left(2\pi\right)^\frac{N}{2}}\int d\vect{k}
e^{i\vect{k}\cdot\vect{x}} (ik^j) \tilde{a}(\vect{k}).
\label{eq:l_in_k}
\end{align}

This method can equally be applied to discrete data. Firstly a
discrete Fourier Transform is applied to the field data $a_k$, then
the data is multiplied point by point by the corresponding $ik^j$
values, and then the discrete inverse Fourier Transform is applied to
obtain the $a'_k$ data. Extending this method for higher order
derivatives is trivial.  \eqn{eq:l_maps} illustrates some
example mappings:
\begin{align}
f(k_x,k_y,...) & \mapsto f(k_x,k_y,...); \nonumber \\
\frac{\partial}{\partial x} & \mapsto ik_x; \nonumber \\
\frac{\partial^{2}}{\partial x^{2}} & \mapsto (ik_x)^2 = -k_x^{2};
\nonumber \\ 
\frac{\partial^{2}}{\partial x \partial y} & \mapsto
(ik_x)(ik_y) = -k_x k_y. 
\label{eq:l_maps}
\end{align}

This method becomes computationally more efficient than the matrix
method when a number of higher order derivatives are desired, for
example when a Laplacian operator is being applied in multiple
transverse dimensions. The main limitations of this method are that it
implies periodic boundary conditions and that it cannot be implemented
on a non-regular grid.

\subsection{Explicit Picture Methods} 
\label{sec:exMethod}

In the explicit picture, all linear and nonlinear operators are
explicitly calculated and summed to generate the total derivative for
each field component, whereas in the interaction picture linear
operators are dealt with separately to nonlinear operators.

Once a set of partial differential equations has been discretised to a
finite lattice in the transverse dimensions, it is formally equivalent
to a large set of coupled ordinary differential equations, and
therefore it can be solved by higher order methods such as the
semi-implicit method of \Sec{sec:siMethod}, or the fourth
order Runge-Kutta (RK4) method of \Sec{sec:rk4Method}.
Partial derivatives of the fields can be determined by Crank-Nicholson
methods or spectral methods. Thus explicit picture methods are capable
of solving the fully fledged generalised PDE:
\begin{align}
\frac{\partial}{\partial x^0}\vect{a}(\vect{x}) & =
  \vect{\mathcal{N}}\left(\vect{x},\;\vect{a}(\vect{x}),
  \vect{p}(\vect{x}),\;\vect{\xi}(\vect{x})\right),\\ 
  p^i(\vect{x}) & = \mathcal{F}^{-1}\left[\Sigma_j
  \mathcal{L}^{ij}\left(x^0,\vect{k}_{\bot}\right)
  \mathcal{F}\left[a^j(\vect{x})\right]\right].
  \label{eq:exPde}
\end{align}
This is the PDE that is listed in \eqn{eq:xmdsPdeEx} with
the cross propagating vector $\vect{b}$ omitted. Refer to
\Chap{chap:functionality} for a full explanation of the
variables. The important thing to note is that the total derivative of
the vector $\vect{a}$ is expressed in terms of the general nonlinear
functionals $\vect{\mathcal{N}}$, and that this form allows for
nonlinear partial derivatives that may have spatial dependence, for
example $\frac{1}{r}\left(\frac{\partial}{\partial r}\right)^2$.

\subsection{The Semi-Implicit method in the Explicit Picture} 
\label{sec:siexMethod}

Here we use the Fourier Transform (or Spectral) method to determine
the transverse derivatives. The procedure, after
\Sec{sec:siMethod}, is:
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a}_I = \vect{a}$
\item For N-1 iterations do:
  \begin{enumerate}
  \item $\vect{p} =
  \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
  \mathcal{F}\left[a\right]\right]$
  \item $\vect{a} = h \:
  \vect{\mathcal{N}}\left(\vect{x},\vect{a},\vect{p},\vect{\xi}\right)$
  \item $\vect{a} = \vect{a}_I + \frac{1}{2}\vect{a}$
  \end{enumerate}
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[a\right]\right]$
\item $\vect{a} = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a},\vect{p},\vect{\xi}\right)$
\item $\vect{a} = \vect{a}_I + \vect{a}$
\item $x^0 = x^0 + \frac{1}{2}h$
\end{enumerate}
Here an extra copy of the field is needed to store the initial
$\vect{a}$ into $\vect{a}_I$, and a similar size vector is needed for
the derivatives $\vect{p}$.

\subsection{The Fourth Order Runge-Kutta Method in the Explicit Picture} 
\label{sec:rk4exMethod}

Again we use the Fourier Transform (or Spectral) method to determine
the transverse derivatives. The procedure, after
\Sec{sec:rk4Method}, then is:
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
\item $\vect{a}_K = \vect{a}$
\item $\vect{a}_I = \vect{a}$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} + \frac{1}{6}\vect{a}_K$
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a}_K = \vect{a}_I + \frac{1}{2}\vect{a}_K$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} + \frac{1}{3}\vect{a}_K$
\item $\vect{a}_K = \vect{a}_I + \frac{1}{2}\vect{a}_K$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} + \frac{1}{3}\vect{a}_K$
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a}   = \vect{a}_I + \vect{a}$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} + \frac{1}{6}\vect{a}_K$
\end{enumerate}

This method requires two extra copies of the field. One, $\vect{a}_I$,
for the initial field value, and the other, $\vect{a}_K$, to act as a
working field while the original vector $\vect{a}$ collects the
derivative contributions. This saves the final step of copying
$\vect{a}_K$ back into $\vect{a}$ that would had to have been performed
if $\vect{a}_K$ had collected the derivative contributions.  It also
requires a similar size vector for the derivatives $\vect{p}$. On most
systems where the calculation is not memory-limited, this method tends
to achieve a desired accuracy many times faster than the semi-implicit
method.

\subsection{The Fourth/ Fifth Order adaptive Runge-Kutta Method in the Explicit Picture} 
\label{sec:ark45exMethod}

Once more we use the Fourier Transform (or Spectral) method to determine
the transverse derivatives. In order to minimize memory usage the current contents of the arrays containing the final solution are reused to calculate some of the k-vectors. The optimized procedure, after
\Sec{sec:ark45Method}, then is:
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
\item $\vect{a}_K = \vect{a}$
\item $\vect{a}_I = \vect{a}$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(x^0,\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} +c_{1}\vect{a}_K$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{1}^{*}\vect{a}_K$
\item $x^0 = x^0 + a_{2}h$
\item $\vect{a}_K = \vect{a}_I + b_{21}\vect{a}_K$
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}\left(x^0,\vect{a}_K,\vect{p},\vect{\xi}\right)$
\item $\vect{a}_J =(1-b_{31}/c_{1}) \vect{a}_I + b_{31}/c_{1}\vect{a} +b_{32} \vect{a}_K$
\item $x^0 = x^0 + (a_{3}-a_{2})h$ \\
	(Note: $c_{2}=c_{2}^{*}=0$) 
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_J = h \:
\vect{\mathcal{N}}\left(x^0,\vect{a}_J,\vect{p},\vect{\xi}\right)$
\item $\vect{a}_L =(1-b_{41}/c_{1}) \vect{a}_I + b_{41}/c_{1}\vect{a} +b_{42} \vect{a}_K +b_{43}\vect{a}_{J}$
\item $\vect{a}   = \vect{a} +c_{3}\vect{a}_J$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{3}^{*}\vect{a}_J$
\item $x^0 = x^0 + (a_{4}-a_{3})h$ \\
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_L = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_L,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} +c_{4}\vect{a}_L$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{4}^{*}\vect{a}_L$
\item $\vect{a}_L =(1-b_{51}/c_{1}) \vect{a}_I + b_{51}/c_{1}\vect{a} +b_{52} \vect{a}_K +(b_{53}-b_{51}c_{3}/c_{1})\vect{a}_{J} +(b_{54}-b_{51}c_{4}/c_{1})\vect{a}_{L} $
\item $x^0 = x^0 + (a_{5}-a_{4})h$ \\
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_L = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_L,\vect{p},\vect{\xi}\right)$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{5}^{*}\vect{a}_L$\\
	(Note: $c_{5}=0$) 
\item $\vect{a}_L =f_{1} \vect{a}_I  +f_{2} \vect{a}_K + f_{3}\vect{a}_{J} + f_{4} \vect{a} +f_{5}\vect{a}_{L} + f_{6} \vect{a}^{*}$		
\item $x^0 = x^0 + (a_{6}-a_{5})h$ \\
\item $\vect{p} =
\mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot
\mathcal{F}\left[\vect{a}_K\right]\right]$
\item $\vect{a}_L = h \:
\vect{\mathcal{N}}\left(\vect{x},\vect{a}_L,\vect{p},\vect{\xi}\right)$
\item $\vect{a}   = \vect{a} +c_{6}\vect{a}_J$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{6}^{*}\vect{a}_J$
\end{enumerate}
The $f_{i}$ coefficients for the last function evaluation are derived from the original Cash-Karp coefficients (see table \ref{tab:ck}) in the following way:
%
\begin{align}
g&=c_{1}c_{4}^{*}-c_{1}^{*}c_{4} \nonumber \\
f_{1}&=1 + (b_{64}(c_{1}^{*}-c_{1}) + b_{61}(c_{4}-c_{4}^{*}))/g  \nonumber \\
f_{2}&=b_{62}  \nonumber \\
f_{3}&=b_{63} + (b_{64}(c_{1}^{*}-c_{3}) + b_{61}( c_{3}^{*} c_{4} -  c_{3} c_{4}^{*} ))/g  \nonumber \\
f_{4}&=(b_{61}c_{4}^{*} -b_{64}c_{1}^{*})/g  \nonumber \\
f_{5}&=b_{65} +c_{5}^{*}(b_{61}c_{4} -b_{64}c_{1})/g  \nonumber \\
f_{6}&= (b_{64}c_{1} - b_{61}c_{4})/g \nonumber \\
\end{align}
%
This method requires six extra copies of the field. One, $\vect{a}_I$,
for the initial field value, three $\vect{a}_K$,  $\vect{a}_J$,  $\vect{a}_L$, to act as a
working fields while the original vector $\vect{a}$ and $\vect{a}^{*}$ collect the
derivative contributions. It also requires a similar size vector for the derivatives $\vect{p}$. 

\subsection{The Ninth Order Runge-Kutta Method in the Explicit Picture} 
\label{sec:rk9exMethod}

We continue to use the Fourier Transform (or Spectral) method to 
determine the transverse derivatives. The procedure is
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
% Step 1 
\item $x^0 = x^0 + a_{1}h$
\item $\vect{a}_a = \vect{a}$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_a\right]\right]$
\item $\vect{a}_a = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_a,\vect{p},\vect{\xi}\right)$
% Step 2 
\item $x^0 = x^0 + (a_{2}-a_{1})h$
\item $\vect{a}_b = \vect{a} + b_{2,1}\vect{a}_a$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_b\right]\right]$
\item $\vect{a}_b = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_b,\vect{p},\vect{\xi}\right)$
% Step 3 
\item $x^0 = x^0 + (a_{3}-a_{2})h$
\item $\vect{a}_c = \vect{a} + b_{3,1}\vect{a}_a + b_{3,2}\vect{a}_b$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_c\right]\right]$
\item $\vect{a}_c = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_c,\vect{p},\vect{\xi}\right)$
% Step 4 
\item $x^0 = x^0 + (a_{4}-a_3)h$
\item $\vect{a}_d = \vect{a} + b_{4,1}\vect{a}_a + b_{4,2}\vect{a}_b + b_{4,3}\vect{a}_c$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_d\right]\right]$
\item $\vect{a}_d = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_d,\vect{p},\vect{\xi}\right)$
% Step 5 
\item $x^0 = x^0 + (a_{5}-a_{4})h$
\item $\vect{a}_e = \vect{a} + b_{5,1}\vect{a}_a + b_{5,2}\vect{a}_b + b_{5,3}\vect{a}_c + b_{5,4}\vect{a}_d$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_e\right]\right]$
\item $\vect{a}_e = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_e,\vect{p},\vect{\xi}\right)$
% Step 6 
\item $x^0 = x^0 + (a_{6}-a_{5})h$
\item $\vect{a}_i = \vect{a} + b_{6,1}\vect{a}_a + b_{6,2}\vect{a}_b + b_{6,3}\vect{a}_c + b_{6,4}\vect{a}_d + b_{6,5}\vect{a}_e$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_f\right]\right]$
\item $\vect{a}_i = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_i,\vect{p},\vect{\xi}\right)$
% Step 7 
\item $x^0 = x^0 + (a_{7}-a_{6})h$
\item $\vect{a}_j = \vect{a} + b_{7,1}\vect{a}_a + b_{7,2}\vect{a}_b + b_{7,3}\vect{a}_c + b_{7,4}\vect{a}_d + b_{7,5}\vect{a}_e + b_{7,6}\vect{a}_i$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_j\right]\right]$
\item $\vect{a}_j = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_j,\vect{p},\vect{\xi}\right)$
% Step 8 
\item $x^0 = x^0 + (a_{8}-a_{7})h$
\item $\vect{a}_b = \vect{a} + b_{8,1}\vect{a}_a + b_{8,6}\vect{a}_i + b_{8,7}\vect{a}_j$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_b\right]\right]$
\item $\vect{a}_b = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_b,\vect{p},\vect{\xi}\right)$
% Step 9 
\item $x^0 = x^0 + (a_{9}-a_{8})h$
\item $\vect{a}_c = \vect{a} + b_{9,1}\vect{a}_a + b_{9,6}\vect{a}_i + b_{9,7}\vect{a}_j + b_{9,8}\vect{a}_b$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_c\right]\right]$
\item $\vect{a}_c = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_c,\vect{p},\vect{\xi}\right)$
% Step 10 
\item $x^0 = x^0 + (a_{10}-a_{9}))h$
\item $\vect{a}_d = \vect{a} + b_{10,1}\vect{a}_a + b_{10,6}\vect{a}_i + b_{10,7}\vect{a}_j + b_{10,8}\vect{a}_b + b_{10,9}\vect{a}_c$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_d\right]\right]$
\item $\vect{a}_d = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_d,\vect{p},\vect{\xi}\right)$
% Step 11 
\item $x^0 = x^0 + (a_{11}-a_{10})h$
\item $\vect{a}_e = \vect{a} + b_{11,1}\vect{a}_a + b_{11,6}\vect{a}_i + b_{11,7}\vect{a}_j + b_{11,8}\vect{a}_b + b_{11,9}\vect{a}_c + b_{11,10}\vect{a}_d$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_e\right]\right]$
\item $\vect{a}_e = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_e,\vect{p},\vect{\xi}\right)$
% Step 12 
\item $x^0 = x^0 + (a_{12}-a_{11})h$
\item $\vect{a}_f = \vect{a} + b_{12,1}\vect{a}_a + b_{12,6}\vect{a}_i + b_{12,7}\vect{a}_j + b_{12,8}\vect{a}_b + b_{12,9}\vect{a}_c + b_{12,10}\vect{a}_d + b_{12,11}\vect{a}_e$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_f\right]\right]$
\item $\vect{a}_f = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_f,\vect{p},\vect{\xi}\right)$
% Step 13 
\item $x^0 = x^0 + (a_{13}-a_{12})h$
\item $\vect{a}_g = \vect{a} + b_{13,1}\vect{a}_a + b_{13,6}\vect{a}_i + b_{13,7}\vect{a}_j + b_{13,8}\vect{a}_b + b_{13,9}\vect{a}_c + b_{13,10}\vect{a}_d + b_{13,11}\vect{a}_e + b_{13,12}\vect{a}_f$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_g\right]\right]$
\item $\vect{a}_g = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_g,\vect{p},\vect{\xi}\right)$
% Step 14 
\item $x^0 = x^0 + (a_{14}-a_{13})h$
\item $\vect{a}_h = \vect{a} + b_{14,1}\vect{a}_a + b_{14,6}\vect{a}_i + b_{14,7}\vect{a}_j + b_{14,8}\vect{a}_b + b_{14,9}\vect{a}_c + b_{14,10}\vect{a}_d + b_{14,11}\vect{a}_e + b_{14,12}\vect{a}_f + b_{14,13}\vect{a}_g$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_h\right]\right]$
\item $\vect{a}_h = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_h,\vect{p},\vect{\xi}\right)$
% Step 15 and 16 combined to reduce memory use 
\item $\vect{a}_i = \vect{a} + b_{15,1}\vect{a}_a + b_{15,6}\vect{a}_i + b_{15,7}\vect{a}_j + b_{15,8}\vect{a}_b + b_{15,9}\vect{a}_c + b_{15,10}\vect{a}_d + b_{15,11}\vect{a}_e + b_{15,12}\vect{a}_f + b_{15,13}\vect{a}_g + b_{15,14}\vect{a}_h$
\item $\vect{a}_j = (1-b_{16,6}/b_{15,6})\vect{a}  + (b_{16,1}-b_{15,1}b_{16,6}/b_{15,6})\vect{a}_a + (b_{16,6}/b_{15,6})\vect{a}_i + (b_{16,7}-b_{15,7}b_{16,6}/b_{15,6})\vect{a}_j + (b_{16,8}-b_{16,8}b_{16,6}/b_{15,6})\vect{a}_b + (b_{16,9}-b_{15,9}b_{16,6}/b_{15,6})\vect{a}_c + (b_{16,10}-b_{15,10}b_{16,6}/b_{15,6})\vect{a}_d + (b_{16,11}-b_{15,11}b_{16,6}/b_{15,6})\vect{a}_e + (b_{16,12}-b_{15,12}b_{16,6}/b_{15,6})\vect{a}_f + (b_{16,13}-b_{15,13}b_{16,6}/b_{15,6})\vect{a}_g + (b_{16,14}-b_{15,14}b_{16,6}/b_{15,6})\vect{a}_h$
\item $x^0 = x^0 + (a_{15}-a_{14})h$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_i\right]\right]$
\item $\vect{a}_i = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_i,\vect{p},\vect{\xi}\right)$
\item $x^0 = x^0 + (a_{16}-a_{15})h$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_j\right]\right]$
\item $\vect{a}_j = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_j,\vect{p},\vect{\xi}\right)$
% Take full step
\item $\vect{a} = \vect{a} + c_{1}\vect{a}_a + c_{8}\vect{a}_b + c_{9}\vect{a}_c + c_{10}\vect{a}_d + c_{11}\vect{a}_e + c_{12}\vect{a}_f + c_{13}\vect{a}_g + c_{14}\vect{a}_h + c_{15}\vect{a}_i + c_{16}\vect{a}_j$
\end{enumerate}
This method requires 10 extra copies of the field, $\vect{a}_a$ to $\vect{j}_j$. Many of the $a_i$ $b_{ij}$ and $c_i$ constants are zero in the  Eight/Ninth order Runge Kutta method are zero, meaning not all fields produced at each midstep are required all the time, this allows us to reduce the total number of fields required significantly. The 15th and 16th steps were merged together as well to reduce the total fields needed in memory. Stochastic problems when solved with this methods will achieve a significant performance boost many orders or magnitude greater than extra memory use.

\subsection{The  Eighth/Ninth Order adaptive Runge-Kutta Method in the Explicit Picture} 
\label{sec:rk89exMethod}

We continue to use the Fourier Transform (or Spectral) method to 
determine the transverse derivatives. The procedure is

\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
% Step 1 
\item $x^0 = x^0 + a_{1}h$
\item $\vect{a}_a = \vect{a}$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_a\right]\right]$
\item $\vect{a}_a = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_a,\vect{p},\vect{\xi}\right)$
% Step 2 
\item $x^0 = x^0 + (a_{2}-a_{1})h$
\item $\vect{a}_b = \vect{a} + b_{2,1}\vect{a}_a$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_b\right]\right]$
\item $\vect{a}_b = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_b,\vect{p},\vect{\xi}\right)$
% Step 3 
\item $x^0 = x^0 + (a_{3}-a_{2})h$
\item $\vect{a}_c = \vect{a} + b_{3,1}\vect{a}_a + b_{3,2}\vect{a}_b$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_c\right]\right]$
\item $\vect{a}_c = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_c,\vect{p},\vect{\xi}\right)$
% Step 4 
\item $x^0 = x^0 + (a_{4}-a_3)h$
\item $\vect{a}_d = \vect{a} + b_{4,1}\vect{a}_a + b_{4,2}\vect{a}_b + b_{4,3}\vect{a}_c$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_d\right]\right]$
\item $\vect{a}_d = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_d,\vect{p},\vect{\xi}\right)$
% Step 5 
\item $x^0 = x^0 + (a_{5}-a_{4})h$
\item $\vect{a}_e = \vect{a} + b_{5,1}\vect{a}_a + b_{5,2}\vect{a}_b + b_{5,3}\vect{a}_c + b_{5,4}\vect{a}_d$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_e\right]\right]$
\item $\vect{a}_e = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_e,\vect{p},\vect{\xi}\right)$
% Step 6 
\item $x^0 = x^0 + (a_{6}-a_{5})h$
\item $\vect{a}_i = \vect{a} + b_{6,1}\vect{a}_a + b_{6,2}\vect{a}_b + b_{6,3}\vect{a}_c + b_{6,4}\vect{a}_d + b_{6,5}\vect{a}_e$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_f\right]\right]$
\item $\vect{a}_i = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_i,\vect{p},\vect{\xi}\right)$
% Step 7 
\item $x^0 = x^0 + (a_{7}-a_{6})h$
\item $\vect{a}_j = \vect{a} + b_{7,1}\vect{a}_a + b_{7,2}\vect{a}_b + b_{7,3}\vect{a}_c + b_{7,4}\vect{a}_d + b_{7,5}\vect{a}_e + b_{7,6}\vect{a}_i$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_j\right]\right]$
\item $\vect{a}_j = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_j,\vect{p},\vect{\xi}\right)$
% Step 8 
\item $x^0 = x^0 + (a_{8}-a_{7})h$
\item $\vect{a}_b = \vect{a} + b_{8,1}\vect{a}_a + b_{8,6}\vect{a}_i + b_{8,7}\vect{a}_j$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_b\right]\right]$
\item $\vect{a}_b = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_b,\vect{p},\vect{\xi}\right)$
% Step 9 
\item $x^0 = x^0 + (a_{9}-a_{8})h$
\item $\vect{a}_c = \vect{a} + b_{9,1}\vect{a}_a + b_{9,6}\vect{a}_i + b_{9,7}\vect{a}_j + b_{9,8}\vect{a}_b$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_c\right]\right]$
\item $\vect{a}_c = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_c,\vect{p},\vect{\xi}\right)$
% Step 10 
\item $x^0 = x^0 + (a_{10}-a_{9}))h$
\item $\vect{a}_d = \vect{a} + b_{10,1}\vect{a}_a + b_{10,6}\vect{a}_i + b_{10,7}\vect{a}_j + b_{10,8}\vect{a}_b + b_{10,9}\vect{a}_c$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_d\right]\right]$
\item $\vect{a}_d = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_d,\vect{p},\vect{\xi}\right)$
% Step 11 
\item $x^0 = x^0 + (a_{11}-a_{10})h$
\item $\vect{a}_e = \vect{a} + b_{11,1}\vect{a}_a + b_{11,6}\vect{a}_i + b_{11,7}\vect{a}_j + b_{11,8}\vect{a}_b + b_{11,9}\vect{a}_c + b_{11,10}\vect{a}_d$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_e\right]\right]$
\item $\vect{a}_e = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_e,\vect{p},\vect{\xi}\right)$
% Step 12 
\item $x^0 = x^0 + (a_{12}-a_{11})h$
\item $\vect{a}_f = \vect{a} + b_{12,1}\vect{a}_a + b_{12,6}\vect{a}_i + b_{12,7}\vect{a}_j + b_{12,8}\vect{a}_b + b_{12,9}\vect{a}_c + b_{12,10}\vect{a}_d + b_{12,11}\vect{a}_e$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_f\right]\right]$
\item $\vect{a}_f = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_f,\vect{p},\vect{\xi}\right)$
% Step 13 
\item $x^0 = x^0 + (a_{13}-a_{12})h$
\item $\vect{a}_g = \vect{a} + b_{13,1}\vect{a}_a + b_{13,6}\vect{a}_i + b_{13,7}\vect{a}_j + b_{13,8}\vect{a}_b + b_{13,9}\vect{a}_c + b_{13,10}\vect{a}_d + b_{13,11}\vect{a}_e + b_{13,12}\vect{a}_f$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_g\right]\right]$
\item $\vect{a}_g = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_g,\vect{p},\vect{\xi}\right)$
% Step 14 
\item $x^0 = x^0 + (a_{14}-a_{13})h$
\item $\vect{a}_h = \vect{a} + b_{14,1}\vect{a}_a + b_{14,6}\vect{a}_i + b_{14,7}\vect{a}_j + b_{14,8}\vect{a}_b + b_{14,9}\vect{a}_c + b_{14,10}\vect{a}_d + b_{14,11}\vect{a}_e + b_{14,12}\vect{a}_f + b_{14,13}\vect{a}_g$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_h\right]\right]$
\item $\vect{a}_h = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_h,\vect{p},\vect{\xi}\right)$
% Step 15 and 16 combined to reduce memory use 
\item $\vect{a}_i = \vect{a} + b_{15,1}\vect{a}_a + b_{15,6}\vect{a}_i + b_{15,7}\vect{a}_j + b_{15,8}\vect{a}_b + b_{15,9}\vect{a}_c + b_{15,10}\vect{a}_d + b_{15,11}\vect{a}_e + b_{15,12}\vect{a}_f + b_{15,13}\vect{a}_g + b_{15,14}\vect{a}_h$
\item $\vect{a}_j = (1-b_{16,6}/b_{15,6})\vect{a}  + (b_{16,1}-b_{15,1}b_{16,6}/b_{15,6})\vect{a}_a + (b_{16,6}/b_{15,6})\vect{a}_i + (b_{16,7}-b_{15,7}b_{16,6}/b_{15,6})\vect{a}_j + (b_{16,8}-b_{16,8}b_{16,6}/b_{15,6})\vect{a}_b + (b_{16,9}-b_{15,9}b_{16,6}/b_{15,6})\vect{a}_c + (b_{16,10}-b_{15,10}b_{16,6}/b_{15,6})\vect{a}_d + (b_{16,11}-b_{15,11}b_{16,6}/b_{15,6})\vect{a}_e + (b_{16,12}-b_{15,12}b_{16,6}/b_{15,6})\vect{a}_f + (b_{16,13}-b_{15,13}b_{16,6}/b_{15,6})\vect{a}_g + (b_{16,14}-b_{15,14}b_{16,6}/b_{15,6})\vect{a}_h$
\item $x^0 = x^0 + (a_{15}-a_{14})h$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_i\right]\right]$
\item $\vect{a}_i = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_i,\vect{p},\vect{\xi}\right)$
\item $x^0 = x^0 + (a_{16}-a_{15})h$
\item $\vect{p} = \mathcal{F}^{-1}\left[\mathcal{L}\left(x^0,\vect{k}_{\bot}\right) \cdot \mathcal{F}\left[\vect{a}_j\right]\right]$
\item $\vect{a}_j = \vect{\mathcal{N}}\left(\vect{x},\vect{a}_j,\vect{p},\vect{\xi}\right)$
% Take full step
\item $\vect{a}_{Int} = \vect{a}$
\item $\vect{a} = \vect{a} + c_{1}\vect{a}_a + c_{8}\vect{a}_b + c_{9}\vect{a}_c + c_{10}\vect{a}_d + c_{11}\vect{a}_e + c_{12}\vect{a}_f + c_{13}\vect{a}_g + c_{14}\vect{a}_h + c_{15}\vect{a}_i + c_{16}\vect{a}_j$
\item $\vect{a}_a = \vect{a}_{Int} + c_{1}^*\vect{a}_a + c_{8}^*\vect{a}_b + c_{9}^*\vect{a}_c + c_{10}^*\vect{a}_d + c_{11}^*\vect{a}_e + c_{12}^*\vect{a}_f + c_{13}^*\vect{a}_g + c_{14}^*\vect{a}_h + c_{15}^*\vect{a}_i + c_{16}^*\vect{a}_j$
\end{enumerate}

This method requires 11 extra copies of the field. An additional copy is required for the adaptive time-step (compared to the non adaptive ARK9EX) since the approximate error in the step must be calculated. This is achieved using the field stored in $\vect{a}_a$ which holds the 8th order solution, analogous to the $\vect{a}^*$ in the ARK45 method. Once again we should note the adaptive step implemented in this method is stochastically safe, making it the best performing solver for stochastic problems.

\subsection{Interaction Picture Methods} 
\label{sec:ipMethod}

Interaction picture methods can only be applied to equations of the
general form:
\begin{equation}
\frac{\partial}{\partial x^0}a^i(x^0,\vect{x}_\bot) =
 \mathcal{F}^{-1}\left[\mathcal{L}^i(x^0,\vect{k}_{\bot})
 \mathcal{F}[a^i(\vect{x})] \right] +
 \mathcal{N}^i\left(\vect{x}, \vect{a}(\vect{x}),
 \vect{\xi}(\vect{x})\right).
 \label{eq:ipPde}
\end{equation}
Although it is possible to define interaction picture algorithms in
which the linear operators may have the general matrix form
$\mathcal{L}^{ij}$, here they are restricted to the diagonal form
$\mathcal{L}^i$. Non-diagonal operators require eigen-vector-value
decomposition each time step and at each lattice point, which becomes
computationally expensive. These $\mathcal{L}^i$ operators may still
have coefficients that depend on the the propagation dimension and an
the Fourier space coordinates $\vect{k}_\bot$.

The evolution of the field $\vect{a}$ is carried out in both normal
space {\em and} Fourier space in alternating steps, hence this method
is also known as a {\em Split-Step} method, or a {\em Split-Operator}
method. The exact implementation depends on the main integration
method employed, as is detailed in the following Sections,
\ref{sec:siipMethod} and \ref{sec:rk4ipMethod}.

As a general rule the interaction picture algorithms are faster and
more stable than their Explicit picture counterparts, but they are
also more restrictive with regard to the allowable linear operators.

\subsection{The Semi-Implicit method in the Interaction Picture} 
\label{sec:siipMethod}

The original derivation of this method was performed by
Drummond~\cite{Drummond:1983:1}, and we repeat it here since it is the
most efficient and stable method for solving stochastic PDEs. Consider
now expressing the field $a$ as the transform:
\begin{equation}
a^i(x^0,\vect{x}_\bot) =
e^{\left(x^0-z\right)\mathcal{L}^i\left(x^0,\vect{k}_{\bot}\right)}
\left[b^i(x^0,\vect{x}_\bot)\right].
\end{equation}

The time derivative of $\vect{a}$ then becomes:
\begin{equation}
\frac{\partial}{\partial x^0} a^i(x^0,\vect{x}_\bot) =
  \mathcal{L}^i\left(x^0,\vect{k}_{\bot}\right) \left[a^i\right] +
  e^{\left(x^0-z\right)f^i(x^0,\vect{k}_\bot)\mathcal{L}^i}
  \left[\dot{b}^i\right],
\end{equation}
since the operators $\mathcal{L}$ are linear. Equating this result
with \eqn{eq:ipPde} we get:
\begin{equation}
\dot{b}^i(x^0,\vect{x}_\bot) =
   e^{-\left(x^0-z\right)\mathcal{L}^i\left(x^0,\vect{k}_{\bot}\right)}
   \left[\mathcal{N}^i\left(\vect{x},\vect{a},\vect{\xi}\right)\right].
\end{equation}

Now solving for $\vect{b}$ using the semi-implicit method yields:
\begin{align}
  b^i\left(x^0 + \frac{h}{2},\vect{x}_\bot\right) & =
  b^i\left(x^0,\vect{x}_\bot\right)\\ 
  & + \frac{1}{2} \: h \:
  e^{-\left(x^0+\frac{h}{2}-z\right)\mathcal{L}^i\left(x^0,\vect{k}_{\bot}\right)}
  \left[ \mathcal{N}^i\left(x^0 + \frac{h}{2}, \vect{x}_\bot,
    \vect{a}\left(x^0 + \frac{h}{2},\vect{x}_\bot\right),\vect{\xi}\right)
    \right].
\end{align}

Finally, if we define $z=x^0+\frac{h}{2}$ then $\vect{a}$ and $\vect{b}$
become identical at the middle of the time step, and we obtain
\begin{equation}
  b^i\left(x^0 + \frac{h}{2},\vect{x}_\bot\right) =
  b^i\left(x^0,\vect{x}_\bot\right) + \frac{h}{2}
  \mathcal{N}^i\left(x^0 + \frac{h}{2}, \vect{x}_\bot,
  \vect{b}\left(x^0 +
  \frac{h}{2},\vect{x}_\bot\right),\vect{\xi}\right).
\label{eq:siip1}
\end{equation}
Hence the $\vect{a}$ and $\vect{b}$ fields are transformed from one to
the other at the beginning and end of the time step using the linear
operators, and the $\vect{b}$ field is evolved over the time step using
the nonlinear operators and the semi-implicit algorithm.  Thus the
algorithm becomes quite simply:
\begin{enumerate}
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item Calculate $\vect{\xi}$ if required.
\item $x^0 = x^0 + \frac{1}{2}h$
\item For each point in space do:
  \begin{enumerate}
  \item $\vect{a_I}=\vect{a}$
  \item For N-1 iterations do:
    \begin{enumerate}
    \item $\vect{a} = \vect{a}_{I} + \frac{1}{2} h \:
      \vect{\mathcal{N}}(\vect{x},\vect{a},\vect{\xi})$
    \end{enumerate}
  \item $\vect{a} = \vect{a}_{I} + h \:
  \vect{\mathcal{N}}(\vect{x},\vect{a},\vect{\xi})$
  \end{enumerate}
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\end{enumerate}

The first and the last steps are easily performed with the field in
Fourier space since the linear operators reduce to $k_\bot$
multipliers. Note that the dependence of the linear operators
$\mathcal{L}^i(x^0,\vect{k}_\bot)$ on the propagation dimension ought
to be weak in comparison to the size of the time step. If they do vary
with time then the exponentials will have to be calculated at every
time step, which can be computationally expensive with some CPU
architectures. Otherwise the exponentials may be pre-calculated and
tabulated for reference.

Further advantages of this method are that the forward and backward
Fourier transforms may be performed {\em in place} (i.e. within the
current memory block), as can the point-by-point multiplication of the
field by the exponentials. Thus no extra copies of the field are
required.

\subsection{The Fourth Order Runge-Kutta Method in the Interaction Picture} 
\label{sec:rk4ipMethod}

One way of reducing the memory overhead on the RK4 algorithm is to
move into an interaction picture, exactly as was performed in the
split-step semi-implicit method above.  This method was derived by Rob
Ballagh's BEC group at the University of Otago.  It is described in
detail in the PhD thesis of
B.M. Caradoc-Davies~\cite{web:Caradoc-Davies}.

We will present only an optimised recipe for implementing this
algorithm:
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
\item $\vect{a}_K =  \vect{a}$
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item $\vect{a}_I =  \vect{a}$
\item $\vect{a}_K = h \: \vect{\mathcal{N}}(\vect{x},\vect{a}_K,\vect{\xi})$
\item $\vect{a}_K =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_K$
\item $\vect{a} = \vect{a}+\frac{1}{6}\vect{a}_K$
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a}_K = \vect{a}_I + \frac{1}{2}\vect{a}_K$
\item $\vect{a}_K = h \: \vect{\mathcal{N}}(\vect{x},\vect{a}_K,\vect{\xi})$
\item $\vect{a} = \vect{a}+\frac{1}{3}\vect{a}_K$
\item $\vect{a}_K = \vect{a}_I + \frac{1}{2}\vect{a}_K$
\item $\vect{a}_K = h \: \vect{\mathcal{N}}(\vect{x},\vect{a}_K,\vect{\xi})$
\item $\vect{a} = \vect{a}+\frac{1}{3}\vect{a}_K$
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a}_K = \vect{a}_I + \vect{a}_K$
\item $\vect{a}_K =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_K$
\item $\vect{a}_K = h \:
\vect{\mathcal{N}}(\vect{x},\vect{a}_K,\vect{\xi})$
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item $\vect{a} = \vect{a}+\frac{1}{6}\vect{a}_K$
\end{enumerate}

Similarly to the semi-implicit algorithm in the interaction picture,
the first and the last steps are easily performed with the field in
Fourier space since the linear operators reduce to $k_\bot$
multipliers. And again note that the dependence of the linear
operators $\mathcal{L}^i(x^0,\vect{k}_\bot)$ on the propagation
dimension ought to be weak in comparison to the size of the time step.

This algorithm has a similar memory overhead as its Explicit picture
counterpart, but it no longer needs the memory for the $\vect{p}$
derivatives vector, nor does it have to do the work to calculate this
vector.

The interaction picture method relies on the derivative Operator being independent of the propagation dimension, otherwise it cannot function as a fourth order algorithm.

\subsection{The Fourth/ Fifth Order adaptive Runge-Kutta Method in the Interaction Picture} 
\label{sec:ark45ipMethod}

In the Fourth/ Fifth Order adaptive Runge-Kutta algorithm it seems we would require too many Fourier transforms for this to be efficient. Following \cite{web:MatthewDavis} the function can however primarily be evolved in Fourier space and transformed into normal space for calculation of the $\vect{\mathcal{N}}(\vect{x},\vect{a}_i,\vect{\xi})$ only. The use of the interaction picture does not allow a reduction of the computational effort to the same extent as in the RK4IP algorithm, but this method can still be vastly superior over the adaptive explicit picture method as it allows larger step sizes for equations containing certain derivative operators. In the following the vectors $\vect{a}$, $\vect{a}^{*}$ are supposed to be initially in Fourier space. Note that for the calculation of the $\vect{\mathcal{N}}()$ these are transformed into normal space. 

\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
\item $\vect{a}_K =\vect{a}$
\item $\vect{a} =
e^{a_{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item $\vect{a}_{I} =\vect{a}$
\item $\vect{a}^{*} = \vect{a}$
\item $\vect{a}_K = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0, \mathcal{F}^{-1}\left[\vect{a}_K\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{K} =
e^{-a_{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{K}$
\item $\vect{a}   = \vect{a} +c_{1}\vect{a}_K$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{1}^{*}\vect{a}_K$
\item $x^0 = x^0 + a_{2}h$
\item $\vect{a}_K = \vect{a}_I + b_{21}\vect{a}_K$
\item $\vect{a}_K = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_K\right],\vect{\xi}\right)\right]$
\item $\vect{a}_J =(1-b_{31}/c_{1}) \vect{a}_I + b_{31}/c_{1}\vect{a} +b_{32} \vect{a}_K$
\item $x^0 = x^0 + (a_{3}-a_{2})h$ \\
	(Note: $c_{2}=c_{2}^{*}=0$) 
\item $\vect{a}_{J} =
e^{(a_{3}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{J}$	
\item $\vect{a}_J = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_J\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{J} =
e^{-(a_{3}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{J}$	
\item $\vect{a}_L =(1-b_{41}/c_{1}) \vect{a}_I + b_{41}/c_{1}\vect{a} +b_{42} \vect{a}_K +b_{43}\vect{a}_{J}$
\item $\vect{a}   = \vect{a} +c_{3}\vect{a}_J$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{3}^{*}\vect{a}_J$
\item $x^0 = x^0 + (a_{4}-a_{3})h$ \\
\item $\vect{a}_{L} =
e^{(a_{4}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$	
\item $\vect{a}_L = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(\vect{x},\mathcal{F}^{-1}\left[\vect{a}_L\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{L} =
e^{-(a_{4}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$	
\item $\vect{a}   = \vect{a} +c_{4}\vect{a}_L$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{4}^{*}\vect{a}_L$
\item $\vect{a}_L =(1-b_{51}/c_{1}) \vect{a}_I + b_{51}/c_{1}\vect{a} +b_{52} \vect{a}_K +(b_{53}-b_{51}c_{3}/c_{1})\vect{a}_{J} +(b_{54}-b_{51}c_{4}/c_{1})\vect{a}_{L} $
\item $x^0 = x^0 + (a_{5}-a_{4})h$ \\
\item $\vect{a}_{L} =
e^{(a_{5}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$	
\item $\vect{a}_L = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(\vect{x},\mathcal{F}^{-1}\left[\vect{a}_L\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{L} =
e^{-(a_{5}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$	
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{5}^{*}\vect{a}_L$\\
	(Note: $c_{5}=0$) 
\item $\vect{a}_L =f_{1} \vect{a}_I  +f_{2} \vect{a}_K + f_{3}\vect{a}_{J} + f_{4} \vect{a} +f_{5}\vect{a}_{L} + f_{6} \vect{a}^{*}$		
\item $x^0 = x^0 + (a_{6}-a_{5})h$ \\
\item $\vect{a}_{L} =
e^{(a_{6}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$
\item $\vect{a}_L = h \:
\mathcal{F}\left[\vect{\mathcal{N}}\left(\vect{x},\mathcal{F}^{-1}\left[\vect{a}_L\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{L} =
e^{-(a_{6}-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}_{L}$
\item $\vect{a}   = \vect{a} +c_{6}\vect{a}_J$
\item $\vect{a}^{*}   = \vect{a}^{*} +c_{6}^{*}\vect{a}_J$
\item $x^0 = x^0 + (a_{6}-a_{5})h$ \\
\item $\vect{a} =
e^{(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item $\vect{a}^{*}=
e^{(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}^{*}$
\end{enumerate}
The $f_{i}$ coefficients are listed in section \ref{sec:ark45exMethod}.
Unlike its fixed step counterpart (RK4IP), this algorithm requires significant memory for the fastest calculation of the Fourier space propagation, as the arguments of the exponentials contain different factors of $a_{i}-a_{2}$ for each application. If this turns out to be a problem, XMDS offers the possibility to sacrifice speed for less memory consumption by the use of an appropriate flag. 

As the two computed solutions are of fourth and fifth order \emph{only} if the derivative operators are independent of the propagation dimension, this algorithm cannot be used at all for problems where this is not the case. 


\subsection{The  Ninth Order Runge-Kutta Method in the Interaction Picture} 
\label{sec:rk9ipMethod}

The ninth order Runge-Kutta does require a large amount of Fourier transforms, since there is no temporal Symmetry in the algorithm. However the improvement in convergence may offset the extra time required, particularly in stochastic problems.
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
% Step 1 
\item $x^0 = x^0 + a_{1}h$
\item $\vect{a}_a = \vect{a}$
\item $\vect{a} = e^{h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}$
\item $\vect{a}_a = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_a\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{a} = e^{h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{a}$
% Step 2 
\item $x^0 = x^0 + (a_{2}-a_{1})h$
\item $\vect{a}_b = \vect{a} + b_{2,1}\vect{a}_a$
\item $\vect{a}_{b} = e^{-(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
\item $\vect{a}_b = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_b\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{b} = e^{(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
% Step 3 
\item $x^0 = x^0 + (a_{3}-a_{2})h$
\item $\vect{a}_c = \vect{a} + b_{3,1}\vect{a}_a + b_{3,2}\vect{a}_b$
\item $\vect{a}_{c} = e^{-(1-a_{3})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
\item $\vect{a}_c = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_c\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{c} = e^{(1-a_{3})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
% Step 4 
\item $x^0 = x^0 + (a_{4}-a_3)h$
\item $\vect{a}_d = \vect{a} + b_{4,1}\vect{a}_a + b_{4,2}\vect{a}_b + b_{4,3}\vect{a}_c$
\item $\vect{a}_{d} = e^{-(1-a_{4})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
\item $\vect{a}_d = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_d\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{d} = e^{(1-a_{4})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
% Step 5 
\item $x^0 = x^0 + (a_{5}-a_{4})h$
\item $\vect{a}_e = \vect{a} + b_{5,1}\vect{a}_a + b_{5,2}\vect{a}_b + b_{5,3}\vect{a}_c + b_{5,4}\vect{a}_d$
\item $\vect{a}_{e} = e^{-(1-a_{5})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
\item $\vect{a}_e = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_e\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{e} = e^{(1-a_{5})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
% Step 6 
\item $x^0 = x^0 + (a_{6}-a_{5})h$
\item $\vect{a}_i = \vect{a} + b_{6,1}\vect{a}_a + b_{6,2}\vect{a}_b + b_{6,3}\vect{a}_c + b_{6,4}\vect{a}_d + b_{6,5}\vect{a}_e$
\item $\vect{a}_{i} = e^{-(1-a_{6})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{i}$
\item $\vect{a}_i = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_i\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{i} = e^{(1-a_{6})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{i}$
% Step 7 
\item $x^0 = x^0 + (a_{7}-a_{6})h$
\item $\vect{a}_j = \vect{a} + b_{7,1}\vect{a}_a + b_{7,2}\vect{a}_b + b_{7,3}\vect{a}_c + b_{7,4}\vect{a}_d + b_{7,5}\vect{a}_e + b_{7,6}\vect{a}_i$
\item $\vect{a}_{j} = e^{-(1-a_{7})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{j}$
\item $\vect{a}_j = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_j\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{j} = e^{(1-a_{7})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{j}$
% Step 8 
\item $x^0 = x^0 + (a_{8}-a_{7})h$
\item $\vect{a}_b = \vect{a} + b_{8,1}\vect{a}_a + b_{8,6}\vect{a}_i + b_{8,7}\vect{a}_j$
\item $\vect{a}_{b} = e^{-(1-a_{8})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
\item $\vect{a}_b = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_b\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{b} = e^{(1-a_{8})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
% Step 9 
\item $x^0 = x^0 + (a_{9}-a_{8})h$
\item $\vect{a}_c = \vect{a} + b_{9,1}\vect{a}_a + b_{9,6}\vect{a}_i + b_{9,7}\vect{a}_j + b_{9,8}\vect{a}_b$
\item $\vect{a}_{c} = e^{-(1-a_{9})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
\item $\vect{a}_c = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_c\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{c} = e^{(1-a_{9})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
% Step 10 
\item $x^0 = x^0 + (a_{10}-a_{9}))h$
\item $\vect{a}_d = \vect{a} + b_{10,1}\vect{a}_a + b_{10,6}\vect{a}_i + b_{10,7}\vect{a}_j + b_{10,8}\vect{a}_b + b_{10,9}\vect{a}_c$
\item $\vect{a}_{d} = e^{-(1-a_{10})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
\item $\vect{a}_d = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_d\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{d} = e^{(1-a_{10})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
% Step 11 
\item $x^0 = x^0 + (a_{11}-a_{10})h$
\item $\vect{a}_e = \vect{a} + b_{11,1}\vect{a}_a + b_{11,6}\vect{a}_i + b_{11,7}\vect{a}_j + b_{11,8}\vect{a}_b + b_{11,9}\vect{a}_c + b_{11,10}\vect{a}_d$
\item $\vect{a}_{e} = e^{-(1-a_{11})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
\item $\vect{a}_e = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_e\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{e} = e^{(1-a_{11})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
% Step 12 
\item $x^0 = x^0 + (a_{12}-a_{11})h$
\item $\vect{a}_f = \vect{a} + b_{12,1}\vect{a}_a + b_{12,6}\vect{a}_i + b_{12,7}\vect{a}_j + b_{12,8}\vect{a}_b + b_{12,9}\vect{a}_c + b_{12,10}\vect{a}_d + b_{12,11}\vect{a}_e$
\item $\vect{a}_{f} = e^{-(1-a_{12})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{f}$
\item $\vect{a}_f = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_f\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{f} = e^{(1-a_{12})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{f}$
% Step 13 
\item $x^0 = x^0 + (a_{13}-a_{12})h$
\item $\vect{a}_g = \vect{a} + b_{13,1}\vect{a}_a + b_{13,6}\vect{a}_i + b_{13,7}\vect{a}_j + b_{13,8}\vect{a}_b + b_{13,9}\vect{a}_c + b_{13,10}\vect{a}_d + b_{13,11}\vect{a}_e + b_{13,12}\vect{a}_f$
\item $\vect{a}_{g} = e^{-(1-a_{13})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{g}$
\item $\vect{a}_g = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_g\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{g} = e^{(1-a_{13})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{g}$
% Step 14 
\item $x^0 = x^0 + (a_{14}-a_{13})h$
\item $\vect{a}_h = \vect{a} + b_{14,1}\vect{a}_a + b_{14,6}\vect{a}_i + b_{14,7}\vect{a}_j + b_{14,8}\vect{a}_b + b_{14,9}\vect{a}_c + b_{14,10}\vect{a}_d + b_{14,11}\vect{a}_e + b_{14,12}\vect{a}_f + b_{14,13}\vect{a}_g$
\item $\vect{a}_{h} = e^{-(1-a_{14})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{h}$
\item $\vect{a}_h = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_h\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{h} = e^{(1-a_{14})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{h}$
% Step 15 and 16 combined to reduce memory use 
\item $\vect{a}_i = \vect{a} + b_{15,1}\vect{a}_a + b_{15,6}\vect{a}_i + b_{15,7}\vect{a}_j + b_{15,8}\vect{a}_b + b_{15,9}\vect{a}_c + b_{15,10}\vect{a}_d + b_{15,11}\vect{a}_e + b_{15,12}\vect{a}_f + b_{15,13}\vect{a}_g + b_{15,14}\vect{a}_h$
\item $\vect{a}_j = (1-b_{16,6}/b_{15,6})\vect{a}  + (b_{16,1}-b_{15,1}b_{16,6}/b_{15,6})\vect{a}_a + (b_{16,6}/b_{15,6})\vect{a}_i + (b_{16,7}-b_{15,7}b_{16,6}/b_{15,6})\vect{a}_j + (b_{16,8}-b_{16,8}b_{16,6}/b_{15,6})\vect{a}_b + (b_{16,9}-b_{15,9}b_{16,6}/b_{15,6})\vect{a}_c + (b_{16,10}-b_{15,10}b_{16,6}/b_{15,6})\vect{a}_d + (b_{16,11}-b_{15,11}b_{16,6}/b_{15,6})\vect{a}_e + (b_{16,12}-b_{15,12}b_{16,6}/b_{15,6})\vect{a}_f + (b_{16,13}-b_{15,13}b_{16,6}/b_{15,6})\vect{a}_g + (b_{16,14}-b_{15,14}b_{16,6}/b_{15,6})\vect{a}_h$
\item $x^0 = x^0 + (a_{15}-a_{14})h$
\item $\vect{a}_i = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_i\right],\vect{\xi}\right)\right]$
\item $x^0 = x^0 + (a_{16}-a_{15})h$
\item $\vect{a}_j = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_j\right],\vect{\xi}\right)\right]$
% Take full step
\item $\vect{a} = \vect{a} + c_{1}\vect{a}_a + c_{8}\vect{a}_b + c_{9}\vect{a}_c + c_{10}\vect{a}_d + c_{11}\vect{a}_e + c_{12}\vect{a}_f + c_{13}\vect{a}_g + c_{14}\vect{a}_h + c_{15}\vect{a}_i + c_{16}\vect{a}_j$
\end{enumerate}
The method requires 10 copies of the field like its explicit picture counterpart. Note no rotations are required for steps 15 and 16 since they both occur at the end of the time step. 

\subsection{The  Eighth / Ninth Order adaptive Runge-Kutta Method in the Interaction Picture} 
\label{sec:rk89ipMethod}

The Eighth/Ninth order Runge-Kutta also requires a large number of Fourier transforms since there is no temporal symmetry in the algorithm. However, the improvement in convergence may offset the extra time required, particularly in stochastic problems.
	
\begin{enumerate}
\item Calculate $\vect{\xi}$ if required.
% Step 1 
\item $x^0 = x^0 + a_{1}h$
\item $\vect{a}_a = \vect{a}$
\item $\vect{a} = e^{h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}$
\item $\vect{a}_a = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_a\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{a} = e^{h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{a}$
% Step 2 
\item $x^0 = x^0 + (a_{2}-a_{1})h$
\item $\vect{a}_b = \vect{a} + b_{2,1}\vect{a}_a$
\item $\vect{a}_{b} = e^{-(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
\item $\vect{a}_b = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_b\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{b} = e^{(1-a_{2})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
% Step 3 
\item $x^0 = x^0 + (a_{3}-a_{2})h$
\item $\vect{a}_c = \vect{a} + b_{3,1}\vect{a}_a + b_{3,2}\vect{a}_b$
\item $\vect{a}_{c} = e^{-(1-a_{3})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
\item $\vect{a}_c = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_c\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{c} = e^{(1-a_{3})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
% Step 4 
\item $x^0 = x^0 + (a_{4}-a_3)h$
\item $\vect{a}_d = \vect{a} + b_{4,1}\vect{a}_a + b_{4,2}\vect{a}_b + b_{4,3}\vect{a}_c$
\item $\vect{a}_{d} = e^{-(1-a_{4})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
\item $\vect{a}_d = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_d\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{d} = e^{(1-a_{4})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
% Step 5 
\item $x^0 = x^0 + (a_{5}-a_{4})h$
\item $\vect{a}_e = \vect{a} + b_{5,1}\vect{a}_a + b_{5,2}\vect{a}_b + b_{5,3}\vect{a}_c + b_{5,4}\vect{a}_d$
\item $\vect{a}_{e} = e^{-(1-a_{5})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
\item $\vect{a}_e = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_e\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{e} = e^{(1-a_{5})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
% Step 6 
\item $x^0 = x^0 + (a_{6}-a_{5})h$
\item $\vect{a}_i = \vect{a} + b_{6,1}\vect{a}_a + b_{6,2}\vect{a}_b + b_{6,3}\vect{a}_c + b_{6,4}\vect{a}_d + b_{6,5}\vect{a}_e$
\item $\vect{a}_{i} = e^{-(1-a_{6})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{i}$
\item $\vect{a}_i = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_i\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{i} = e^{(1-a_{6})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{i}$
% Step 7 
\item $x^0 = x^0 + (a_{7}-a_{6})h$
\item $\vect{a}_j = \vect{a} + b_{7,1}\vect{a}_a + b_{7,2}\vect{a}_b + b_{7,3}\vect{a}_c + b_{7,4}\vect{a}_d + b_{7,5}\vect{a}_e + b_{7,6}\vect{a}_i$
\item $\vect{a}_{j} = e^{-(1-a_{7})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{j}$
\item $\vect{a}_j = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_j\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{j} = e^{(1-a_{7})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{j}$
% Step 8 
\item $x^0 = x^0 + (a_{8}-a_{7})h$
\item $\vect{a}_b = \vect{a} + b_{8,1}\vect{a}_a + b_{8,6}\vect{a}_i + b_{8,7}\vect{a}_j$
\item $\vect{a}_{b} = e^{-(1-a_{8})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
\item $\vect{a}_b = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_b\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{b} = e^{(1-a_{8})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{b}$
% Step 9 
\item $x^0 = x^0 + (a_{9}-a_{8})h$
\item $\vect{a}_c = \vect{a} + b_{9,1}\vect{a}_a + b_{9,6}\vect{a}_i + b_{9,7}\vect{a}_j + b_{9,8}\vect{a}_b$
\item $\vect{a}_{c} = e^{-(1-a_{9})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
\item $\vect{a}_c = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_c\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{c} = e^{(1-a_{9})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{c}$
% Step 10 
\item $x^0 = x^0 + (a_{10}-a_{9}))h$
\item $\vect{a}_d = \vect{a} + b_{10,1}\vect{a}_a + b_{10,6}\vect{a}_i + b_{10,7}\vect{a}_j + b_{10,8}\vect{a}_b + b_{10,9}\vect{a}_c$
\item $\vect{a}_{d} = e^{-(1-a_{10})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
\item $\vect{a}_d = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_d\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{d} = e^{(1-a_{10})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{d}$
% Step 11 
\item $x^0 = x^0 + (a_{11}-a_{10})h$
\item $\vect{a}_e = \vect{a} + b_{11,1}\vect{a}_a + b_{11,6}\vect{a}_i + b_{11,7}\vect{a}_j + b_{11,8}\vect{a}_b + b_{11,9}\vect{a}_c + b_{11,10}\vect{a}_d$
\item $\vect{a}_{e} = e^{-(1-a_{11})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
\item $\vect{a}_e = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_e\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{e} = e^{(1-a_{11})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{e}$
% Step 12 
\item $x^0 = x^0 + (a_{12}-a_{11})h$
\item $\vect{a}_f = \vect{a} + b_{12,1}\vect{a}_a + b_{12,6}\vect{a}_i + b_{12,7}\vect{a}_j + b_{12,8}\vect{a}_b + b_{12,9}\vect{a}_c + b_{12,10}\vect{a}_d + b_{12,11}\vect{a}_e$
\item $\vect{a}_{f} = e^{-(1-a_{12})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{f}$
\item $\vect{a}_f = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_f\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{f} = e^{(1-a_{12})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{f}$
% Step 13 
\item $x^0 = x^0 + (a_{13}-a_{12})h$
\item $\vect{a}_g = \vect{a} + b_{13,1}\vect{a}_a + b_{13,6}\vect{a}_i + b_{13,7}\vect{a}_j + b_{13,8}\vect{a}_b + b_{13,9}\vect{a}_c + b_{13,10}\vect{a}_d + b_{13,11}\vect{a}_e + b_{13,12}\vect{a}_f$
\item $\vect{a}_{g} = e^{-(1-a_{13})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{g}$
\item $\vect{a}_g = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_g\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{g} = e^{(1-a_{13})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{g}$
% Step 14 
\item $x^0 = x^0 + (a_{14}-a_{13})h$
\item $\vect{a}_h = \vect{a} + b_{14,1}\vect{a}_a + b_{14,6}\vect{a}_i + b_{14,7}\vect{a}_j + b_{14,8}\vect{a}_b + b_{14,9}\vect{a}_c + b_{14,10}\vect{a}_d + b_{14,11}\vect{a}_e + b_{14,12}\vect{a}_f + b_{14,13}\vect{a}_g$
\item $\vect{a}_{h} = e^{-(1-a_{14})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{h}$
\item $\vect{a}_h = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_h\right],\vect{\xi}\right)\right]$
\item $\vect{a}_{h} = e^{(1-a_{14})h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot \vect{a}_{h}$
% Step 15 and 16 combined to reduce memory use 
\item $\vect{a}_i = \vect{a} + b_{15,1}\vect{a}_a + b_{15,6}\vect{a}_i + b_{15,7}\vect{a}_j + b_{15,8}\vect{a}_b + b_{15,9}\vect{a}_c + b_{15,10}\vect{a}_d + b_{15,11}\vect{a}_e + b_{15,12}\vect{a}_f + b_{15,13}\vect{a}_g + b_{15,14}\vect{a}_h$
\item $\vect{a}_j = (1-b_{16,6}/b_{15,6})\vect{a}  + (b_{16,1}-b_{15,1}b_{16,6}/b_{15,6})\vect{a}_a + (b_{16,6}/b_{15,6})\vect{a}_i + (b_{16,7}-b_{15,7}b_{16,6}/b_{15,6})\vect{a}_j + (b_{16,8}-b_{16,8}b_{16,6}/b_{15,6})\vect{a}_b + (b_{16,9}-b_{15,9}b_{16,6}/b_{15,6})\vect{a}_c + (b_{16,10}-b_{15,10}b_{16,6}/b_{15,6})\vect{a}_d + (b_{16,11}-b_{15,11}b_{16,6}/b_{15,6})\vect{a}_e + (b_{16,12}-b_{15,12}b_{16,6}/b_{15,6})\vect{a}_f + (b_{16,13}-b_{15,13}b_{16,6}/b_{15,6})\vect{a}_g + (b_{16,14}-b_{15,14}b_{16,6}/b_{15,6})\vect{a}_h$
\item $x^0 = x^0 + (a_{15}-a_{14})h$
\item $\vect{a}_i = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_i\right],\vect{\xi}\right)\right]$
\item $x^0 = x^0 + (a_{16}-a_{15})h$
\item $\vect{a}_j = h \:\mathcal{F}\left[\vect{\mathcal{N}}\left(x^0,\mathcal{F}^{-1}\left[\vect{a}_j\right],\vect{\xi}\right)\right]$
% Take full step
\item $\vect{a}_{Int} = \vect{a}$
\item $\vect{a} = \vect{a} + c_{1}\vect{a}_a + c_{8}\vect{a}_b + c_{9}\vect{a}_c + c_{10}\vect{a}_d + c_{11}\vect{a}_e + c_{12}\vect{a}_f + c_{13}\vect{a}_g + c_{14}\vect{a}_h + c_{15}\vect{a}_i + c_{16}\vect{a}_j$
\item $\vect{a}_a = \vect{a}_{Int} + c_{1}^*\vect{a}_a + c_{8}^*\vect{a}_b + c_{9}^*\vect{a}_c + c_{10}^*\vect{a}_d + c_{11}^*\vect{a}_e + c_{12}^*\vect{a}_f + c_{13}^*\vect{a}_g + c_{14}^*\vect{a}_h + c_{15}^*\vect{a}_i + c_{16}^*\vect{a}_j$
\end{enumerate}
The method requires 11 copies of the field and has a variable time step which is stochastically safe like its explicit picture counterpart. Once again no rotations are required for steps 15 and 16 since they both occur at the end of the time step.

\subsection{Adding a Cross Vector} 
\label{sec:crossVector}

The nonlinear functionals of Equations~(\ref{eq:exPde})
and~(\ref{eq:ipPde}) may also include variables that are a function of
the same transverse space, but are which are governed by ODEs in one
of the transverse dimensions. These components need never be
transformed to Fourier space, and so ought to be separated from the
main vector components for efficiency reasons. In \eqn{eq:xmdsPdeEx}
they are referred to as the vector $\vect{b}$, which propagates in
dimension $x^c$. This vector has a definite value once the main vector
$\vect{a}$ has a definite value. This means that for most integration
algorithms a current value for the cross vector may be calculated
immediately prior to calculating the $\vect{\mathcal{N}}$ functionals.

In the semi-implicit method in the explicit picture the cross vector
$\vect{b}$ is calculated using the semi-implicit method and using
linear interpolation in the $x^c$ dimension to create the mid-lattice
point values for $\vect{a}$ that are required.

Similarly, with the RK4 method, in both the explicit and interaction
pictures, the calculation of $\vect{b}$ is performed using the RK4
method, and again using linear interpolation in the $x^c$ dimension to
create the mid-lattice point values for $\vect{a}$ that are
required. Note that this use of simple linear interpolation causes a
loss of order when implementing a higher order method such as the RK4,
but generally the transverse lattice is sufficiently fine for this to
be negligible.

However, as mentioned earlier the semi-implicit method in the
interaction picture need only sweep through the memory space of main
field vector once with each nonlinear step. Therefore it would be
desirable to use an algorithm that can calculate the cross vector
simultaneously in the same sweep. The semi-implicit method in the
interaction picture does this by employing the standard semi-implicit
algorithm in the propagation direction {\em simultaneously} with
another semi-implicit algorithm in the transverse propagation
dimension. In pseudo-code it looks like this:
\begin{enumerate}
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\item Calculate $\vect{\xi}$ if required.
\item $x^0 = x^0 + \frac{1}{2}h$
\item For each point in space do:
  \begin{enumerate}
  \item $\vect{a}_{I}=\vect{a}$
  \item $\vect{b}_{I}=\vect{b}$
  \item For N-1 iterations do:
    \begin{enumerate}
    \item $\vect{a} = \vect{a}_{I} + \frac{1}{2} h \:
    \vect{\mathcal{N}}(\vect{x},\vect{a},\vect{b},\vect{\xi})$
    \item $\vect{b} = \vect{b}_{I} + \frac{1}{2} \Delta x^c \:
    \vect{\mathcal{H}}(\vect{x},\vect{a},\vect{b})$
    \end{enumerate}
  \item $\vect{a} = \vect{a}_{I} + h \:
  \vect{\mathcal{N}}(\vect{x},\vect{a},\vect{b},\vect{\xi})$
  \item $\vect{b} = \vect{b}_{I} + \frac{1}{2} \Delta x^c \:
  \vect{\mathcal{H}}(\vect{x},\vect{a},\vect{b})$
  \item $\vect{b}(x^c+\Delta x^c) = \vect{b}_{I} + \Delta x^c \:
  \vect{\mathcal{H}}(\vect{x},\vect{a},\vect{b})$
  \end{enumerate}
\item $x^0 = x^0 + \frac{1}{2}h$
\item $\vect{a} =
e^{\frac{1}{2}h\vect{\mathcal{L}}\left(x^0,\vect{k}_{\bot}\right)} \cdot
\vect{a}$
\end{enumerate}
This routine assumes that the cross vector $b$ was already initialised
at the first lattice point in $x^c$.

\section{Discretisation and Sampling Errors} 
\label{sec:discretisationErrors}

In order to produce an accurate result, any numerical integration must
proceed with a sufficiently small time step $h$. The difference
between the correct result and the numerical result with a particular
time step is known as the {\em discretisation error}. For
deterministic equations, it is usually (see next paragraph) sufficient
to reduce $h$ until the result converges to the desired accuracy.  For
stochastic equations this is not precisely correct, as a greater
number of time steps will induce a different choice ``noise'' from the
random number generators. This means that the two trajectories will be
different members of the ensemble of possible trajectories, so they
should not be expected to converge to each other.  In order to verify
the convergence of a stochastic integration, it is necessary to ensure
that the same trajectory, or {\em path}, is examined with time steps
of differing size.  This can be done by averaging the noise
contributions of a fine-grained path to generate the noise
contributions to a coarse-grained path.

Many problems represent some localised interaction in a psuedo-infite
space. In other words any radiation shed from the interaction should
propagate away without ever being reflected back into the region of
interest. In such cases boundaries must be sufficiently far away from
the important area of the field to have no effect on the result, and
when periodic boundary conditions are used (these are in fact enforced
with \xmds) damping must be applied near the boundaries so as to
absorb any such radiation. The lattice used in each transverse
dimension must be sufficiently fine to be able to exhibit the details,
yet not so fine as to cause memory problems - or cause the solution to
take forever to compute. Further still, field details that are
significantly smaller than those present at initialisation may evolve
during the simulation, and the evolution of such detail will
ultimately be affected by the finesse of the transverse
lattice. Therefore simply reducing the time step alone is not a
complete guarantee of accuracy. The results must be inspected to
ensure that the transverse lattice was fine enough to contain the
detail. If unsure then re-run the simulation with a finer lattice, as
the evolution may be attempting to generate singularities.

Usually, meaningful results only come from stochastic integration by
averaging moments of the fields over large numbers of paths. This
leads to a second form of error for stochastic equations due to the
standard error in the mean. This is known as the {\em sampling error},
and scales as the inverse square root of the number of paths taken.

