NSSC/Exercise_03/report.tex

\documentclass[a4paper, 10pt]{article}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{fullpage}
\usepackage{amsmath, amssymb, amstext, amsthm}
\usepackage[pdftex]{hyperref}
\usepackage{xcolor, graphicx}
\usepackage{tikz}
\usepackage{listings}

% Document meta into
\title{NSSC 2 - Assignement 3}
\author{Bianchi Riccardo, Kapla Daniel, Kuen Jakob, Müller David}
\date{November 24, 2021}
% Set PDF title, author and creator.
\AtBeginDocument{
    \hypersetup{
        pdftitle = {NSSC 2 - Exercise 3},
        pdfauthor = {Group 1},
        pdfcreator = {\pdftexbanner}
    }
}


\makeindex

% Setup environments
% Theorem, Lemma
\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{example}{Example}
% Definition
\theoremstyle{definition}
\newtheorem{defn}{Definition}
% Remark
\theoremstyle{remark}
\newtheorem{remark}{Remark}

\DeclareMathOperator*{\argmin}{{arg\,min}}
\DeclareMathOperator*{\argmax}{{arg\,max}}
\renewcommand{\t}[1]{{#1^T}}
\newcommand{\todo}[1]{{\color{red}TODO: #1}}


% Default fixed font does not support bold face
\DeclareFixedFont{\ttb}{T1}{txtt}{bx}{n}{10} % for bold
\DeclareFixedFont{\ttm}{T1}{txtt}{m}{n}{10}  % for normal

% % Custom colors
\definecolor{deepblue}{rgb}{0,0,0.5}
\definecolor{deepred}{rgb}{0.6,0,0}
\definecolor{deepgreen}{rgb}{0,0.5,0}

% Python style for highlighting
\newcommand\pythonstyle{\lstset{
    language=Python,
    basicstyle=\ttm,
    morekeywords={self},              % additional keywords
    keywordstyle=\ttb\color{deepblue},
    stringstyle=\color{deepgreen},
    frame=,                           % t, b, or tb for top, bottom env. lines
    showstringspaces=false,
    numbers=left,
    stepnumber=1,
    numbersep=6pt,
    numberstyle=\color{gray}\tiny
}}

% Python environment
\lstnewenvironment{python}[1][] {
    \pythonstyle
    \lstset{#1}
}{}

\begin{document}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                Exercise 1                                %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical solution of the diffusion equation in a finite domain}
We are given the 1D unsteady diffusion equation
\begin{equation}\label{eq:ex1}
    \frac{\partial C}{\partial t} - D \frac{\partial^2 C}{\partial x^2} = 0
\end{equation}
with a scalar diffusion coefficient $D = 10^{-6}$. The domain size is $h$ which is discretized with $N_x$ points to give the grid space $\Delta x$ with a time step $\Delta t$. The initial condition is $C = 0$ inside the domain.

Furthermore, denote with $N_t\in\mathbb{N}$ the number of time steps of the simulation where $t_0 = 0$ is the initial time. All $N_t$ discretized time points $t_n = n \Delta t$ leading to the final time point of the simulation as $t_{N_t} = N_t \Delta t$. The discretized $x$ points of the domain are $x_i = (i - 1)\Delta x = (i - 1)h / (N_x - 1)$ for $i = 1, ..., N_x$. This gives exactly $N_x$ grid points with equal distance between adjacent points on the domain $[0, h]$ with the first point $x_1 = 0$ and the last point $x_{N_x} = h$.

In the following we will also use the short hand notation
\begin{displaymath}
    C_i^n = C(x_i, t_n), \qquad
    \partial_x^k C_i^n = \left.\frac{\partial^k C(x, t)}{\partial x^k}\right|_{x = x_i, t = t_n}, \qquad
    \partial_t^k C_i^n = \left.\frac{\partial^k C(x, t)}{\partial t^k}\right|_{x = x_i, t = t_n}
\end{displaymath}
for $k\in\mathbb{N}$ as well as $i$ and $n$ are the space and time discretization indices, respectively. To distringuish between the exact solution and the approximation we add a ``hat'' $\widehat{\phantom{C}}$ to the approximation. For example $C_i^n$ is the true value and $\widehat{C}_i^n$ is the discretized solution with discretization errors at the evaluation point $(x_i, t_n)$.

Furthermore, define
\begin{displaymath}
    d = \frac{D\Delta t}{\Delta x^2}.
\end{displaymath}

\subsection{Explicit scheme with Dirichlet/Neumann BC}\label{sec:task01_1}
Given the Dirichlet boundary conditions $C(0, t) = 0$ and the Neumann BC $\partial_x C(h, t) = 0$. First we derive a $2^{nd}$ order \emph{explicit} finite difference approach with a $1^{st}$ order discretization in time (see Figure~\ref{fig:ex1}).

\begin{figure}[h!]
    \centering
    \begin{tikzpicture}[>=latex]
        \begin{scope}
            \draw[->] (0.8, 0) -- (5.5, 0) node[anchor = west] {$x$};
            \draw[->] (1, -0.2) -- (1, 3.5) node[anchor = south] {$t$};
            \draw[dashed] (0.8, 2) node[anchor = east] {$t_n$} -- (5.5, 2);
            \draw[dashed] (3, -0.2) node[anchor = north] {$x_i$}
                -- (3, 3.5) node[anchor = south] {explicit};
            \foreach \x in {1, ..., 5} {
                \foreach \y in {0, ..., 3} {
                    \node[circle, draw, fill = white,
                        inner sep = 0pt, outer sep = 0pt,
                        minimum size = 4pt] at (\x, \y) {};
                }
            }
            \node[circle, fill = black,
                inner sep = 0pt, outer sep = 0pt,
                minimum size = 4pt] at (3, 2) {};
            \foreach \x/\y in {2/2, 4/2, 3/3} {
                \node[circle, fill = gray,
                    inner sep = 0pt, outer sep = 0pt,
                    minimum size = 4pt] at (\x, \y) {};
            }
        \end{scope}
        \begin{scope}[xshift = 7cm]
            \draw[->] (0.8, 0) -- (5.5, 0) node[anchor = west] {$x$};
            \draw[->] (1, -0.2) -- (1, 3.5) node[anchor = south] {$t$};
            \draw[dashed] (0.8, 2) node[anchor = east] {$t_n$} -- (5.5, 2);
            \draw[dashed] (3, -0.2) node[anchor = north] {$x_i$}
                -- (3, 3.5) node[anchor = south] {implicit};
            \foreach \x in {1, ..., 5} {
                \foreach \y in {0, ..., 3} {
                    \node[circle, draw, fill = white,
                        inner sep = 0pt, outer sep = 0pt,
                        minimum size = 4pt] at (\x, \y) {};
                }
            }
            \node[circle, fill = black,
                inner sep = 0pt, outer sep = 0pt,
                minimum size = 4pt] at (3, 2) {};
            \foreach \x/\y in {2/3, 4/3, 3/3} {
                \node[circle, fill = gray,
                    inner sep = 0pt, outer sep = 0pt,
                    minimum size = 4pt] at (\x, \y) {};
            }
        \end{scope}
    \end{tikzpicture}
    \caption{\label{fig:ex1}Dependency relation for the \emph{explicit} (left) and \emph{implicit} (right) finite difference approximation at a grid point indexed $(i, n)$.}
\end{figure}

To derive a second order scheme (assuming $C$ is three time continuously differentiable as a function from $[0, h]\times\mathbb{R}^+\to\mathbb{R}$) in space we first considure the Taylor expantion of $C$ with respect to $x$ given by
\begin{displaymath}
    C(x + \Delta x, t) = C(x, t) + \partial_x C(x, t)\Delta x + \partial_x^2 C(x, t)\frac{\Delta x^2}{2} + \partial_x^3 C(x, t)\frac{\Delta x^3}{6} + \mathcal{O}(\Delta x^4).
\end{displaymath}
Replace $\Delta x$ with $\pm\Delta x$ and add the two equation together, then
\begin{displaymath}
    C(x + \Delta x, t) + C(x - \Delta x, t) = 2 C(x, t) + \partial_x^2 C(x, t)\Delta x^2 + \mathcal{O}(\Delta x^4).
\end{displaymath}
The first and third order terms drop out due to different signs. Finally, this results in
\begin{displaymath}
    \partial_x^2 C(x, t) = \frac{C(x + \Delta x, t) - 2 C(x, t) + C(x - \Delta x, t)}{\Delta x^2} + \mathcal{O}(\Delta x^2)
\end{displaymath}
which is an approximation of the second derivative with error proportional to $\Delta x^2$. For the first order approximation in time the same trick can be used except that only the first two terms of the Taylor expantion need to be considured.

Therefore, the $2^{nd}$ order scheme for space and $1^{st}$ order in time at $(i, n)$ using the explicit scheme derives as
\begin{displaymath}
    \partial_x^2 \widehat{C}_i^n = \frac{\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}}{\Delta x^2},
    \qquad
    \partial_t \widehat{C}_i^n = \frac{\widehat{C}_{i}^{n+1} - \widehat{C}_{i}^{n}}{\Delta t}.
\end{displaymath}
Substitution into \eqref{eq:ex1} yields after rearranging the update rule for internal points as
\begin{equation}\label{eq:task01_1_update}
    \widehat{C}_i^{n+1} = \widehat{C}_i^n + \frac{D\Delta t}{\Delta x^2}(\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}).
\end{equation}
This holds at $x_i$ where $i = 2, ..., N_x - 1$, or in other words everywhere inside the space domain excluding the boundary. The boundary needs to be handled seperately. The left boundary condition is a Dirichlet constraint which is simply constanct giving $C_1^n = 1$ for all time while the Neumann condition on the right requires an additional discretization step for computing the next value. To achieve a $2^{nd}$ order discretization scheme it's usualy required to take the sum/difference of an Taylor expantion into the left and right direction at a given point. As we are on the boundary there is no point at the right, which can be simply solved by adding a ghost point $x_{N_x+1}$ to the system (which will not show up in the final scheme). This leads to the two $3^{rd}$ order expantions evaluated at the right boundary as
\begin{align*}
    C_{N_x-1}^n &= C_{N_x}^n - \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3), \\
    C_{N_x+1}^n &= C_{N_x}^n + \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3).
