\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{May 17, 2022}
% 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}$ are 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 boundary conditions $\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 which gives
\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 derive 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 except the boundary. The boundary needs to be handled seperately. The left boundary condition is a Dirichlet constraint which is simply a constant 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 is usualy required to take the sum/difference of an Taylor expantion to the left and right direction at a given point.However, 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}
This results in an system of equations containing 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}

\begin{figure}[!h]
    \centering
    \includegraphics[width=0.48\textwidth]{plots/task01_1_00023.png}
    \includegraphics[width=0.48\textwidth]{plots/task01_x_00023.png}
    \includegraphics[width=0.48\textwidth]{plots/task01_1_00199.png}
    \includegraphics[width=0.48\textwidth]{plots/task01_x_00199.png}
    \caption{\label{fig:task01_stable}Comparison for stable (left, $d = 0.48 \leq 0.5$) vs. unstable (right, $d = 0.55 > 0.5$). At the top the simulation snapshot after 120 time steps and at the bottom after 1000 timesteps.}
\end{figure}



\subsection{Explicit scheme with Dirichlet BC at both boundaries}\label{sec:task01_2}
In this task two Dirichlet boundary conditions are given. 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 to 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}

See Figure~\ref{fig:task01_inf} for $\lim_{t\to\infty}C(x, t)$ simulation (many timesteps) for subtask 1 and 2. In comparison both lead to a straight line solution (last energy) such that the Dirichlet or Neumann BC are fullfilled.

\begin{figure}[!h]
    \centering
    \includegraphics[width=0.48\textwidth]{plots/task01_1_Cinf.png}
    \includegraphics[width=0.48\textwidth]{plots/task01_2_Cinf.png}
    \caption{\label{fig:task01_inf}Solution from Section~\ref{sec:task01_1} (left) and Section~\ref{sec:task01_2} (right) for $t\to\infty$ simulated by setting number of time steps to $N_t = 20000$.}
\end{figure}

\subsection{Implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_3}
The following computations are similar in nature to the explicit scheme. The implicit discretization for the derivatives are (see: Figure~\ref{fig:ex1})
\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}$.

\begin{figure}[!h]
    \centering
    \includegraphics[width = 0.8\textwidth]{plots/task01_all.png}
    \caption{\label{fig:task01_all}Comparison of all four subtasks. As expected 1.1, 1.3 and 1.4 solve the same probelm and as such have (basically) identical solution while the subtask 1.2 solves for different BC (Dirichlet) which is fullfilled.}
\end{figure}

\subsection{Second order in time for implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_4}
\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 one uses the Crank-Nicolson approach to derive the descritization at the time midpoints $t_{n+1/2} = (n + 1/2)\Delta t$ which are indexed with $n+1/2$.

By expantion into $\pm\Delta t / 2$ at the midpoints one 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 one need to ensure the evaluation to be done at the known grid points instead of the unknown time midpoints. Therefore, replacing the evaluation at the time midpoints by the 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 consideration of the Dirichlet and Neumann BC one 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 one defines 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}$.

For comparison against the other methods see Figure~\ref{fig:task01_all}.

\subsection{Notes on Implementation}
All four tasks follow the same basic sructure which can be written as a generalized linear system of the form
\begin{displaymath}
    L C^{n+1} = R C^n
\end{displaymath}
for updating the solution by one time step. The two matrices are dependent on the used scheme and the boundary conditions, see Table~\ref{tab:ex1}.

\begin{table}[h!]
    \centering
    \begin{tabular}{l | c c}
        Task & $L$ & $R$ \\ \hline
        1.1  & $I$ & $A_1$ \\
        1.2  & $I$ & $A_2$ \\
        1.3  & $A_3$ & $I$ \\
        1.4  & $I + A_4$ & $I - A_4$ \\
    \end{tabular}
    \caption{\label{tab:ex1}Left and right hand side matrices of the general update rule $L C^{n+1} = R C^n$ with $A_i$ depending on the subtask as defined in Section~\ref{sec:task01_1} till \ref{sec:task01_4} and $I$ the identity matrix.}
