\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts}

\newcommand{\x}{\boldsymbol{x}}
\newcommand{\y}{\boldsymbol{y}}
\newcommand{\f}{\boldsymbol{f}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\r}{\boldsymbol{r}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\renewcommand{\k}{\boldsymbol{k}}
\newcommand{\z}{\boldsymbol{z}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\renewcommand{\d}{\mathrm{d}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}

\def\nl{$\displaystyle\vphantom\int$\\}

\begin{document}

\title{Numerical Methods II}
\author{Lab Session 1}
\date{February 9, 2004}
\maketitle


\begin{enumerate}

\item Implement an adaptive Runge--Kutta method using the embedded
Prince--Dormand pair of order 4,5 with coefficient table
\begin{center}
\begin{tabular}{|c|c|ccccccc|}
\hline
$i$ & $a_i$ & \multicolumn{7}{c|}{$b_{ij}$ for $j=1,\dots,i-1$} \nl 
\hline
\hline
$1$ & $0$ &&&&&&&\nl
$2$ & $\tfrac15$ & $\tfrac15$ &&&&&&\nl
$3$ & $\tfrac3{10}$ & $\tfrac3{40}$ & $\tfrac{9}{40}$
      &&&&&\nl
$4$ & $\tfrac45$ & $\tfrac{44}{45}$ & $-\tfrac{56}{15}$ &
      $\tfrac{32}{9}$ &&&&\nl
$5$ & $\tfrac89 $ & $\tfrac{19372}{6561}$ & $-\tfrac{25360}{2187}$ & 
      $\tfrac{64448}{6561}$ & $-\tfrac{212}{729}$
      &&&\nl
$6$ & $1$ & $\tfrac{9017}{3168}$ & $-\tfrac{355}{33}$ &
      $\tfrac{46732}{5247}$ & $\tfrac{49}{176}$ &
      $-\tfrac{5103}{18656}$ &&\nl
$7$ & $1$ & $\tfrac{35}{384}$ & ${0}$ & $\tfrac{500}{1113}$ &
      $\tfrac{125}{192}$ & $-\tfrac{2187}{6784}$ &
      $\tfrac{11}{84}$ &\nl
\hline
\hline
\multicolumn{2}{|c|}{$c_j$} &
      $\tfrac{5179}{57600}$ & ${0}$ & $\tfrac{7571}{16695}$ &
      $\tfrac{393}{640}$ & $-\tfrac{92097}{339200}$ &
      $\tfrac{187}{2100}$ & $\tfrac{1}{40}$ \nl
\hline
\multicolumn{2}{|c|}{$\tilde c_j$} &
      $\tfrac{35}{384}$ & ${0}$ & $\tfrac{500}{1113}$ &
      $\tfrac{125}{192}$ & $-\tfrac{2187}{6784}$ &
      $\tfrac{11}{84}$ & ${0}$ \nl 
\hline
\end{tabular}
\end{center}
where
\[
  \k_i = \f\bigl( t_n + a_i h, \y_n + h (b_{i1} \, \k_1 + \dots +
                  b_{i,i-1} \, \k_{i-1}) \bigr) 
\]
for $i=1,\dots,7$, and
\begin{gather*}
  \y_{n+1} = \y_n + h \,
             \bigl( c_1 \, \k_1 + \dots + c_7 \, \k_7 \bigr) \,, \\
  \tilde \y_{n+1} = \y_n + h \, \bigl( \tilde c_1 \, \k_1 +
             \dots + \tilde c_7 \, \k_7 \bigr) \,.          
\end{gather*}

The solver should take arguments of the form
\begin{verbatim}
  [y,t] = ode_rk54 ('f', [t0,t1], y0, tol, Nmax)
\end{verbatim}
where \texttt{f} is the function $f(t, y)$ on the right hand side of
the differential equation, \texttt{t0} is the initial time,
\texttt{t1} is the final time, \texttt{y0} the initial condition,
\texttt{tol} the local error tolerance and \texttt{Nmax} the maximum
number of steps.  The return valus \texttt{[y,t]} are the matrix of
solution values $(\y_0, \dots, \y_N)$ and the row vector of evaluation
times $(t_0, \dots, t_N)$, respectively.

Test your method with the van der Pol equation
\begin{align*}
  \dot x & = y \,, \\
  \dot y & = \mu \, (1-x^2) \, y - x \,,
\end{align*}
with initial data $x (0) = 2$ and $y(0) = 0$, and $\mu=10$.

\item Modify your code for the implicit BDF4 solver to use Broyden's
method rather than simple iteration to solve the system of nonlinear
equations that arises at each time step.

Recall that Broyden's method for iteratively solving $\f(\x)=0$
constructs---in addition to the sequence $\x_k$ that (hopefully)
converges to the solution $\x$---a sequence of approximate Jacobians
of $\f$ via
\begin{gather*}
  B_k \, \Delta \x_k = - \f(\x_k) \,, \\
  \x_{k+1} = \x_k + \Delta \x_k \,, \\
  \Delta \f_k = \f(\x_{k+1}) - \f(\x_k) \,,\\
  B_{k+1} = B_k + \frac{(\Delta \f_k - B_k \, \Delta \x_k) \, \Delta
                        \x_k^T}{\Delta \x_k^T \, \Delta \x_k} \,.
\end{gather*}

Test your code with the van der Pol equation.  How large time steps
can you take in the initial part of the computation?  How should your
initialize your multi-step method?

\end{enumerate}

\end{document}