\end{align*}
Taking the difference of both equations gives
\begin{displaymath}
    0 = \partial_x C_{N_x}^n\Delta x = \frac{C_{N_x-1}^n - C_{N_x+1}^n}{2\Delta x} + \mathcal{O}(\Delta x^2).
\end{displaymath}
which yields a $2^{nd}$ order equation for the Neumann BC on the right as
\begin{displaymath}
    C_{N_x-1}^n = C_{N_x+1}^n.
\end{displaymath}
Now where ready to give an equation system of the entire system (including the boundaryies) which reads
\begin{align*}
    \widehat{C}_1^{n+1} &= 1, \\
    \widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n},
        && i = 2, ..., N_x - 1 \\
    \widehat{C}_{N_x}^{n+1} &=  2d\widehat{C}_{N_x-1}^{n} + (1 - 2d)\widehat{C}_{N_x}^{n}.
\end{align*}
or in matrix form
\begin{displaymath}
    C^{n+1} = A_1 C^{n}
\end{displaymath}
for a tridiagonal matrix $A_1$ of dimensions $N_x\times N_x$
\begin{displaymath}
    A_1 = \begin{pmatrix}
        1 & 0 \\
        d & 1-2d & d \\
        & d & 1-2d & d \\
        & & d & 1-2d & \ddots \\
        & & &  \ddots & \ddots & d \\
        & & & & d & 1-2d & d \\
        & & & & & 2d & 1-2d
    \end{pmatrix}.
