\documentclass[12pt]{article}
%\usepackage{ccfonts,euler}
\usepackage{fullpage,amsmath,amsfonts,graphicx}

\newcommand{\x}{\boldsymbol{x}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\r}{\boldsymbol{r}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}

\begin{document}

\title{Numerical Methods I}
\author{Final Exam}
\date{December 19, 2005}
\maketitle

% \vspace*{\fill}
% {\baselineskip=40pt
% \begin{tabbing}
% Student ID Number: \= \kill
% Last Name: \> \rule[-2pt]{10cm}{0.3pt} \\
% First Name:  \> \rule[-2pt]{10cm}{0.3pt} \\
% Signature: \> \rule[-2pt]{10cm}{0.3pt} \\
% %Code: \> \fbox{\rule{0pt}{1.5em}\rule{1.5ex}{0pt}}
% %\fbox{\rule{0pt}{1.5em}\rule{1.5ex}{0pt}}
% %\fbox{\rule{0pt}{1.5em}\rule{1.5ex}{0pt}}
% %\fbox{\rule{0pt}{1.5em}\rule{1.5ex}{0pt}}
% \end{tabbing}}


% \vspace*{\fill}

For $N$ even, the discrete Fourier transform of an $N$-tuple of
numbers $u_0,\dots,u_{N-1}$, is defined
\begin{equation*}
  \hat u_k = \frac1N \sum_{j=0}^{N-1} \e^{-\i k x_j} \, u_j
\end{equation*}
with $x_j = jh$ and $h=2\pi/N$.  The inverse discrete Fourier
transform is given by
\begin{equation*}
  u_j = \sum_{k=-N/2}^{N/2-1} \e^{\i k x_j} \, \hat u_k \,.
\end{equation*}


\begin{enumerate}
\item
Consider the following \textsf{Octave} program.
\begin{quote}
\begin{verbatim}
format long;
epsilon = 1e-14;

A = vandermonde (10); % vandermonde(n) gives the n x n Vandermonde matrix
b = sum (A,2);

[L,U,P] = lu(A);
y = L\(P*b);
x = U\y; 
r = b-A*x; 
[x r]

while norm(r)>epsilon
  v = L\(P*r);
  w = U\v;
  x = x + w; 
  r = b-A*x;
  [x r]
end
\end{verbatim}
\end{quote}

\begin{enumerate}
\item What does this code do?  Explain.


\item If you run the code in Octave, you'll find the following:

\begin{verbatim}
ans =

   0.999661641799626   0.000000000000000
   1.000619472137607   0.000000016163540
   0.999521918673179  -0.000000093248673
   1.000205001753909  -0.000000229105353
   0.999945946146878  -0.000000005587935
   1.000009127129884  -0.000000104308128
   0.999999009142851  -0.000000208616257
   1.000000066933674  -0.000000119209290
   0.999999997438762   0.000000000000000
   1.000000000042420   0.000000000000000

ans =

  0.999992500612691  0.000000000000000
  1.000014060201136  0.000000000000000
  0.999988837483866  0.000000000000000
  1.000004942746758  0.000000000000000
  0.999998650834743  0.000000000000000
  1.000000235987559  0.000000000000000
  0.999999973495598  0.000000000000000
  1.000000001846354  0.000000000000000
  0.999999999927497  0.000000000000000
  1.000000000001225  0.000000000000000
\end{verbatim}
What can you say about the condition of the problem?  Is the
residual a good error indicator?


\item Will this code always terminate?  Explain.
\end{enumerate}
\hfill(10+10+10)


\item Use Newton's method for solving the quadratic equation $x^2=q$
for a given positive real number $q$.

\begin{enumerate}
\item Show that, in this case, the fixed point iteration reads
\begin{equation*}
  x_{k+1} = \frac{x_k}2 + \frac{q}{2 \, x_k} \,.
\end{equation*}
\item Let $x$ denote the exact solution.  Show that
\begin{equation*}
  x_{k+1} - x = \frac{(x_k-x)^2}{2 \, x_k} \,.
\end{equation*}
\item What can you say about the order of convergence?
\end{enumerate}
\hfill(10+5+5)