\end{table}

As a boild down solver the following is sufficient to compute the solution given an initial solution $C_0$ (argument \texttt{C}) with the appropriately defined matrices $L, R$ (arguments \texttt{L, R}) which incorporate the scheme, space and time step size $\Delta x, \Delta t$ as well as the dispersion parameter $D$ and the boundary conditions) at time $T = N_t \Delta$ (implicitly defined through the parameter \texttt{nx}).
\begin{python}
def solve(C, L, R, nt):
    for t in range(nt):
        rhs = C if R is None else R @ C
        C = rhs if L is None else solve_banded((1, 1), L.data, rhs)
    return C
\end{python}
This code assumes \texttt{C} to be a 1D \texttt{numpy} array as well as $L, R$ to by \texttt{scipy.sparse.spmatrix} sparse tridiagonal matrices or \texttt{None} which is equivalent to the identity matrix.

The additional code in the script \texttt{task01.py} contains setting up the needed matrices and plotting.

\newpage
\section{Numerical solution of the linear advection equation in a periodic domain}
First a few words about the upwind scheme. With a positive $U$ the upwind means that it is discretized by ``going left'', meaning into the opposite direction of the current position. In a similar faschion as above this leads to the discretization
\begin{align*}
    \partial_x C_i^{n} &= \frac{C_i^n - C_{i-1}^n}{\Delta x} + \mathcal{O}(\Delta x), \\
    \partial_t C_i^{n+1} &= \frac{C_i^{n+1} - C_i^n}{\Delta t} + \mathcal{O}(\Delta t).
\end{align*}
One ends up with a first order discretized equation
\begin{displaymath}
    0 = \frac{\widehat{C}_i^{n+1} - \widehat{C}_i^n}{\Delta t} + U \frac{\widehat{C}_i^n - \widehat{C}_{i-1}^n}{\Delta x}
\end{displaymath}
leading to the update rule
\begin{align*}
    C_i^{n+1} &= (1 - C_0)C_i^n + C_oC_{i-1}^n, && i = 2, ..., N_x
\end{align*}
with the current $C_0 = \frac{U \Delta x}{\Delta t}$. The boundary condition enforces periodizity which means $C_1^n = C_{N_x}^n$ and therefore the update for the left most point is
\begin{displaymath}
    C_1^{n+1} = (1 - C_0)C_1^n + C_oC_{N_x}^n.
\end{displaymath}

\begin{figure}[!h]
    \centering
    \includegraphics[width = 0.3\textwidth]{plots/task02_gauss_C0_0-7_00150.png}
    \includegraphics[width = 0.3\textwidth]{plots/task02_gauss_C0_1_00150.png}
    \includegraphics[width = 0.3\textwidth]{plots/task02_gauss_C0_1-1_00150.png}
    \includegraphics[width = 0.3\textwidth]{plots/task02_square_C0_0-7_00150.png}
    \includegraphics[width = 0.3\textwidth]{plots/task02_square_C0_1_00150.png}
    \includegraphics[width = 0.3\textwidth]{plots/task02_square_C0_1-1_00150.png}
    \caption{\label{fig:task02_all}For initial Gauss wavelet (top) and Square wavelet (bottom) after 150 iterations (one leap around the boundary) with different current $C_0$. Left: $C_0 = 0.7$, Center: $C_0 = 1$ (perfect) and Right: $C_0 = 1.1$ (unstable).}
\end{figure}

For implementation one gets a very simple solver which computes the wave form at time $T = N_t \Delta t$ with $\Delta t = U \Delta x / C_0$ as the following;
\begin{python}
def solve(C, C0, nt):
    nx = C.shape[0]
    im1 = np.array([nx - 1] + list(range(nx - 1)))
    for t in range(nt):
        C = (1 - C0) * C + C0 * C[im1]
    return C
\end{python}
where \texttt{C} is the initial wave form, the constant \texttt{C0} is the current $C_0$ and \texttt{nt} is the number of time steps to be performed. Again, the attached script \texttt{task02.py} implements this with additional code for setting up the initial conditions and plotting.

\end{document}