\documentclass[12pt]{article}
\setlength{\textwidth}{7.5in}
\setlength{\evensidemargin}{-.5in}
\setlength{\oddsidemargin}{-.5in}

\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bd}{\begin{displaymath}}
\newcommand{\ed}{\end{displaymath}}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\bn}{\begin{enumerate}}
\newcommand{\en}{\end{enumerate}}
\newcommand{\pa}{\partial}
\newcommand{\f}{\frac}
\newcommand{\ci}{\cite}
\newcommand{\eps}{\epsilon}
\newcommand{\del}{\delta}

\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{epic}
\usepackage{eepic}
\usepackage{amssymb}
\usepackage{afterpage}

\setlength{\unitlength}{1.0cm}


\begin{document}

\begin{center}
Homework 4 \\
MA/CS 375, Fall 2005 \\
Due November 11 \\
\end{center}

This homework will count as part of your grade so you must work independently.
It is permissible to discuss it with your instructor, the TA, fellow students,
and friends. However, the programs/scripts and report must be done only
by the student doing the project. Please follow the guidelines in the syllabus
when preparing your solutions.

\begin{enumerate}

\item (25 pts.) In this problem we approximate the derivative of a function
defined on a uniform grid on the interval $\left(a, b\right)$,
i.e. for the grid
\[x_k = a + k\frac{(b-a)}{n}\ ,\ k=0,1,\cdots,n \ \]
\begin{enumerate}
\item Show (using Taylor series) that the formula
\[
f'(x_k) - \frac{f(x_{k+1}) - f(x_{k-1})}{2h} =
          \frac{h^2}{6}f'''(\xi_k)   \ ,
          \  \xi_k \in \left(x_{k-1},x_{k+1}\right)\]
is valid for all interior grid points, i.e. $k = 1,2,\cdots,n-1$.
\item
For the end points, use Taylor series to show that
\begin{eqnarray*}
f'(x_0) - \frac{-3f(x_{0}) +4f(x_{1}) - f(x_{2})}{2h} & = &
     -\frac{h^2}{3}f'''(\xi_0) \ ,\ \xi_0 \in \left(x_{0},x_{2}\right) \\
f'(x_n) + \frac{-3f(x_{n}) +4f(x_{n-1}) - f(x_{n-2})}{2h} & = &
     +\frac{h^2}{3}f'''(\xi_n) \ ,\ \xi_n \in \left(x_{n-2},x_{n}\right)
\end{eqnarray*}
\item Write a code that uses the above formulas to produce a second order
accurate approximation for the derivative of the function
\[ f(x) = \cos{(10 x)} - 3\sin{(13 \pi x)} \]
on the interval $-1 \le x \le 1$, with an error no larger than
$tol = 10^{-6}$. You must estimate the third derivative of your function
using a worst case senario, and choose the stepsize $h$ based on your
estimates.
You must compute the exact derivative at the same grid points and produce
a plot of your absolute error.
\end{enumerate}

\item (25 pts.) Modify one of the programs in the text (midpointc,
trapc or simpsonc) to use the composite Newton's $3/8$-rule on a uniform
partition to approximate an integral. On a subinterval $(x_{k-1},x_k)$
of width $h$ the $3/8$-rule is defined by:
\bd
\int_{x_{k-1}}^{x_k} f(x) dx \approx \f {h}{8} \left( f(x_{k-1}) + 3
f(x_{k-1} + h/3) + 3 f(x_{k-1} + 2h/3) + f(x_k) \right) .
\ed
Show that the rule is exact for polynomials of degree $3$ or less and use this
fact to predict the order of the method. Use the method with
$n=[5 \ 10 \ 20 \ 40 \ 60 \ 80 \ 100]$ subintervals to approximate:
\bd
\int_{0}^{1} \sin{ (\pi x)} e^{2x} dx .
\ed
Using the fact that the exact value of the integral is $\f {\pi}{4+\pi^2}
(e^2 +1)$, compute the errors in your approximations and make a loglog plot
of the error versus $n$. Is the result approximately a straight line? Use
polyfit to compute a linear approximation, $\log{(\rm error)} = a \log{(n)} +
b$ to the data. Is the slope, $a$, what you expect from the error analysis?

\item (25 pts.) The four point Gauss-Lobatto formula for approximating
integrals takes the form:
\bd
\int_{x_{k-1}}^{x_k} f(x) dx \approx h \left( w_1 f(x_{k-1}) + w_2
f(x_{k-1} + c_1 h) + w_3 f(x_{k-1} + c_2 h) + w_4 f(x_k) \right) , \ \
h=x_k - x_{k-1} .
\ed
The weights $w_j$, $j=1, \ldots ,4$ and constants $c_1$, $c_2$ can be chosen
so that the formula is exact for polynomials of degree $5$. Find them.
What is the order of the resulting method? \\
Hint: solve the problem on
the standard interval $[-1,1]$ and then translate the results to an arbitrary
interval. The weights $w_j$ can be determined for arbitrary values of $c_j$
using cubic interpolation and will lead to methods exact for cubic polynomials.
You can then determine the two free node locations so that the formula is
exact for $x^4$ and $x^5$. You will need to solve an algebraic equation -
you can do this either symbolically or using fzero. Solve it accurately
enough so that you can give us the $w_j$ and $c_j$ to $14$ digits.

\item (25 pts.) The adaptive quadrature code given in the text, estimates the
error in Simpson's rule by performing two integrations over the interval
$A$. The first uses the simple Simpson rule,
with the nodes $\left[a, (a+b)/2, b\right]$, i.e. with stepsize $H_1 = b-a$
to compute the approximation $I_s$.
The second uses the composite Simpson's rule  with 2 subintervals, i.e.
nodes $\left[a, (a+b)/4, (a+b)/2; (a+b)/2, 3(a+b)/4, b\right]$
and stepsize $H_2 = (b-a)/2$, to compute the (more accurate)
approximation $I_s^c$.

The text uses the two approximations to estimate the error and, if that
turns out to be less than a predefined tolerance, it uses the more
accurate result of the two, i.e. $I_s^c$, as the actual estimate for the
integral in the interval $A$. This problem is about improving this
calculation through the judicious use of both results, $I_s$ and $I_s^c$.

It can be shown that the error in Simpson's
method for integrating a function that is differentiable at least
five times can be actually written as
\[ I[f] - I_s^c[f;M] =
    c_1 *\left(\frac{b-a}{M}\right)^4 + c_2*\left(\frac{b-a}{M}\right)^5
+\cdots
 \ \]
where $c_1$, $c_2$, etc are constants, independent of M
(for functions $f$ that have more
derivatives in the interval of integration, we can add more terms to that
expression, but that is not important here).
The leading term in this expression is of order 4, which is the order of
the method. In this problem we learn how to get a more accurate result
using the work already performed. To that purpose:
\begin{enumerate}
\item
Find a linear
combination of the two results that has leading error term of
order 5; that is you must find constants $A_1$ and $A_2$ such that
\[ I[f] - \left(A_1 I_s^c + A_2 I_s\right) = c_2^\prime \left(\frac{b-a}{M}\right)^5 +\cdots \ \]
The idea is to combine the two results in such a way so that the leading
error term is cancelled.
\item Now alter the program {\bf simpadpt} in the text, p.98, to use
this improved estimate of the integral over $A$ instead of the (less
accurate) estimate $I_s^c$. Demonstrate your code by producing a more
accurate computation of the integral in {\bf example 4.3}. Your computation
should employ the same tolerance ($ tol = 10^{-5}$) and minimum stepsize
($hmin = 10^{-3}$) as the values in the text, and therefore end up using the
same number of subintervals. Nevertheless, your result ought to be closer to
the exact value
\[ \int_{-1}^1 e^{-(10)(x-1)^2} dx = 0.28024956081990 \ ,\]
than the one found by program {\bf simpadpt}.

\end{enumerate}

\end{enumerate}
\end{document}
