\documentclass[11pt]{article}
\usepackage{fullpage,amsmath}
\newcommand{\T}{^{\mathsf{T}}} % matrix transposition
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\bff}{\mathbf{f}}
\newcommand{\bfy}{\mathbf{y}}
\begin{document}
\begin{center}
{\Large COMPUTER SCIENCE 614}\\
{\Large Numerical Solution of Ordinary Differential Equations} \\
{\large SPRING 2012} \\
ASSIGNMENT \#6 (29 points) \\
{\small April 2}
\end{center}

\paragraph{Announcement}
\begin{itemize}
\item[] The second examination is Friday, April 20 in class
 (from  12:30 pm to 1:20 pm).  It is a closed book exam, and only pens,
 pencils, and erasers are permitted.
 It covers Chapters 6--9 of the class notes.
 Questions are based on Assignments \#4 through \#6.
\end{itemize}

\paragraph{Due Friday, April 13 in class}
This assignment covers Chapter 8 and Sections  9.3--9.5 of the class notes,
which corresponds to Sections 9.1--9.8, 10.1--10.3, 11.0, 11.1, 11.3, 11.4,
15.0--15.5 of the textbook.

\begin{enumerate}

\item (2 points)  % 8.2
 Do Exercise 14.2 from Griffiths \& Higham,
 but only for the case of Runge-Kutta methods.

\item (8 points)  % 8.6
 For a system of ODEs $\bfy' = \bff(\bfy)$ in autonomous form,
 the theta method is given by
 \[
 \bfy_n =\bfy_{n-1} + h_n(\theta\bff(\bfy_n) + (1 -\theta)\bff(\bfy_{n-1}))
 \]
 where $\theta$ is a method parameter in the range $0\le\theta\le 1$.
 \begin{enumerate}
 \item (6 points)
  For which values of $\theta$ is the method A-stable?
  It follows from the maximum modulus theorem (of complex variables)
  that if a function is analytic in the left half plane,
  it attains its maximum modulus on the imaginary axis.
 \item (2 points)
  Determine the Butcher tableau (Butcher array) for the theta method
  (as an implicit Runge-Kutta method).
 \end{enumerate}

\item (19 points)
 For the exercises below,
 you are expected to abide by the stated specifications.
 Deviations should be explained and may still be penalized.
 Hand in your code for parts (a)--(c), your plots for part (d), and
 your answers for parts (e) and (f).
 \begin{enumerate}
 \item (4 points)
  Consider the Hamiltonian $\frac12 M^{-1} p^2 + U(q)$ where $U(q)$ is
  the Morse potential given by
  \[ U(q) = D(1 -\mathrm{e}^{-S(q - q_0)})^2 \]
  with $D = 90.5\times 0.0004814$, $S = 1.814$, $q_0 = 1.41$,
  and $M = 0.9953$.
  Here $q$ represents the distance between a pair of bonded atoms.
  Write a MATLAB function {\it Ham}, which given a 2-dimensional array
  $[q, p]\T$ returns the value of the Hamiltonian $H(q,p)$
  (for the Morse potential).
  Write a {\it rhs} function,
  which given a scalar position $q$ and momentum $p$,
  returns the velocity and the force.
 \item (4 points)
  Write a {\it RKstep} function which takes one step of Kutta's Simpson
  method. It will take as arguments a 2-dimensional array
  $[q_{n-1}, p_{n-1}]\T$ and a stepsize $h$ and return a 2-dimensional array
  $[q_n, p_n]\T$. 
 \item (4 points)
  In the class notes, find a formula for the leapfrog method that
  expresses $q_n$, $p_n$ in terms of $q_{n-1}$, $p_{n-1}$.
  Implement it by writing an {\it LFstep} function
  which takes one step of the leapfrog method.
  It will take as arguments a 2-dimensional array
  $[q_{n-1}, p_{n-1}]\T$ and a stepsize $h$ and return a 2-dimensional array
  $[q_n, p_n]\T$.
  Do not use half-step values as arguments or return values.
 \item (5 points)
  Using initial values $q(0) = 1.4155$, $p(0) = 1.545 M/48.888$, 
  run {\it RKstep} for 250 steps of size $h = 8$, 
  and plot the deviation in the Hamiltonian from its initial value.
  Run {\it LFstep} for 1000 steps of size $h = 2$, 
  $h = 2.3684$, $h = 2.3685$,
  and, for each of these, plot the deviation in the Hamiltonian.
 \item (1 point) How does the leapfrog method compare to Kutta's Simpson method
  for equal work?
 \item (1 point) What do you observe when the step size changes
 for the leapfrog method from $h = 2.3684$ to $h = 2.3685$?
 \end{enumerate}

\end{enumerate}

\end{document}