\end{displaymath}

\subsubsection{Stability} \todo{Prove that the update scheme is unstable for $d = \frac{D\Delta x^2}{\Delta t} > 0.5$}

\subsection{Explicit scheme with Dirichlet BC at both boundaries}\label{sec:task01_2}
This is identical to the previous task except the right boundary condition. The new boundary condition is a Dirichlet BC with $C_{N_x}^n = 0$ which alters the equation system from Section~\ref{sec:task01_1} in the last two equations to
\begin{align*}
    \widehat{C}_1^{n+1} &= 1, \\
    \widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n},
        && i = 2, ..., N_x - 2 \\
    \widehat{C}_{N_x-1}^{n+1} &= d\widehat{C}_{N_x-2}^{n} + (1 - 2d)\widehat{C}_{N_x-1}^{n}, \\
    \widehat{C}_{N_x}^{n+1} &= 0.
\end{align*}
which yields a new tridiagonal matrix $A_2$ of size $N_x\times N_x$ for the update rule $C^{n+1} = A_2 C^{n}$ as
\begin{displaymath}
    A_2 = \begin{pmatrix}
        1 & 0 \\
        d & 1-2d & d \\
        & d & 1-2d & d \\
        & & d & 1-2d & \ddots \\
        & & &  \ddots & \ddots & d \\
        & & & & d & 1-2d & 0 \\
        & & & & & 0 & 0
    \end{pmatrix}.
