\documentclass[11pt]{article}
\usepackage{fullpage,amsmath}
\newcommand{\T}{^{\mathsf{T}}} % matrix transposition
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\bfg}{\mathbf{g}}
\newcommand{\bfx}{\mathbf{x}}
\newcommand{\bflambda}{\mathbf{\lambda}}
\newcommand{\sfE}{\mathsf{E}}
\begin{document}
\begin{center}
{\Large COMPUTER SCIENCE 614}\\
{\Large Numerical Solution of Ordinary Differential Equations} \\
{\large SPRING 2012} \\
ASSIGNMENT \# 4 (34 points) \\
{\small February 24}
\end{center}

\paragraph{Announcement}
\begin{itemize}
\item[] The second examination is Friday, March 2 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 1--5 of the class notes.
 Questions are based on Assignments \#1 through \#3.
\end{itemize}

\paragraph{Due Friday, March 9 in class}
\begin{enumerate}

\item (3 points) \label{it:init}
 For the pendulum example given in Section~5.1 of the class notes,
 write it as a first order system and solve for $\lambda$.
 Use your solution to determine consistent initial values,
 given $X(0) = 0.6$, $X'(0) = 0$, $Y(0) < 0$.

\item (11 points)  \label{it:bdf}
 This problem is an (unsuccessful) attempt to solve the pendulum  example
 using direct discretization.
 Use $m = 1$ and $g = 32$.
 \begin{enumerate}
 \item
 Implement a MATLAB/Octave function
  {\small\begin{verbatim}
function ynew = bdf(h, y, F_handle)
\end{verbatim}}\noindent
  which applies one step of a 2-step BDF method with step size {\tt h}
  to a DAE $F(y, y') = 0$.
  If {\tt y} has 2 columns containing $y_{n}$ and $y_{n+1}$,
  {\tt ynew} will contain $y_{n+2}$.
  If {\tt y} has a single column, the backward Euler method is used.
  The argument {\tt F\_handle} is a handle for a function having a calling
  sequence of the form
  {\small\begin{verbatim}
function [F, dFdy, dFdyp] = F_etc(y, yp)
 \end{verbatim}}\noindent
  which returns not only the value of $F$ but also the square Jacobian
  matrices $\partial F/\partial y$ and $\partial F/\partial y'$.
  Use a Newton-Raphson solver as described in Problem~2(a) of Assignment \#~3.
  The error estimate mentioned there
  is the max norm of the difference between two successive iterates.
  Have your function print the error estimates as well as $y_n$.
  Include your code with your submission.
 \item Write a function {\tt F\_etc} for the pendulum DAE.
  Include your code with your submission.
 \item Apply {\tt bdf} to {\tt F\_etc} with step size 0.002 for five steps.
  Use your solution to Problem~\ref{it:init} as initial values.
  Include the output with your submission.
 \end{enumerate}

\newpage
\item (3 points)
 Following is the output for Problem~\ref{it:bdf} for other step sizes
 using code given in the appendix:
 {\small\begin{verbatim}
h =  0.010000
err =  22.400
err =  0.016415
err =  4.0782e-07
err =  2.8023e-13
solution =

   0.59795  -0.80154  -0.20539  -0.15363   9.58358

err =  22.531
err =  0.016651
err =  1.0333e-07
err =  7.2742e-13
solution =

   0.59451  -0.80409  -0.41252  -0.30583   9.41903

err =  0.020004
err =  0.029514
not converging
-----------------
h =  0.0050000
err =  22.400
err =  0.0040980
err =  6.3890e-09
solution =

   0.599488  -0.800384  -0.102474  -0.076804   9.595902

err =  22.433
err =  0.0041126
err =  1.5898e-09
solution =

   0.59863  -0.80102  -0.20517  -0.15343   9.55490

err =  0.0048348
err =  0.0073059
not converging
\end{verbatim}}\noindent
What anomalies do you observe in the output?
How does this correlate with the very large initial {\tt err}?
What essential theoretical property does the numerical method seem to lack?

\newpage
\item (6 points)
 Let $0 = t_0 < t_1 < \cdots < t_{k-1} < t_{k+1} < \cdots < t_N$.
 Let $W_i$, $i\ne k$ and $0\le i\le N$, be Gaussian
 random variables satisfying $\sfE[W_i] = 0$ and
 $\sfE[W_iW_j] = \min\{t_i, t_j\}$.
 Let $Z$ be a standard Gaussian random variable independent of
 $W_i$, $i\ne k$, $0\le i\le N$.
 Define
 \begin{equation}  \label{eq:one}
 W_k = aW_{k-1} + bW_{k+1} + \sqrt{ab}\sqrt{t_{k+1} - t_{k-1}} Z
 \end{equation}
 where $a = (t_{k+1} - t_k)/(t_{k+1} - t_{k-1})$
 and $b = (t_k - t_{k-1})/(t_{k+1} - t_{k-1})$.
 Determine $\sfE[W_k]$ and $\sfE[W_kW_j]$, $j = 0, 1, \ldots, N$.
 And of course simplify!
 What does formula~(\ref{eq:one}) do for us?

\item (1 point)
 Express the Riemann-Stieltjes integral $\int_0^T f(t)\rmd W(t)$
 as a Riemann integral by integrating by parts under the assumption
 that $f\in\mathrm{C}^1[0,T]$.

\item (5 points)
 How would one generate a random number having the same distribution
 as $\int_0^T W(t)\rmd t$ given the ability to generate a
 standard Gaussian $Z$?  Hint: integrate by parts to get something
 more familiar.

\item (5 points)
 Do Exercise~16.7 of the textbook.

\end{enumerate}

\subsection*{Appendix}

Here is an Octave example of a driver for Problem~2:
{\small\begin{verbatim}
y = [0.6; -0.8; 0; 0; -12.8];
N = 5;
h = 0.01
try
   for n = 2:N+1
      ynew = bdf(h, y, @F_etc);
      y = [y(:, end), ynew];
   end
catch
   disp(lasterror().message)
end
y = [0.6; -0.8; 0; 0; -12.8];
disp('-----------------')
h = 0.005
try
   for n = 2:N+1
      ynew = bdf(h, y, @F_etc);
      y = [y(:, end), ynew];
   end
catch
   disp(lasterror().message)
end
\end{verbatim}}\noindent

\end{document}