         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         % 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 6 - Solutions}\\*[.8cm]
\end{center}

\section*{Theory}

\begin{enumerate}
\item Explain the shooting method for BVPs.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  The boundary value problem for $y(x)$ with boundary data at $x=a, b$
  is converted to an \emph{initial} value problem for $y(x)$ by, at
  first, guessing the additional (initial) boundary data $z$ at $x=a$
  that isrequired for a properly posed (i.e., completely specified)
  IVP. The IVP can then be solved using any appropriate solver to get
  some solution $y(x;z)$ that depends on the guessed initial data
  $z$. By comparing against the required boundary data at $y(b)$ we
  can check if we have the correct solution of the original BVP. To be
  precise, we can write
  \begin{equation*}
    f(z) = y(x;z)|_{x=b} - y(b),
  \end{equation*}
  a nonlinear equation for $z$. At the root where $f(z) = 0$ we have
  the appropriate initial data $z$ such that the solution of the IVP
  is also a solution of the original BVP. The root of this nonlinear
  equation can be found using any standard method such as bisection or
  the secant method.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Give a \emph{complete} algorithm for solving the BVP
  \begin{equation*}
    y'' - 3 y' + 2 y = 0, \quad y(0) = 0, \quad y(1) = 1
  \end{equation*}
  using the finite difference method. Include the description of the
  grid, the grid spacing, the treatment of the boundary conditions,
  the finite difference operators and a description of the linear
  system to be solved. You do not need to say which method would be
  used to solve the linear system, but should mention any special
  properties of the system that might make it easier to solve.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  We first choose the grid. We will use $N+2$ point to cover the
  domain $x \in [0, 1]$; this implies that we have a grid spacing $h =
  1 / (N+1)$ and we can explicitly write the coordinates of the grid
  points as
  \begin{equation*}
    x_i = h i, \qquad i = 0, 1, \dots, N+1.
  \end{equation*}
  We denote the value of the (approximate finite difference) solution
  at the grid points as $y_i (\approx y(x_i))$. We will impose the
  boundary conditions using
  \begin{align*}
    y_0 & = y(0) & y_{N+1} & = y(1) \\
    & = 0 & & = 1.
  \end{align*}
  We will use central differencing which gives
  \begin{align*}
    y'(x)|_{x = x_i} & \approx \frac{y_{i+1} - y_{i-1}}{2 h}, \\
    y''(x)|_{x = x_i} & \approx \frac{y_{i+1} + y_{i-1} - 2
      y_i}{h^2}.
  \end{align*}

  We can then substitute all of these definitions into the original
  definition to find the finite difference equation that holds for the
  interior points $i = 1, \dots, N$:
  \begin{equation*}
    y_{i+1} \left( 1 - \frac{3}{2} h \right) + y_i \left( -2 + 2 h^2 \right)
    + y_{i-1} \left( 1 + \frac{3}{2} h \right) = 0.
  \end{equation*}
  This defines a linear system for the unknowns $y_i, i = 1, \dots, N$
  of the form
  \begin{equation*}
    T {\bf y} = {\bf f}.
  \end{equation*}
  We can see that the matrix $T$ is tridiagonal and has the form
  \begin{equation*}
    T =
    \begin{pmatrix}
      -2 + 2 h^2 & 1 - \tfrac{3}{2} h & 0 & 0 & 0 & \dots & 0 \\
      1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h^2 & 0 & 0 & \dots
      & 0 \\
      0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h & 0 & \dots
      & 0 \\
      0 & 0 & \ddots & \ddots & \ddots & \dots & 0 \\
      0 & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h
      & 0 \\
      0 & \dots & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2}
      h \\
      0 & \dots & \dots & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2
    \end{pmatrix}
  \end{equation*}
  The right hand side vector results from the boundary data and is
  \begin{align*}
    {\bf f} & =
    \begin{pmatrix}
      - \left( 1 + \tfrac{3}{2} h \right) y_0 \\
      0 \\
      \vdots \\
      0 \\
      - \left( 1 - \tfrac{3}{2} h \right) y_{N+1}
    \end{pmatrix} \\
    & =
    \begin{pmatrix}
      0 \\
      0 \\
      \vdots \\
      0 \\
      - \left( 1 - \tfrac{3}{2} h \right)
    \end{pmatrix}.
  \end{align*}
  As the system is given by a tridiagonal matrix it is simple and
  cheap to solve using, e.g., the Thomas algorithm.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Explain how your algorithm would have to be modified to solve
  the BVP where the boundary condition at $x=1$ becomes the Neumann
  condition 
  \begin{equation*}
    y'(1) = 1 + \frac{e}{e-1}.
  \end{equation*}
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  First a finite difference representation of the boundary condition
  is required. A first order representation would be to use backward
  differencing
  \begin{equation*}
    \frac{y_{N+1} - y_N}{h} = 1 + \frac{e}{e-1}.
  \end{equation*}
  This can be rearranged to give
  \begin{equation*}
    y_{N+1} = y_N + h \left( 1 + \frac{e}{e-1} \right).
  \end{equation*}
  So now whenever we replaced $y(1)$ as represented by $y_{N+1}$ by
  the boundary value in the previous algorithm we must instead replace
  it with the above equation which uses the known boundary data and
  unknown interior values. 

  Explicitly, this modifies the matrix $T$ to
  \begin{equation*}
    T =
    \begin{pmatrix}
      -2 + 2 h^2 & 1 - \tfrac{3}{2} h & 0 & 0 & 0 & \dots & 0 \\
      1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h & 0 & 0 & \dots
      & 0 \\
      0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h & 0 & \dots
      & 0 \\
      0 & 0 & \ddots & \ddots & \ddots & \dots & 0 \\
      0 & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2} h
      & 0 \\
      0 & \dots & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 & 1 - \tfrac{3}{2}
      h \\
      0 & \dots & \dots & \dots & 0 & 1 + \tfrac{3}{2} h & -2 + 2 h^2 +
      \left(1 - \tfrac{3}{2} h \right)
    \end{pmatrix}
  \end{equation*}
  and the right hand side vector ${\bf f}$ to
  \begin{equation*}
    {\bf f}  =
    \begin{pmatrix}
      0 \\
      0 \\
      \vdots \\
      0 \\
      - \left( 1 - \tfrac{3}{2} h \right) h \left( 1 + \frac{e}{e-1}
        \right)
    \end{pmatrix}.
  \end{equation*}
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\end{enumerate}

\end{document}