\end{displaymath}

\todo{include pictures and compare with Section~\ref{sec:task01_1}}

\subsection{Implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_3}
The following computations are similar in nature to the explicit scheme, therefore well keep it short. The implicit discretization for the derivatives is (see: Figure~\ref{fig:ex1}) then
\begin{displaymath}
    \partial_x^2 \widehat{C}_{i}^{n+1} = \frac{\widehat{C}_{i-1}^{n+1} - 2\widehat{C}_{i}^{n+1} + \widehat{C}_{i+1}^{n+1}}{\Delta x^2},
    \qquad
    \partial_t \widehat{C}_i^n = \frac{\widehat{C}_{i}^{n+1} - \widehat{C}_{i}^{n}}{\Delta t}.
\end{displaymath}
Substitution into \eqref{eq:ex1} gives the implicit (inverse) update rules
\begin{align*}
    \widehat{C}_1^{n} &= \widehat{C}_1^{n+1} = 1, \\
    \widehat{C}_i^{n} &= -d \widehat{C}_{i-1}^{n+1} + (1 + 2 d) \widehat{C}_i^{n+1} - d \widehat{C}_{i+1}^{n+1} && i = 2, ..., N_x - 1 \\
    \widehat{C}_{N_x}^{n} &= -2d \widehat{C}_{N_x-1}^{n+1} + (1 + 2 d) \widehat{C}_{N_x}^{n+1}
\end{align*}
which leads again to an matrix equation
\begin{displaymath}
    C^n = A_3 C^{n+1}
\end{displaymath}
where the $N_x\times N_x$ tridiagonal matrix $A_3$ has the form
\begin{displaymath}
    A_3 = \begin{pmatrix}
        1 & 0 \\
        -d & 1+2d & -d \\
        & -d & 1+2d & -d \\
        & & -d & 1+2d & \ddots \\
        & & &  \ddots & \ddots & -d \\
        & & & & -d & 1+2d & -d \\
        & & & & & -2d & 1+2d
      \end{pmatrix}
\end{displaymath}
and the update is performed by solving for $C^{n+1}$.

\todo{include pictures and compare with Section~\ref{sec:task01_1}}

\subsection{Second order in time for implicit scheme with Dirichlet/Neumann BC}
\begin{figure}[h!]
    \centering
    \begin{tikzpicture}[>=latex]
        \draw[->] (0.8, 0) -- (5.5, 0) node[anchor = west] {$x$};
        \draw[->] (1, -0.2) -- (1, 3.5) node[anchor = south] {$t$};
        \draw[dashed] (0.8, 1.5) node[anchor = east] {$t_{n+1/2}$} -- (5.5, 1.5);
        \draw[dashed] (3, -0.2) node[anchor = north] {$x_i$}
            -- (3, 3.5) node[anchor = south] {Crank-Nicolson};
        \foreach \x in {1, ..., 5} {
            \foreach \y in {0, ..., 3} {
                \node[circle, draw, fill = white,
                    inner sep = 0pt, outer sep = 0pt,
                    minimum size = 4pt] at (\x, \y) {};
            }
        }
        \node[circle, fill = black,
            inner sep = 0pt, outer sep = 0pt,
            minimum size = 4pt] at (3, 1.5) {};
        \foreach \x/\y in {2/1, 3/1, 4/1, 2/2, 3/2, 4/2} {
            \node[circle, fill = gray,
                inner sep = 0pt, outer sep = 0pt,
                minimum size = 4pt] at (\x, \y) {};
        }
    \end{tikzpicture}
    \caption{\label{fig:ex4}.Crank-Nicolson dependency relations.}
