\documentclass[10pt,letterpaper]{article}
\usepackage{graphicx,amssymb}
\usepackage[colorlinks=true]{hyperref}

\topmargin -0.5in
\textheight 9in
\oddsidemargin=-0.05in
\evensidemargin=-0.05in
\textwidth 6.5in
\pagestyle{empty}

\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand\ben{\begin{enumerate}}
\newcommand\een{\end{enumerate}}
\newcommand\vg{\vspace{2ex}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}
\newcommand{\ie}{{\it i.e.\ }}
\newcommand{\eg}{{\it e.g.\ }}
\newcommand{\eps}{\varepsilon}
\newcommand{\mbf}[1]{{\mathbf #1}}
\newcommand{\tbox}[1]{{\mbox{\tiny #1}}}
\newcommand{\emach}{\eps_\tbox{mach}}


\begin{document}
\title{Math 56 Compu \& Expt Math, Spring 2014: Homework 3}
\date{\vspace{-2ex} due 10am Thursday April 17th}
\maketitle
\ben

\item Here you learn how to ``roll your own'' finite difference formulae.
Let's say you have access to $f$ at only $x$, $x+h$, and $x+2h$,
and want a 2nd-order accurate approximation for $f'(x)$.
Note that this is at the leftmost point of the three; e.g.\ at the extreme
end of a grid of values.
\ben
\item
We want $f'(x) \approx a f(x) + b f(x+h) + c f(x+2h)$.
Our goal is to solve for the coefficients $a$, $b$, $c$,
by setting up a 3-by-3 linear system for them, as follows.
Expand the right-hand side via Taylor series about $x$.
Since you want this to hold at all $h$,
you may extract just the constant terms to give one equation, just the
coefficients of $h$ to give a 2nd equation, and just the coefficients
of $h^2$ for the 3rd equation.

%[\LaTeX\ hint:
%\verb+ \left[\begin{array}{lll} x & y & z \\ w ...  \end{array}\right] +
%]

\item Solve the system either by hand or computer, hence write
your new finite difference formula. How do you know the solution is unique?

\item Give a {\em rigorous} upper bound on the error of this formula
(in exact arithmetic, i.e.\ ignore rounding).
\een


\item Stability.
\ben
\item Show whether subtraction $x_1-x_2$ is backwards stable
(with respect to the two input data) under the
rules of floating point.
\item
Cosine, like anything else done by machine, cannot be more accurate
than a relative error of $\emach$ in its output (ie ``forward error'').
Explain if a machine implementation of cos could be backwards stable near
$x=0$ ? Near $\pi/2$ ?
%In a worksheet you found that $1+x$ as done by the rules of
%floating point arithmetic is not backward stable.
%Show whether $1+x$ is {\em stable} or not.

% decided to cut stability!
\een


\item Here's a new formula for matrix 2-norm:
$$\|A\| = \sqrt{\lambda_\tbox{max}(A^T A)}, \qquad \mbox{
where $\lambda_\tbox{max}(A^T A)$
is the largest eigenvalue of the matrix $A^T A$.}$$
\ben
\item
Let $A = \left[\begin{array}{ll}1&-1\\2&2\end{array}\right]$.
Use the new formula to compute by hand $\|A\|$.
How does it compare to the size of the largest eigenvalue of $A$?
(for which you can use {\tt eig})
\item Use this to compute the matrix condition number $\kappa(A)$.
Is it well-conditioned?
\item Take 100 points $\mbf{x}\in\mathbb{R}^2$
equi-spaced on the unit circle, and plot them, and $A\mbf{x}$ for each.
What geometric property does $\kappa(A)$ measure of the ellipse produced?
\een



\item
Download the two $100\times100$ matrices {\tt A1} and {\tt A2} from the HW page,
and use {\tt textread} to read them into Matlab (you will need to {\tt
reshape} them).
\ben
\item Compare their matrix 2-norms and condition numbers.
What worst-case relative errors do you expect for solving linear systems
with matrix {\tt A1}? With {\tt A2}?
(Use our backward stability theorem, and assume standard double precision.)

\item Let's focus on $A = $ {\tt A1}, and load in the RHS
$\mbf{b} = $ {\tt bvec} from the HW page.
Solve $A\mbf{x} = \mbf{b}$.
Then perturb $\mbf{b}$ by a random vector of norm $\eps_\tbox{mach}$
to get $\tilde{\mbf{b}}$
(this emulates rounding error applied to the RHS),
and solve again $A\tilde{\mbf{x}} = \tilde{\mbf{b}}$.
What relative norm change $\|\tilde{\mbf{x}}-\mbf{x}\|/\|\mbf{x}\|$ results?
Does this match your prediction from (a)?
\item
Repeat (b) except using the RHS $\mbf{c} = $ {\tt cvec} from the HW page.
Surprising?
Is it consistent with (a)? Repeat for random unit-norm RHS vectors---do they
behave more like $\mbf{b}$ or like $\mbf{c}$?

\item [BONUS] Explain the different behaviors [hint: $\|\mbf{x}\|$],
deducing how the
directions of $\mbf{b}$ and $\mbf{c}$ relate to long and short axes of
the ellipse of the image of the unit sphere under $A$.

\item Given $A\in\mathbb{R}^{M\times P}$ and $B\in\mathbb{R}^{P\times N}$,
prove a bound on $\|AB\|$ in terms of the norms of the individual matrices.
[Hint: HW2 6(c).]

\een


\item
\ben
\item {[Trefethen and Bau, Ex.\ 13.3]}
On the same axis, plot $p(x) = (x-2)^9$ around its root two ways. First plot it by evaluating $p(x)$ via its factored form, then plot it via its expanded form $p(x) = x^9-18 x^8+144x^7-672x^6+2016x^5-4032x^4+5376x^3-4608x^2+2304x-512$. Evaluate it at the points $x = 1.920, 1.921, 1.922,...,2.08$. Explain the discrepancy.
\item
Compute the roots of $q(x) = x^2-\frac{100000001}{10000}x+1$ using the quadratic formula (in \verb|Matlab|). Notice that $q(x) = (x-10000)(x-\frac{1}{10000})$. What is the relative error of the computed roots? Explain why one root is good to only 9 digits of accuracy.
Can you find a way to get this root
to close to machine accuracy using only the input coefficients $a$, $b$ and $c$
and the other root?
\een

\een
\end{document}
