FEM1D.tex

\documentclass[a4paper]{article}

\author{Hidde-Jan Jongsma}
\title{Finite Element Method}
\date{\today}

% Extra symbols and environments
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{mathtools}

% Use full page
\usepackage{fullpage}

% Create graphics in LaTex
\usepackage{tikz}

% Extra commands
\newcommand{\dd}{\mathrm{d}}

\newcommand{\Reals}{\mathbb{R}}

\newcommand{\dux}{\frac{\partial u}{\partial x}}
\newcommand{\duxx}{\frac{\partial^2 u}{\partial x^2}}
\newcommand{\duxxx}{\frac{\partial^3 u}{\partial x^3}}

\newcommand{\duy}{\frac{\partial u}{\partial y}}
\newcommand{\duyy}{\frac{\partial^2 u}{\partial y^2}}
\newcommand{\duyyy}{\frac{\partial^3 u}{\partial y^3}}
\newcommand{\lapu}{\nabla^2 u}

\newcommand{\vct}{\mathbf}

\providecommand{\abs}[1]{\lvert#1\rvert}
\providecommand{\norm}[1]{\lVert#1\rVert}

\newcommand{\LO}{\ensuremath{L_2(\Omega)}}
\newcommand{\HO}{\ensuremath{H^1(\Omega)}}
\newcommand{\HOzero}{\ensuremath{H^1_{(0}(\Omega)}}


\begin{document}

Introduction

\section{Finite element method}

For the simple, one dimensional case, the exact solution to the
Helmholtz problem is known. Finding the exact solution for more
complex and higher dimensional problems turns out to be hard and
impractical.  Often, an approximation to the exact solution suffices
for engineering purposes. A widely used method to find such
approximations is known as the \emph{Galerkin finite element method}.

Using this method, an approximation to the exact solution is found by
transforming the problem in a system linear equations.  This results
in a sparse linear system for which many solving techniques have been
studied. Although the linear systems resulting from this
method are sparse, for an accurate approximation of solution to higher
dimensional or heavily oscillating problems a very large system can be
needed.

%When searching for an approximate solution to the Helmholtz problem
%using a Finite Element approach, we aim to approximate a function
%$u$ by constructing a finite linear system and solving the linear
%equations. For this approximation to be possible, we want to
%restrict our search space for $U(x) \approx u(x)$ to a finite dimensional
%function space.
%
%In this section, we will show how an approximation for $u$ can be found
%using the \emph{Galerkin finite element method} in both
%one and higher dimensions. We will give error estimates in both
%cases and we will address the problems that arise when
%dealing with heavily oscillating functions.

%\subsection{One dimensional Helmholtz problem}
\subsection{Galerkin finite element method} \label{seq:Galerkin}

As an example, we will consider a simple one dimensional wave problem.
Suppose that for $u(x)$ on $[0, 1]$ the conditions in \eqref{helmi} -
\eqref{eq:impbd} hold. 
%\begin{align}
%    - \frac{\dd^2u}{\ \dd x^2} - k^2 u
%      & = f, \quad \text{on}\ \Omega = (0, 1), \\
%    u(0) & = 0, \label{eq:dirbc} \\
%    \frac{\dd u}{\ \dd x}(1) - i k u(1) & = 0. \label{eq:impbc} &
%\end{align}
%It can be shown that the exact solution to is given by
%\begin{equation}
%  u(x) = \frac{e^{ikx}}{k} \int^x_0 \sin(ks)f(s) \ \dd s
%          + \frac{\sin(kx)}{k} \int^1_x e^{iks} f(s) \ \dd s,
%\end{equation}
The exact solution is given by \eqref{eq:1dexact}
which is periodic with wavelength $\lambda = \frac{2\pi}{k}$.
%
%While this specific problem has a simple solution, in general, higher
%dimensional problems can only be solved numerically. For these complex
%problems, no exact solution is known.
%
%\subsection{Galerkin finite element method} \label{seq:Galerkin}
Although the exact solution for this problem is known, for higher
dimensional and more complex problems no exact solution is known. For
these problems, we want to find an approximation to the exact solution
using a numerical approach.

Since the original problem is infinite dimensional, we simplify the
problem by restricting the approximation of the exact solution to a
finite dimensional search space. Suppose we have a bounded domain
$\Omega \subset \Reals^n$, $n = 1, 2, 3$.  We define the function
space $\LO$ of square integrable function on $\Omega$ by saying that
$f \in \LO$ if
\begin{equation}
  \lVert f \rVert := \left(
    \int_\Omega | f(x) |^2 \ \dd\Omega \right)^{\frac{1}{2}} < \infty.
\end{equation}
Furthermore, we say $f \in \HOzero$ if $f(0) = 0$ and $f \in \HO$,
which implies
\begin{equation}
  \lVert \nabla f \rVert^2 + \lVert f \rVert^2 < \infty.
\end{equation}

For solving the above system with a finite element method, we will
first rewrite the problem in its weak form.  We multiply both sides
with a test function $v \in \HOzero$ and integrate afterwards to
obtain
\begin{equation}
  - \int^1_0 u''(x)v(x) - k^2u(x)v(x) \ \dd x = \int^1_0 f(x)v(x) \ \dd x.
\end{equation}
By integrating the first term by parts and substituting boundary
condition \eqref{eq:impbc} and taking the fact that $v(0) = 0$ into
consideration, we obtain
\begin{equation} \label{eq:weakprob}
  \int^1_0 u'(x)v'(x)\ \dd x - k^2 \int^1_0 u(x)v(x) \ \dd x - iku(1)v(1) = \int^1_0 f(x)v(x) \ \dd x.
\end{equation}
Now we have the weak formulation of our problem; Find $u(x) \in
\HOzero$ such that \eqref{eq:weakprob} holds for all $v(x) \in
\HOzero$.  To solve this problem, first we define a finite element
mesh $X_h$ on $\Omega$, by
\begin{equation}
  X_h := \{ x_i \ | \ x_i = ih, \ i = 0, 1, \ldots, N \},
\end{equation}
where $h = 1/N$.
We limit our search space to $S_h(0,1) \subset \HOzero$, the space of
piecewise continuous linear functions with nodal values at the points
in $X_h$, satisfying \eqref{eq:dirbc}. This function space is spanned
by a set of \emph{hat functions} defined as
\begin{equation}
  \chi_j(x) = \begin{dcases*}
    \frac{1}{h} (x - x_{j - 1}), & $x \in [x_{j - 1}, x_j],$ \\
    \frac{1}{h} (x_{j + 1} - x), & $x \in [x_{j}, x_{j + 1}],$ \\
    0 & elsewhere,
  \end{dcases*}
\end{equation}
for $j = 1, 2, \ldots, N - 1$ and for $j = N$ by
\begin{equation}
  \chi_N(x) = \begin{dcases*}
    \frac{1}{h} (x - x_{j - 1}), & $x \in [x_{j - 1}, 1],$ \\
    0 & elsewhere.
  \end{dcases*}
\end{equation}
Since this set of function spans $S_h(0,1)$, they are known as
\emph{basis functions}. The intervals $\tau_i = (x_{i - 1}, x_i)$ are
called \emph{finite elements}. In our case, every $\tau_i$ has length
$h$, hence $X_h$ is called a \emph{uniform} mesh. Now, if we require
that both $U$ and $v$ are in  $S_h(0,1)$, then we can write
\begin{equation} \label{eq:Usum}
  U(x) = \sum^N_{j = 1} u_j \chi_j(x),
\end{equation}
and \eqref{eq:weakprob} transforms to
\begin{equation} \label{eq:linsys}
  \sum^N_{j = 1} \left[ \int^1_0 \chi_j'(x) \chi_m'(x) \ \dd x
    - k^2 \int^1_0 \chi_j(x) \chi_m(x) \ \dd x \right] u_j
    - i k u_N \chi_m(1)
  =
    \int^1_0 f(x) \chi_m(x) \ \dd x,
\end{equation}
for $m = 1, 2, \ldots, N$.

\subsection{Linear system}

Since each $\chi_m$ is known, we can easily compute the integrals in
\eqref{eq:linsys}. They are given by
\begin{equation*}
  \int^1_0 \chi_j'(x) \chi_m'(x) \ \dd x
    = \begin{dcases*}
        0,              & if $|j - m| > 1$,   \\
        - {1}/{h},  & if $|j - m| = 1$,   \\
        {2}/{h},    & if $j = m \neq N$,  \\
        {1}/{h},    & if $j = m = N$,
      \end{dcases*}
\end{equation*}
\begin{equation*}
  \int^1_0 \chi_j(x) \chi_m(x) \ \dd x
    = \begin{dcases*}
        0,    & if $|j - m| > 1$,   \\
        h/6,  & if $|j - m| = 1$,   \\
        2h/3, & if $j = m \neq N$,  \\
        h/3,  & if $j = m = N$.
      \end{dcases*}
\end{equation*}
Also note that
\begin{equation*}
  \chi_m(1)
  = \begin{dcases*}
      1, & if $m = N$, \\
      0  & otherwise.
    \end{dcases*}
\end{equation*}
If we substitute these equations back into \eqref{eq:linsys}, the
result is a system of linear equations with unknowns $u_i$, $i = 1, 2,
\ldots, N$. This linear system is given by
\begin{equation}
  (A - k^2 B - i k C)\vct{u} = \vct{f}.
\end{equation}
Here, the elements of matrices $A$ and $B$ are determined by
\begin{equation}
  A_{i j} = \int^1_0 \chi_i'(x) \chi_j'(x) \ \dd x, \quad
  B_{i j} = \int^1_0 \chi_i(x) \chi_j(x) \ \dd x.
\end{equation}
and
\begin{equation*}
  A = \begin{pmatrix}
        \frac{2}{h}   & - \frac{1}{h} & 0             & \hdots        & 0             \\
        - \frac{1}{h} & \frac{2}{h}   & - \frac{1}{h} &  \ddots       & \vdots        \\
        0             & - \frac{1}{h} & \frac{2}{h}   & \ddots        & 0             \\
        \vdots        & \ddots        & \ddots        & \ddots        & - \frac{1}{h} \\
        0             & \hdots        & 0             & - \frac{1}{h} & \frac{1}{h}
      \end{pmatrix},
  \quad
  B = \begin{pmatrix}
        \frac{2h}{3}  & \frac{h}{6}   & 0             & \hdots      & 0           \\
        \frac{h}{6}   & \frac{2h}{3}  & \frac{h}{6}   & \ddots      & \vdots      \\
        0             & \frac{h}{6}   & \frac{2h}{3}  & \ddots      & 0           \\
        \vdots        & \ddots        & \ddots        & \ddots      & \frac{h}{6} \\
        0             & \hdots        & 0             & \frac{h}{6} & \frac{h}{3}
      \end{pmatrix},
\end{equation*}
\begin{equation*}
  C = \begin{pmatrix}
        0       & 0       & \hdots  & 0       \\
        0       & 0       & \ddots  & \vdots  \\
        \vdots  & \ddots  & \ddots  & 0       \\
        0       & \hdots  & 0       & 1
      \end{pmatrix},
  \quad
  \vct{u} = \begin{pmatrix}
              u_1 \\
              u_2 \\
              \vdots \\
              u_N
            \end{pmatrix},
  \quad
  \vct{f} = \begin{pmatrix}
              \int^1_0 f(x) \chi_1(x) \ \dd x \\
              \int^1_0 f(x) \chi_2(x) \ \dd x \\
              \vdots \\
              \int^1_0 f(x) \chi_N(x) \ \dd x \\
            \end{pmatrix}.
\end{equation*}

We see that the system matrix $(A - k^2 B - i k C)$ is sparse, its
entries are mostly zero. Furthermore, it is also \emph{tridiagonal},
which means only its main, lower and upper diagonal are non-zero. For
this type of systems, very efficient solvers are available (see
Section \ref{seq:numfem}). After solving the above system, we plug the values
$u_i$ back into \eqref{eq:Usum} and obtain our estimate $U(x)$ to
$u(x)$.

\subsection{Error estimates}

For small $h$, the Galerkin finite element method has an error estimate
of the form
\begin{equation} \label{eq:ErrEst}
  \frac{\lVert u - U \rVert}{\lVert u \rVert} \leq C h^2,
\end{equation}
where the constant $C$ is independent of $h$ \cite[p.\ 10]{femnotes}.
Although we see that the error goes to zero as $h$ decreases (i.e. $N$
increases), it is unclear what $h$ is required for a certain level of
accuracy. Furthermore, th $C$ in \eqref{eq:ErrEst} can depend multiple
factors. For instance, it could depend on the forcing term $f$, the
exact solution $u$ or the wavenumber $k$. The last turns out to be
problematic in multiple ways.

A higher wavenumber means higher oscillatory behavior of the exact
solution. Since the wavelength $\lambda$ is defined as
$\frac{2\pi}{k}$, the wavelength decreases as $k$ increases. To keep
the number of grid points per wavelength constant and the approximation
of $u$ accurate, the step size $h$ has to decrease as $k$. A common
idea is to keep the number of grid points per wavelength constant, but
this leads to bigger systems that have to be solved and a great deal
of additional computation costs.

Another problem that accompanies higher wave numbers is the
introduction of pollution errors into the estimate $U$. Errors caused
by the wavelength not being modelled accurately accumulate and reduce
the accuracy of our estimation. It was proven \cite{femnotes} that the
relative error satisfies
\begin{equation} \label{eq:ErrEstK}
  \frac{\lVert u - U \rVert}{\lVert u \rVert} \leq C_1 k h + C_2 k^3
  h^2,
\end{equation}
where constants $C_1$, $C_2$ are independent of both $k$ and $h$. We
see that, to keep this bound constant, we have to fix $k^3 h^2$. As a
result, our linear system rapidly increases in size as $k$ increases.

\end{document}