\end{figure}

To derive the second order accurate scheme in time we employ the Crank-Nicolson approach to derive the descritization at the time midpoints $t_{n+1/2} = (n + 1/2)\Delta t$ which we index with $n+1/2$.

By expantion into $\pm\Delta t / 2$ at the midpoints we get
\begin{align*}
    C_i^n &= C_i^{n+1/2} - \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3), \\
    C_i^{n+1} &= C_i^{n+1/2} + \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3).
\end{align*}
Taking the difference yields a $2^{nd}$ order approx at the time midpoints
\begin{displaymath}
    \partial_t C_{i}^{n+1/2} = \frac{C_i^{n+1} - C_i^n}{\Delta t} + \mathcal{O}(\Delta t^2).
\end{displaymath}
while the $2^{nd}$ order space discretization is analog to above but evaluated also at $(x_i, t_{n+1/2})$ to get
\begin{displaymath}
    \partial_x^2 C_{i}^{n+1/2} = \frac{C_{i-1}^{n+1} - 2C_{i}^{n+1} + C_{i+1}^{n+1}}{\Delta x^2} + \mathcal{O}(\Delta x^2)
\end{displaymath}
leading to an second order in time and space discretization scheme
\begin{displaymath}
    \partial_t \widehat{C}_{i}^{n+1/2} - D \partial_x^2 \widehat{C}_{i}^{n+1/2} = \mathcal{O}(\Delta x^2 + \Delta t^2).
\end{displaymath}

Next we need to ensure the evaluation to be done at the known grid points instead of the unknown time midpoints. Therefore, we replace the evaluation at the time midpoints by there in time mean of the two direct neighbouring grid points
\begin{displaymath}
    \widehat{C}_i^{n+1/2} = \frac{\widehat{C}_i^{n} + \widehat{C}_i^{n+1}}{2}.
\end{displaymath}

After a bit of algebra and condideration of the Dirichlet and Neumann BC we end up with
\begin{align*}
    \widehat{C}_1^{n+1} &= \widehat{C}_1^{n} = 1, \\
    \widehat{C}_i^{n+1} - \frac{d}{2}\widehat{C}_{i-1}^{n+1} + dC_{i}^{n+1} - \frac{d}{2}\widehat{C}_{i+1}^{n+1}
        &= C_i^{n} + \frac{d}{2}\widehat{C}_{i-1}^{n} - dC_{i}^{n} + \frac{d}{2}\widehat{C}_{i+1}^{n},
            && i = 2, ..., N_x - 1, \\
    \widehat{C}_i^{n+1} - d\widehat{C}_{N_x-1}^{n+1} + dC_{N_x}^{n+1}
            &= C_i^{n} + d\widehat{C}_{N_x-1}^{n} - dC_{N_x}^{n}.
\end{align*}
To write this again in a compact matrix equation system we define the $N_x\times N_x$ matrix $A_4$ as
\begin{displaymath}
    A_4 = \begin{pmatrix}
        0 & 0 \\
        -d/2 & d & -d/2 \\
        & -d/2 & d & \ddots \\
        & & \ddots & \ddots & -d/2 \\
        & & & -d/2 & d & -d/2 \\
        & & & & -d & d\\
    \end{pmatrix}
\end{displaymath}
then the entire system has the form
\begin{displaymath}
    (I + A_4)C^{n+1} = (I - A_4)C^n
\end{displaymath}
where $I$ is the identity. Again, the update is performed by solving this tridiagonal system for $C^{n+1}$.

\end{document}