\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 $$\label{eq:one} W_k = aW_{k-1} + bW_{k+1} + \sqrt{ab}\sqrt{t_{k+1} - t_{k-1}} Z$$ 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}