         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % Information sheet for the Matlab lab - Maths 6111 %
         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[10pt]{article} 
\input ma_no_html_header

\usepackage{color}
\usepackage{hyperref}
\hypersetup{breaklinks=true,colorlinks=true}
%\input ma_header
\setlength{\parindent}{0pt}
\pagestyle{myheadings}
% \markright{
% \protect {\protect \epsfxsize=0.2 true cm \protect \epsffile {dolph.line.eps}}
% \it Maths 3018/6111 - Numerical methods \hfill}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\thispagestyle{empty}
\begin{center}
\textbf{\Large Maths 3018/6111 - Numerical Methods \\*[8mm]
Worksheet 2 - Solutions}\\*[.8cm]
\end{center}

\section*{Theory}

\begin{enumerate}
\item Perform the $LU$ decomposition of
  \begin{equation*}
    A =
    \begin{pmatrix}
      1 & 3 \\ 4 & 16
    \end{pmatrix}.
  \end{equation*}
  Use both standard factorisation methods.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  Remember the general algorithm from Table 2.1 in the notes.

  We first assume the Doolittle factorization ($\ell_{kk} = 1$). We
  immediately get that
  \begin{align*}
    u_{11} & = \left( a_{11} - \sum_{s=1}^{0} \ell_{1s} u_{s1} \right) /
    \ell_{11} \qquad \text{(using $k=1$)} \\
    & = 1.
  \end{align*}
  We then build the first row of $U$ to get
  \begin{align*}
    u_{12} & = \left( a_{12} - \sum_{s=1}^{0} \ell_{1s} u_{s2} \right) /
    \ell_{11} \qquad \text{(using $k=1$)} \\
    & = 3.
  \end{align*}
  We then build the first column of $L$ to get
  \begin{align*}
    \ell_{21} & = \left( a_{21} - \sum_{s=1}^{0} \ell_{2s} u_{s1} \right) /
    \ell_{11} \qquad \text{(using $k=1$)} \\
    & = 4.
  \end{align*}
  Given the factorization we have that $\ell_{22} = 1$ and hence that
  \begin{align*}
    u_{22} & = \left( a_{22} - \sum_{s=1}^{1} \ell_{2s} u_{s2} \right) /
    \ell_{22} \qquad \text{(using $k=2$)} \\
    & = (16 - 4 \times 3) \\
    & = 4.
  \end{align*}
  Hence we have
  \begin{equation*}
    L =
    \begin{pmatrix}
      1 & 0 \\ 4 & 1
    \end{pmatrix}, \qquad U =
    \begin{pmatrix}
      1 & 3 \\ 0 & 4
    \end{pmatrix}.
  \end{equation*}
  Quickly check that
  \begin{align*}
    L U & = \begin{pmatrix}
      1 & 0 \\ 4 & 1
    \end{pmatrix} 
    \begin{pmatrix}
      1 & 3 \\ 0 & 4
    \end{pmatrix} \\
    & = \begin{pmatrix}
      1 & 3 \\ 4 & 16
    \end{pmatrix} \\
    & = A.
  \end{align*}

  With the Crout factorization ($u_{kk} = 1$), we immediately get
  that
  \begin{align*}
    \ell_{11} & = \left( a_{11} - \sum_{s=1}^{0} \ell_{1s} u_{s1} \right) /
    u_{11} \qquad \text{(using $k=1$)} \\
    & = 1.
  \end{align*}
  We then build the first row of $U$ to get
  \begin{align*}
    u_{12} & = \left( a_{12} - \sum_{s=1}^{0} \ell_{1s} u_{s2} \right) /
    \ell_{11} \qquad \text{(using $k=1$)} \\
    & = 3.
  \end{align*}
  We then build the first column of $L$ to get
  \begin{align*}
    \ell_{21} & = \left( a_{21} - \sum_{s=1}^{0} \ell_{2s} u_{s1} \right) /
    \ell_{11} \qquad \text{(using $k=1$)} \\
    & = 4.
  \end{align*}
  Given the factorization we have that $u_{22} = 1$ and hence that
  \begin{align*}
    \ell_{22} & = \left( a_{22} - \sum_{s=1}^{1} \ell_{2s} u_{s2} \right) /
    u_{22} \qquad \text{(using $k=2$)} \\
    & = (16 - 4 \times 3) \\
    & = 4.
  \end{align*}
  Hence we have
  \begin{equation*}
    L =
    \begin{pmatrix}
      1 & 0 \\ 4 & 4
    \end{pmatrix}, \qquad U =
    \begin{pmatrix}
      1 & 3 \\ 0 & 1
    \end{pmatrix}.
  \end{equation*}
  Quickly check that
  \begin{align*}
    L U & = \begin{pmatrix}
      1 & 0 \\ 4 & 4
    \end{pmatrix} 
    \begin{pmatrix}
      1 & 3 \\ 0 & 1
    \end{pmatrix} \\
    & = \begin{pmatrix}
      1 & 3 \\ 4 & 16
    \end{pmatrix} \\
    & = A.
  \end{align*}

  Note that in this case only one entry has changed; in general the
  different factorizations will give totally different answers.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item{} [Additional] Write out the Thomas algorithm for a tridiagonal
  system.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  We assume that the tridiagonal matrix $A$ has diagonal entries $b_i,
  i = 1, \dots, N$, sub-diagonal entries $a_j, j = 1, \dots, N-1$, and
  super-diagonal entries $c_j, j = 1, \dots, N-1$. We assume that we
  are trying to solve the problem
  \begin{equation*}
    A {\bf x} = {\bf f},
  \end{equation*}
  which has size $N$.

  The Thomas algorithm is essentially Gaussian elimination; forward
  elimination is used to get an upper triangular problem, and then
  back substitution used to find the answer. The tridiagonal nature of
  the problem means that it is sufficiently simple to write it out
  explicitly. 

  The first step is the forward elimination process. This yields the
  system 
  \begin{equation*}
    B {\bf x} = {\bf y},
  \end{equation*}
  where $B$ has entries only on the diagonal, and these entries are
  all 1, and on the super-diagonal, which are written as $c_j / d_j, j
  = 1, \dots, N-1$. The vectors ${\bf d}, {\bf y}$ are given by the
  procedure 
  \begin{enumerate}
    %
  \item At the first step $d_1=b_1$ and $y_1=f_1/d_1$;
    % 
  \item At the $k$-th step $d_k=b_k-a_{k-1}c_{k-1}/d_{k-1}$ and
    $y_k=(f_k-y_{k-1} a_{k-1})/d_k$.
  \end{enumerate}
  We can then solve this simplified system using back substitution to find
  \begin{enumerate}
  \item At the first step $x_N = y_N$;
  \item In reverse order, $k = N-1, \dots, 1$ we have
    $x_{k}=y_{k}-x_{k+1} c_{k}/d_{k}$.
  \end{enumerate}
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Write down the general framework for iterative methods for
  linear systems. Give the convergence matrix. If the linear system
  uses the matrix $A$ above, will an iterative method converge? [Hint:
  remember what to do with the diagonal entries]
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  %
  The system to be solved is
  \begin{equation*}
    A {\bf x} = {\bf b},
  \end{equation*}
  with $\det(A) \ne 0$. \emph{We first scale the problem such that all
    diagonal entries of $A$ are 1}. We split the coefficient matrix
  $A$ into the matrices $N, P$, such that
  \begin{equation*}
    A = N - P
  \end{equation*}
  where $\det(N) \ne 0$. We therefore have that
  \begin{equation*}
    N {\bf x} = P {\bf x} + {\bf b}
  \end{equation*}
  and we use this to write an iteration scheme as
  \begin{equation*}
    N {\bf x}^{(n)} = P {\bf x}^{(n-1)} + {\bf b}, \qquad n = 1, 2, 3,
    \dots
  \end{equation*}
  where we start from some arbitrary initial guess ${\bf x}^{(0)}$.

  The convergence matrix $M$ is given by $M = N^{-1} P$. Convergence
  is guaranteed if the spectral radius of $M$ is less than one. For
  the matrix given in the first part we have the \emph{rescaled} $A$
  is
  \begin{equation*}
    A =
    \begin{pmatrix}
      1 & 3 \\ 1/4 & 1
    \end{pmatrix}.
  \end{equation*}
  Clearly to talk about convergence we need to first decide how to
  split $A$ into $N, P$ matrices. For simplicity we consider Jacobi's
  method for which $N = I$ and hence we have
  \begin{align*}
    N & =
    \begin{pmatrix}
      1 & 0 \\ 0 & 1
    \end{pmatrix}, \\
    P & =
    \begin{pmatrix}
      0 & -3 \\ -1/4 & 0
    \end{pmatrix}, \\
    M & = N^{-1} P \\
      & = P \\
      & =    
    \begin{pmatrix}
      0 & -3 \\ -1/4 & 0
    \end{pmatrix}, \\
    \Rightarrow \qquad \rho(M) & = \max_i | \lambda_i | \\
    & = \max \left\{ \pm \sqrt{-3/4} \right\} \\
    & = \sqrt{3/4} \\
    & < 1.
  \end{align*}
  Therefore Jacobi will converge for this matrix.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Check which of the matrices on this sheet are diagonally
  dominant.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  For (strict) diagonal dominance we need the absolute value of the
  diagonal coefficient to be (strictly) greater than the sum of the
  absolute values of all other coefficients in the row.

  For $A$ this is not true; for the first row we have $1 \ngtr 3$,
  although dominance holds for the second row.

  For $B$ it is not true; for the second row we have $28 \ngtr 24 +
  53$, although dominance holds for the other rows.

  For $C$ it is not true; for the first row it fails, although
  dominance holds for the other rows.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Briefly explain what is meant by quadrature methods based on
  polynomial interpolation.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  {\it Standard exam question; see, e.g., 07/08}.

  The function $f(x)$ to be integrated is replaced by an interpolating
  function, in this case a polynomial of degree $n$, that interpolates
  it at $n+1$ \emph{nodes} $x_j$, $j=0,1,\ldots,n$. In general, when
  constructing a compound quadrature formula, we may split the
  interval into subintervals, each subinterval containing $n+1$ nodes,
  and use a polynomial interpolating function on each subinterval. The
  integral of the function is approximated by the integral of the
  (piecewise polynomial) interpolating function.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item{} [3018 only] Write down the contraction mapping theorem. Check
  that $g(x) = \cos(x)$ is contracting on the unit interval.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  The contraction mapping theorem is
  \begin{quote}
    If $g(x)$ is a contraction mapping in an interval $I=[a,b]$ then
    there exists one and only one fixed point of the map in $[a,b]$.
  \end{quote}
  In order to state we need the definition of a contraction mapping,
  which is
  \begin{quotation}
    \textbf{Definition} - A continuous map $g(x)$ from an interval
    $I=[a,b] \subseteq \bR$ into $\bR$ is \textit{contracting} (or a
    \emph{contraction mapping}) if
    % 
    \begin{enumerate}
      % 
    \item the image of $I$ is contained in $I$:
      % 
      \begin{align*}
        g(I) & \subseteq I \\
        \Leftrightarrow \quad g(x) & \in I \, \forall x \in I .
      \end{align*}
      % 
    \item the function $g(x)$ is Lipschitz continuous in $I$ with
      Lipschitz constant $L < 1$:
      % 
      \begin{equation*} 
        | g(x) - g(y) | \le L | x - y | \quad \forall x,y \in I .
      \end{equation*}
    \end{enumerate}
  \end{quotation}

  Given this we want to check $g(x) = \cos(x)$ on $I = [0, 1]$. We
  have that $g$ is continuous and differentiable, so it must be
  Lipschitz continuous. We know that $L < \max_{x \in I} |g'(x)|$ so
  we look at $g' = -\sin(x)$. We note that the extrema of $g'$ occur
  when $g'' = - \cos(x) = 0$ which is at $x = (n+1/2)\pi$, of which
  none occur in the interval $[0, 1]$. We also note that $g'$ is
  monotonic on the interval. Therefore the extreme values of $g'$ are
  taken at the ends of the interval and are $[0, -\sin(1)]$, which are
  both less than 1 in absolute value. So $L<1$.

  We also need to check that $g(I) \subseteq I$. Again we note that
  $g$ is monotonic on the interval. Hence we only need to check that
  $g(0) = 1 \in I$ and $g(1) = \cos(1) \in I$, both of which
  hold. Therefore the map is contracting.
  %
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\end{enumerate}

\end{document}

