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

\section*{Theory}

\begin{enumerate}
\item Find the linear system that results from using central
  differencing to solve the elliptic equation
  \begin{equation*}
    u_{xx} + u_{yy} = \sin ( \pi y ) \left( 2 - 6 x - (\pi x)^2 (1 -
      x) \right).
  \end{equation*}
  Central differencing, trivial boundary conditions, and a $3 \times
  2$ grid should be used.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  The question is incomplete as it does not specify the spatial
  domain. However, it is reasonably to expect it to be the unit square
  $(x, y) \in [0, 1]^2$. So we first specify the grid using the
  coordinates
  \begin{align*}
    x_i & = i h_x & i & = 0, \dots, 4, & h_x & = \frac{1}{4}, \\
    y_j & = j h_y & j & = 0, \dots, 3, & h_y & = \frac{1}{3}.
  \end{align*}
  This gives an evenly spaced grid with 3 interior points in $x$ and 2
  interior points in $y$. The boundary conditions are imposed as
  \begin{align*}
    u_{0, j} & = 0, & u_{4, j} & = 0, \\
    u_{i, 0} & = 0, & u_{i, 3} & = 0.
  \end{align*}
  We are of course using the notation $u_{i,j} \approx u(x_i, y_j)$.

  We use central differencing at all interior points to get
  \begin{align*}
    u_{xx} & \approx \frac{u_{i+1,j} + u_{i-1,j} - 2 u_{i,j}}{h_x^2}, \\
    u_{yy} & \approx \frac{u_{i,j+1} + u_{i,j-1} - 2 u_{i,j}}{h_y^2}.
  \end{align*}
  We then substitute these approximations into the differential
  equation to find the difference equation satisfied at all interior
  points,
  \begin{equation*}
    \left(u_{i+1,j} + u_{i-1,j} - 2 u_{i,j} \right) + \left(
      \frac{h_x}{h_y} \right)^2 \left( u_{i,j+1} + u_{i,j-1} - 2
      u_{i,j} \right) = h_x^2 \sin ( \pi y_j ) \left( 2 - 6 x_i - (\pi
      x_i)^2 (1 - x_i) \right).
  \end{equation*}
  We will write $h_x / h_y = \alpha$.

  We choose natural ordering by rows to write the unknowns $u_{i,j}$
  as a vector ${\bf u}$ as
  \begin{equation*}
    {\bf u} = \left[ u_{1,1}, u_{1,2}, u_{2,1}, u_{2,2}, u_{3,1},
      u_{3,2} \right]^T.
  \end{equation*}
  It then follows that the matrix that defines the system has the form
  \begin{equation*}
    A =
    \begin{pmatrix}
      -2 ( 1 + \alpha ) & \alpha & 1 & 0 & 0 & 0 \\
      \alpha & -2 ( 1 + \alpha ) & 0 & 1 & 0 & 0 \\
      1 & 0 & -2 ( 1 + \alpha ) & 1 & \alpha & 0 \\
      0 & 1 & \alpha & -2 ( 1 + \alpha ) & 0 & 1 \\
      0 & 0 & 1 & 0 & -2 ( 1 + \alpha ) & \alpha \\
      0 & 0 & 0 & \alpha & 1 &  -2 ( 1 + \alpha ) 
    \end{pmatrix}
  \end{equation*}
  and the right hand side vector is solely in terms of the right hand
  side data
  \begin{equation*}
    h_x^2 \sin ( \pi y_j ) \left( 2 - 6 x_i - (\pi x_i)^2 (1 - x_i) \right)
  \end{equation*}
  as all the boundary data is trivial.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Explain the strategy for solving evolutionary PDEs using finite
  differencing methods.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  The domain is discretized using a set of points in space. These
  points will usually be evenly spaced with a (representative) step
  length $h$. The grid is extended to cover the time domain as well by
  discretizing time using a time step $\delta$. The derivatives
  defining the differential equation are replaced with finite
  difference approximations, the initial data fixes the unknowns at
  the initial time, and the boundary conditions imposed using
  equivalent (or consistent) discrete boundary data. This implies a
  finite difference formula which gives the unknown data at timestep
  $n+1$ in terms of the known data at timestep $n$ and the boundary
  data. 

  Depending on the finite difference formula derived, there may be a
  relation between the time step $\delta$ and space step $h$ for which
  the algorithm is stable and convergent, which can be checked using
  von Neumann analysis.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Write out a FTCS method for the equations
  \begin{align*}
    \partial_t y - \partial_{xx} y & = \sin(x), & x & \in [0,1], &
    y(t, 0) & = 0 = y(t, 1), \\
    \partial_t y + \partial_{x} y & = e^{-(x-1/2)^2}, & x & \in [0,1], &
    y(t, 0) & = 0 = y(t, 1).
  \end{align*}
  In both cases the initial data should be taken to be $y(0, x) =
  f(x)$. Full details of the grid, the discrete initial and boundary
  conditions, and the discrete update algorithm should be given.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  We shall use $N+1$ points in space, denoted $x_i = i h$ with $i = 0,
  \dots, N$ and the space step $h = 1 / N$; this is an evenly spaced
  grid with boundary points $x_0 = 0$ and $x_N = 1$. We will also
  discretize time using $t^n = n \delta$ where $\delta$, the timestep,
  is unconstrained. We will denote our approximate solution at $(x, t)
  = (x_i, t^n)$ as $y_i^n \approx y(x_i, t^n)$.

  The initial data will be imposed by the equivalent discrete
  condition $y_i^0 = f(x_i)$. The (trivial Dirichlet) boundary
  conditions to be imposed in both cases will be given by the discrete
  conditions
  \begin{equation*}
    y_0^n = 0, \qquad y_N^n = 0, \qquad \forall n \ge 0.
  \end{equation*}

  To form the discrete update algorithm we use forward differencing in
  time (FT),
  \begin{equation*}
    y_t|_{x=x_i, t = t^n} \approx \frac{y_i^{n+1} - y_i^{n}}{\delta},
  \end{equation*}
  and central differencing in space (CS), 
  \begin{align*}
    y_x|_{x=x_i, t = t^n} & \approx \frac{y_{i+1}^{n} - y_{i-1}^{n}}{2
      h}, \\
    y_{xx}|_{x=x_i, t = t^n} & \approx \frac{y_{i+1}^{n} + y_{i-1}^{n}
      - 2 y_i^n}{h^2}.
  \end{align*}
  We then replace the derivatives in the PDE with the finite
  difference approximations and rearrange, to get for the heat equation
  \begin{align*}
    \frac{y_i^{n+1} - y_i^{n}}{\delta} - \frac{y_{i+1}^{n} + y_{i-1}^{n}
      - 2 y_i^n}{h^2} & = \sin(x_i) \\
    \Rightarrow \qquad y_i^{n+1} = y_i^n + \delta \sin(x_i) + s \left(
      y_{i+1}^{n} + y_{i-1}^{n} - 2 y_i^n \right),
  \end{align*}
  where $s = \delta / h^2$, and for the advection equation
  \begin{align*}
    \frac{y_i^{n+1} - y_i^{n}}{\delta} + \frac{y_{i+1}^{n} -
      y_{i-1}^{n}}{2 h} & = e^{-(x_i - 1/2)^2} \\
    \Rightarrow \qquad y_i^{n+1} & = y_i^n + \delta e^{-(x_i - 1/2)^2}
    - c \left( y_{i+1}^n - y_{i-1}^n \right)
  \end{align*}
  where $c = \delta / (2 h)$. In both cases these equations hold for
  the interior points $n \ge 0$ and $i = 1, \dots, N-1$. The boundary
  points are set using the above initial and boundary data.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\item Using von Neumann analysis, find when the BTCS method is stable
  for the advection equation.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
  The finite difference formula for the BTCS method for the advection
  equation
  \begin{equation*}
    \partial_t y + v \partial_x y = 0
  \end{equation*}
  is straightforwardly checked to be
  \begin{equation*}
    y_i^{n+1} = y_i^n - c \left( y_{i+1}^{n+1} - y_{i-1}^{n+1} \right)
  \end{equation*}
  where the convection number $c = v \delta / (2 h)$ with $\delta$ the
  timestep and $h$ the space step. 

  We use the standard von Neumann analysis ansatz, which is that
  $y_{\ell}^k = q^k e^{\alpha \ell j h}$ where $q$ is the (unknown)
  growth rate to be found, $\alpha$ is related to the frequency of the
  (generic) initial data (and hence can take any real value), $h$ is
  the space step and $j = \sqrt{-1}$.

  Identifying $\ell \leftrightarrow i$ and $k \leftrightarrow n$ and
  substituting the ansatz into the finite difference equation we see
  that
  \begin{equation*}
    q = 1 - q c \left( e^{\alpha j h} - e^{- \alpha j h} \right).
  \end{equation*}
  which can be written as
  \begin{align*}
    q & = 1 - 2 j q c \sin \left(\alpha h \right) \\
    \Rightarrow \qquad q & = \frac{1}{1 + 2 j c \sin \left( \alpha h
      \right)} \\
    \Rightarrow \qquad |q|^2 & = \left| \frac{1 - 2 j c \sin \left(
          \alpha h \right)}{1 + 4 c^2 \sin^2 \left( \alpha h \right)}
    \right|^2 \\
    & = \frac{1 + 4 c^2 \sin^2 \left( \alpha h \right)}{\left[ 1 + 4
        c^2 \sin^2 \left( \alpha h \right) \right]^2} \\
    & = \frac{1}{1 + 4 c^2 \sin^2 \left( \alpha h \right)} \\
    & \le 1 \qquad \forall \alpha.
  \end{align*}
  Therefore BTCS is unconditionally stable for the advection equation.
  % 
  \begin{center}
    \rule{0.9\textwidth}{.1pt}
  \end{center}
  % 
\end{enumerate}

\end{document}

