\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.

\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 = 0$ for all time while the Neumann condition on the right requires an additional discretization step for computing the next value. Therefore, we take the second order Taylor expantion of $C$ with respect to $x$ in the negative direction given by
\begin{displaymath}
    C(x - \Delta x, t) = C(x, t) - \partial_x C(x, t)\Delta x + \mathcal{O}(\Delta x^2).
\end{displaymath}
As we are interesetd in the boundary value for the Neumann BC, the derivative at $x = h$ is given which leads to
\begin{displaymath}
    C(h, t) = C(h - \Delta x, t) + \partial_x C(h, t)\Delta x + \mathcal{O}(\Delta x^2)
\end{displaymath}
which is second order accurate. Therefore, the right boundary value for our boundary condition $\partial_x C(h, t_n) = 0$ at time $t_n$ is
\begin{displaymath}
    \widehat{C}_{N_x}^n = \widehat{C}_{N_x - 1}^n.
\end{displaymath}
See \texttt{task01\_1-2.py} or \texttt{Mixed\_EA.m} for an implementation of this scheme.

\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}
Not we have the same setting as in Section~\ref{sec:task01_1} except for Dirichlet boundary conditions on both sides given by $C(0, t) = 1$ and $C(h, t) = 0$. This means tha the left and right values are constants and internal nodes follow the same update scheme as in \eqref{eq:task01_1_update}.

\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 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 rule
\begin{displaymath}
    \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}
\end{displaymath}
for $i = 2, ..., N_x - 1$ and $d = \frac{D \Delta x^2}{\Delta t}$. The boundary conditions are ether $\widehat{C}_1^{n+1} = 0$ for the left Dirichlet and $\widehat{C}_{N_x}^{n+1} = \widehat{C}_{N_x-1}^{n+1}$ as the right Neumann BC. By collecting all coefficients of the $i = 2, ..., N_x - 1$ into a $(N_x - 2)\times (N_x - 2)$ trigiagonal matrix
\begin{displaymath}
    A = \begin{pmatrix}
        1+2d & -d &   \\
        -d & 1+2d & -d &   \\
          & -d & 1+2d & \ddots \\
          &   & \ddots & \ddots & -d \\
          & & & -d & 1+2d & -d \\
          & & &   & -d & 1+d
      \end{pmatrix}
\end{displaymath}
\todo{fix $1+d$ which is first order accurate!!!}

where the last entry $A_{N_x, N_x} = 1 - d$ is due to the Neumann BC $\widehat{C}_{N_x}^{n+1} = \widehat{C}_{N_x-1}^{n+1}$ which means
\begin{displaymath}
    \widehat{C}_i^{N_x-1} = -d \widehat{C}_i^{N_x-2} + (1 + 2 d) \widehat{C}_i^{N_x - 1} - d \widehat{C}_i^{N_x}
        = -d \widehat{C}_i^{N_x-2} + (1 + d) \widehat{C}_i^{N_x - 1}.
\end{displaymath}
Finaly, we end up with the implicit update rule
\begin{displaymath}
    \widehat{C}^{n+1} = A \widehat{C}^{n}
\end{displaymath}
which is perfomed by solving the linear system for $\widehat{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}
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$.

The second order im space is identical except the evaluation points
\begin{displaymath}
    \partial_x^2 C_i^{n+1/2} = \frac{C_{i-1}^{n+1/2} - 2C_{i}^{n+1/2} + C_{i+1}^{n+1/2}}{\Delta x^2} + \mathcal{O}(\Delta x^2)
\end{displaymath}
and for the time descritization we take the Taylor expantion


\begin{displaymath}
    \partial_{t}C_i^{n+1/2} = \frac{C_i^{n+1} - C_i^{n}}{2\Delta t} + \mathcal{O}(\Delta t^2)
\end{displaymath}
\todo{check if this follows from $\exists \partial_x^4 C$ beeing continuous}
\begin{displaymath}
    \partial_x^2 C_i^{n+1/2} = \frac{\partial_x^2 C_i^{n} + \partial_x^2 C_i^{n+1}}{2} + \mathcal{O}(\Delta x^2)
\end{displaymath}

Meaning
\begin{align*}
    0 = \partial_{t}C_i^{n+1/2}-D\partial_x^2 C_i^{n+1/2} &= \frac{C_i^{n+1} - C_i^{n}}{2\Delta t}
     -D\frac{\partial_x^2 C_i^{n} + \partial_x^2 C_i^{n+1}}{2} + \mathcal{O}(\Delta t^2 + ?) \\
     &= \frac{C_i^{n+1} - C_i^{n}}{2\Delta t} -
         D\frac{C_{i-1}^{n} + C_{i-1}^{n+1} - 2(C_{i}^{n} + C_{i}^{n+1}) + C_{i+1}^{n} + C_{i+1}^{n+1}}{2\Delta x^2}
\end{align*}
rearranging yields
\begin{displaymath}
    -d C_{i-1}^{n+1} + (1 + 2d) C_i^{n+1}-dC_{i+1}^{n+1}
        = d C_{i-1}^n + (1 - 2d)C_i^n + dC_{i+1}^n
        = -(-d C_{i-1}^n + (1 + 2d)C_i^n - dC_{i+1}^n) + 2C_i^n
\end{displaymath}
Now with second order approximation of the Neumann BC
\begin{displaymath}
    0\overset{!}{=} \partial_x C_{N_x}^{n+1} = \frac{C_{N_x-2}^{n+1} - C_{N_x}^{n+1}}{2\Delta x} + \mathcal{O}(\Delta x^2)
\end{displaymath}
which means that (\todo{check why?! should be at $N_x - 1$}) for index $i = N_x - 1$ 
\begin{align*}
    -d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}-dC_{N_x}^{n+1}
        &= -(-d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n - dC_{N_x}^n) + 2C_{N_x-1}^n \\
    -d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}-d C_{N_x-2}^{n+1}
        &= -(-d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n - d C_{N_x-2}^{n}) + 2C_{N_x-1}^n \\
    -2d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}
        &= -(-2d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n) + 2C_{N_x-1}^n
\end{align*}

\begin{displaymath}
    A C^{n+1} = (2 I - A) C^n
\end{displaymath}
with 
\begin{displaymath}
    A = \begin{pmatrix}
        1 & 0 & 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}

\end{document}