\item Compute the location of the quadrature points for Gauss
quadrature on the interval $[-1,1]$ with two quadrature points.

\emph{Hint:} you may use the fact that Gauss quadrature integrates
polynomials up to a certain degree exactly.
\hfill(10)


\item \label{spline} \textbf{Note:} The remaining questions are all
connected, but can be worked on independently.

In the following, we consider the cubic interpolating spline on $N$
equidistant grid points $x_j = jh$ for $j=0,\dots,N-1$.

Let $y_j$ denote the given value of the spline on the grid nodes, and
$y_j''$ the (unknown) value of the second derivative.

Show that the following system of linear equations is satisfied:
\begin{equation*}
  y_{j-1}'' + 4 \, y_j'' + y_{j+1}''
  = \frac6{h^2} \, (y_{j-1} - 2 \, y_{j} + y_{j+1}) \,.
  \tag{*}
\end{equation*}

\emph{Hint:} Use the notation from class: write
\begin{equation*}
  s_j(x) = a_j \, (x-x_j)^3 + b_j \, (x-x_j)^2 + c_j \, (x-x_j) + d_j 
\end{equation*}
to denote the spline function on the $j$th interval.  Show that $2\,
b_j = y_j''$ and $d_j=y_j$ and eliminate $a_j$ and $c_j$ by using the
matching conditions at the interpolation nodes.
\hfill(10)


\item Assume that the cubic spline from question~\ref{spline} is
defined on a periodic grid.  In other words, interpolation node $x_0$
is identified with $x_N$.  

Explain why system (*) has the same number of equations as there are
unkowns, i.e., why we do not need to impose additional conditions as in
the case of non-periodic splines.
\hfill(10)


\item Let $u_j$ with $j=0,\dots,N-1$ be a given tuple of numbers on an
$N$-periodic equidistant grid with grid spacing $h=2\pi/N$.  Let $\hat
u_k$ with $k=-N/2,\dots,N/2-1$ denote its discrete Fourier transform.
Further, let $\tau_\ell$ denote translation by $\ell$ grid points,
i.e.\
\begin{equation*}
  (\tau_\ell u)_j = u_{j+\ell} \,.
\end{equation*}

Show that
\begin{equation*}
  (\widehat{\tau_\ell u})_k = \e^{\i hk\ell} \, \hat u_k \,.
\end{equation*}
\hfill(10)


\item \label{dft} The linear system (*) for computing the periodic
spline in question~\ref{spline} can be solved by taking the discrete
Fourier transform on both sides of the equality.

\begin{enumerate}
\item
Show that in the special case that the right hand side
is non-zero on a single interpolation node only, 
\begin{equation*}
  y_{j-1}'' + 4 \, y_j'' + y_{j+1}'' = \delta_{0j}
\end{equation*}
the solution of the linear system is given by
\begin{equation*}
  y_j'' = \frac1N \sum_{k=-N/2}^{N/2-1} 
          \frac{\e^{\i hkj}}{4 + 2 \, \cos(kh)} \,.
\end{equation*}

\emph{Note:} You are required to use the discrete Fourier transform
and inverse discrete Fourier transform.  The given solution is only
for your convenience---direct substitution into the linear system will
not earn credit.

\item Show that the result from part (a) implies that at the
gridpoint $j=N/2$,
\begin{equation*}
  \lvert y_{N/2}'' \rvert \leq \operatorname{const} \cdot \frac1N \,.
\end{equation*}

\emph{Hint:} First show that
\begin{equation*}
  y_{N/2}'' = \frac1N \sum_{k=1}^{N/2} \frac{(-1)^k}{2 + \cos(kh)} \,.
\end{equation*}
\end{enumerate}
\hfill(10+10)


\item Write out the linear system (*) of question~\ref{spline} in
matrix form \emph{for the periodic case}.

You now know three methods for numerically solving this system:
\begin{enumerate}
\item Gaussian elimination (or $LU$-decomposition)
\item Iterative methods 
\item The procedure outlined in question~\ref{dft}, implemented in
terms of the FFT and IFFT
\end{enumerate}
Comment on the computational complexity and efficiency of each for
this particular matrix.  
\end{enumerate}
\hfill(10)

\end{document}
