\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{\vphantom\vdots\\}

\begin{document}

\title{Numerical Methods II}
\author{Problem Sets 1/2}
\date{due in class, February 23, 2004}
\maketitle


\begin{enumerate}

\item (From Gautschi, 1997, p.\ 324.)
\begin{enumerate}
\item Show that any one-step method of order $p$, when applied to the
linear model problem $\y' = A\y$, yields
\[
  \y_{n+1} = \phi(hA) \, \y_n \,,
\]
where
\[
  \phi(z) = 1 + z + \tfrac1{2!} \, z^2 + \dots + \tfrac1{p!} \, z^p
  + O(z^{p+1}) \,.
\]

\item Show, in particular, that the $O(z^{p+1})$ terms vanish for a
$p$-stage Runge--Kutta method of order $p=1,\dots,4$, and for the
Taylor method of order $p \geq 1$.
\end{enumerate}

\item (From Gautschi, 1997, p.\ 328.) For an analytic function $f$,
show that the level curve $|f(z)|=1$ obeys the differential equation
\[
  \frac{\d z}{\d s} = \i \, \frac{f(z) \, |f'(z)|}{f'(z)} \,,
\]
where $s$ is the arclength parameter along the level curve.

\emph{Hint:} Write $f(z) = r \, \e^{\i \theta}$ and take $\theta$ as
the independent variable.  In a second step, change the independent
variable from $\theta$ to $s$.

\item \textbf{Project:} Use the results from the previous two
exercises to plot the boundary of the region of absolute stability in
the complex $z=\lambda h$ plane for the Taylor methods of order
$p=1,\dots,5$.

\emph{Notes:} You can use, for example, the integrator
\verb+ode_rk4+ from last semester.  The origin $z=0$ is always on
the boundary of the region of absolute stability.  (Explain!)

\item \emph{Shift of the Taylor polynomial.}
\begin{enumerate}
\item Let $p(t)$ be a polynomial of degree $m$, and let $\v$ be the
vector of coefficients of its Taylor expansion about $t=0$, scaled
such that
\[
  p(t) = v_0 + v_1 \, \frac{t}{h} + \dots + v_m \,
         \Bigl( \frac{t}{h} \Bigr)^m \,.
\]
Show that the scaled Taylor coefficients of $p$ about the point $t=h$,
i.e.\ the coefficients defined via
\[
  p(t) = w_0 + w_1 \, \frac{t-h}{h} + \dots + w_m \,
         \Bigl( \frac{t-h}{h} \Bigr)^m \,,
\]
are given by $\w = S\v$, where the components of $S$ are
\[
  S_{ij} =
    \begin{cases}
    0 & \text{for } i>j \\
    \binom{j}{i} & \text{for } i \leq j
    \end{cases}
\]
and $i,j = 0,\dots,m$.

\item
Show that the inverse shift has coefficients
\[
  (S^{-1})_{ij} =
    \begin{cases}
    0 & \text{for } i>j \\
    (-1)^{j-i} \, \binom{j}{i} & \text{for } i \leq j \,.
    \end{cases}
\]
\end{enumerate}

\item \emph{BDF method in Nordsieck's representation.}  Recall that
BDF methods are implicit linear multistep methods obtained through
approximating the left side of the equation $y'=f(t,y)$ by a
polynomial interpolant.

More precisely, for the $m$-step, order $m$ BDF method in the scalar
case, let $p_n$ denote the interpolating polynomial at timestep $t_n$.
Then
\begin{gather*}
  p_n(t_{n-j}) = y_{n-j} \quad \text{for } j=0,\dots,m \,, \\
  p_n'(t_n)= f_n \equiv f(t_n,y_n) \,.
\end{gather*}
Thus, computing the numerical solution $y_n$ at time $t=t_n$ requires
knowledge of $m$ previous values $y_{n-m}, \dots, y_{n-1}$.

Nordsieck's idea is that instead of storing previous values, one may as
well store the scaled Taylor coefficients of the interpolating
polynomial, namely
\[
  \z_n = 
  \begin{pmatrix}
    p_n (t_n) \\
    p_n'(t_n) \, h \\
    \vdots \\
    p_n^{(m)}(t_n) \, h^m/{m!}
  \end{pmatrix} \,.
\]

\begin{enumerate}
\item Show that
\[
  p_{n+1}(t) - p_n(t)
  = (y_{n+1}-p_n(t_{n+1})) \,
    \phi \Bigl( \frac{t-t_{n+1}}h \Bigr) \,,
\]
where
\[
  \phi(z) = \prod_{j=1}^m \Bigl( \frac{z}{j}+1 \Bigr) \,.
\]

\item Show that the BDF method of order $m$ is equivalent to
\[
  \z_{n+1} = S \z_n + \v \, \bigl( h \, f_{n+1}
             - \boldsymbol{e}_1^T S \z_n \bigr) \,,
\]
where $S$ is the shift matrix from Question 4, $\v$ is the vector of
(unscaled) Taylor coefficients of the polynomial $\phi(z)/\phi'(0)$,
and $\boldsymbol{e}_1 = (0,1,0,\dots,0)^T$.

\item Describe how you would change the step size in Nordsieck's
representation of the BDF method.

\end{enumerate}

\item \textbf{Project:} Implement BDF4 with step size
control for a system of linear equations in Nordsieck's representation
as an \textsf{Octave} function.  The solver should take arguments of
the form
\begin{verbatim}
  ode_bdf4_var ('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.

Estimate the local error by using a predictor-corrector pair, where
the predictor is the extrapolated value $p_n(t_{n+1})$, and the
corrector is given by the result of the BDF step.  Use Broyden's
method to solve the system of nonlinear equations that arises at each
time step.

\item \textbf{Project:} Use your variable time-step
BDF method to solve 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$.  Can you get up to
$\mu=1000$ within a reasonable computation time?

\end{enumerate}

\end{document}
