{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TD5: Solving the advection equation with finite-difference methods\n",
    "\n",
    "The goal of this TD is to solve the following advection equation with finite-difference methods:\n",
    "$$\n",
    "    \\begin{cases}\n",
    "        \\displaystyle{\\frac{\\partial u}{\\partial t} + V \\frac{\\partial u}{\\partial x} = 0} \\text{ with } (t,x)\\in\\mathbb{R}_*^+\\times]-L, L[ \\\\\n",
    "        u (t, x + 2L) = u (t, x) \\text { with } (t,x) \\in \\mathbb{R}_*^+\\times[-L, L] \\\\\n",
    "        u (0, x) = u_0 (x) \\text{ with } x \\in [-L, L]\n",
    "   \\end{cases}\n",
    "$$\n",
    "Given an initial function $u_0$, a uniform velocity field $V$ and periodic boundary conditions on the 1D spatial domain $[-L,L]$, we want to solve the advection equation using **three different numerical schemes**: an explicit centered scheme, the Lax-Friedrichs scheme and the Lax-Wendroff scheme. We will focus on the **stability analysis** of these schemes and on some **unwanted numerical behaviors** that stem from them.\n",
    "\n",
    "We will note $\\Delta t$ and $\\Delta x$ the time step and the spatial step respectively. For a given numerical solution $u (t,x)$, $u_j^n$ will refer to the computed value $u (t_n, x_j)$ where $(t_n, x_j) = (n\\Delta t, -L + j\\Delta x)$. With $N = 2L / \\Delta x$, periodic boundary conditions read: $\\forall n, u_{j+N}^n = u_j^n$. In practice we will be interested in computing the solution at a given time $T$.\n",
    "\n",
    "We will need in the following the usual Python packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To draw matplotlib plots within this notebook.\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the following values of the parameters (but of course we are free to use other values in this TD!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 4.0\n",
    "V = 1.0\n",
    "T = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will test our codes with the initial function $u_0(x) = \\max(0, 1-|x|)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u0(x):\n",
    "    return np.maximum(0, 1 - np.abs(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 0: Stability notions\n",
    "\n",
    "In this part we recall a few notions concerning the **stability of a numerical scheme**.\n",
    "\n",
    "### Definition of the stability of a linear numerical scheme\n",
    "\n",
    "For a given numerical solution $u(t,x)$, we will note $u^n = (u_j^n)_{0\\leq j < N}$ the vector in $\\mathbb{R}^n$ space of the computed values of the solution at a time $t_n$. A linear (one-step) scheme implies the existence of a matrix $A$ of shape $N\\times N$ such as:\n",
    "$$\n",
    "    \\forall n\\in \\mathbb{N}, u^{n+1} = A u^n\n",
    "$$\n",
    "Hence, we have (be careful: although there have similar notations, $n$ is an index in the left hand side while $n$ is a power in the right-hand side):\n",
    "$$\n",
    "    \\forall n\\in \\mathbb{N}, u^n = A^n u^0\n",
    "$$\n",
    "\n",
    "From a numerical point of view, what we **DO NOT** want is **divergence**! Thus, for a given norm $||\\text{ . }||$ in $\\mathbb{R}^N$, we will say a linear scheme is stable if and only if there exists $K > 0$ independent of $\\Delta t$ and $\\Delta x$ such as: $\\forall n\\in \\mathbb{N}, ||A^n||\\leq K$ (where the matrix norm is the norm induced from $\\mathbb{R}^n$ norm $||\\text{ . }||$). A stable scheme would then satisfy: $\\forall n\\in \\mathbb{N}, ||u^n|| \\leq K||u^0||$, meaning no divergence!\n",
    "\n",
    "### Von Neumann stability analysis\n",
    "\n",
    "For the specific case of $L^2$ norm stability, we can reason in Fourier space: this is called **Von Neumann stability analysis**. We will admit and apply the following method:\n",
    "1. We write $u_j^n = z^ne^{ikx_j}$ (with $z \\in \\mathbb{C}$ playing the role of a Fourier amplification factor for a given mode $k$)\n",
    "2. We inject the previous expression in the numerical scheme equation.\n",
    "3. The numerical scheme will be stable (in $L^2$ norm) if and only if: $\\forall k, |z| \\leq 1$.\n",
    "\n",
    "The aim of Von Neumann stability analysis is therefore to find some conditions on the numerical parameters ($\\Delta t, \\Delta x, V, ...$) so that $\\forall k, |z| \\leq 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: An explicit centered scheme\n",
    "\n",
    "In this part, we will work on the following **explicit centered scheme**:\n",
    "$$\n",
    "    \\frac{u_j^{n+1} - u_j^n}{\\Delta t} + V\\frac{u_{j+1}^n - u_{j-1}^n}{2\\Delta x} = 0\n",
    "$$\n",
    "\n",
    "**1)** Numerically solve the advection equation using this explicit centered scheme on the spatial domain $[-L, L]$. We will plot on the same figure the initial function $u_0$ and the values of the solution at the time $T$. We can choose several values of $\\Delta t$ and $\\Delta x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have observed the unstability of this scheme.\n",
    "\n",
    "**2)** Show through a **Von Neumann stability analysis** that this scheme is unconditionally unstable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: The Lax-Friedrichs scheme: a diffusive scheme\n",
    "\n",
    "In this part, we will work on the **Lax-Friedrichs scheme** which is defined as follows:\n",
    "$$\n",
    "    \\frac{2u_j^{n+1} - u_{j+1}^n - u_{j-1}^n}{2\\Delta t} + V\\frac{u_{j+1}^n - u_{j-1}^n}{2\\Delta x} = 0\n",
    "$$\n",
    "\n",
    "**3)** Is this scheme stable? We will perform a **Von Neumann stability analysis** of this scheme and show that this scheme is stable if the following CFL condition is respected:\n",
    "$$\n",
    "    |V|\\Delta t \\leq \\Delta x\n",
    "$$\n",
    "\n",
    "**4)** Numerically solve the advection equation using the Lax-Friedrichs scheme on the spatial domain $[-L, L]$. We will first choose values of $V, \\Delta t$ and $\\Delta x$ that meet the previous CFL condition, then we will try values that do not meet this condition to observe an unstability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5)** Perform a **Taylor expansion** of the Lax-Friedrichs scheme around a point $(t_n, x_j)$ for an exact solution of the advection equation $u(t,x)$ to the first non-vanishing order. We should get (when the CFL condition is verified):\n",
    "$$\n",
    "    \\frac{2u (t_{n+1}, x_j) - u (t_{n}, x_{j+1}) - u (t_{n}, x_{j-1})}{2\\Delta t} + V\\frac{u (t_{n}, x_{j+1}) - u (t_{n}, x_{j-1})}{2\\Delta x} = - \\frac{(\\Delta x)^2}{2\\Delta t}\\left(1-\\frac{(V\\Delta t)^2}{(\\Delta x)^2}\\right)\\frac{\\partial^2u}{\\partial x^2}(t_n, x_j) + \\mathcal{O}\\left((\\Delta x)^2 + \\frac{(\\Delta x)^4}{\\Delta t}\\right)\n",
    "$$\n",
    "\n",
    "\n",
    "This shows that the Lax-Friedrichs scheme is also a good approximation of the following convection-diffusion equation:\n",
    "$$\n",
    "    \\frac{\\partial u}{\\partial t} + V \\frac{\\partial u}{\\partial x} - \\nu \\frac{\\partial^2u}{\\partial x^2} = 0 \\text{ with } \\nu = \\frac{(\\Delta x)^2}{2\\Delta t}\\left(1-\\frac{(V\\Delta t)^2}{(\\Delta x)^2}\\right)\n",
    "$$\n",
    "The third term of this equation corresponds to a diffusive term, where $\\nu$ plays the role of a diffusion coefficient. This scheme is therefore a **diffusive scheme**.\n",
    "\n",
    "**6)** Choosing values of $V, \\Delta t$ and $\\Delta x$ so that $\\nu$ is big enough, observe this diffusive behaviour of the Lax-Friedrichs scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: The Lax-Wendroff scheme: a dispersive scheme\n",
    "\n",
    "In this part, we will work on the **Lax-Wendroff scheme** which is defined as follows:\n",
    "$$\n",
    "    \\frac{u_j^{n+1} - u_j^n}{\\Delta t} + V\\frac{u_{j+1}^n - u_{j-1}^n}{2\\Delta x} - \\left(\\frac{V^2\\Delta t}{2}\\right)\\frac{u_{j-1}^n - 2u_j^n + u_{j+1}^n}{(\\Delta x)^2} = 0\n",
    "$$\n",
    "\n",
    "**7)** Perform a **Von Neumann stability analysis** of this scheme and exhibit a CFL condition for its stability.\n",
    "\n",
    "**8)** Numerically solve the Lax-Wendroff scheme on the spatial domain $[-L, L]$ showing at least one stable behavior and one unstable behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**9)** Perform a **Taylor expansion** of the Lax-Wendroff scheme around a point $(t_n, x_n)$ for an exact solution of the advection equation $u(t,x)$ to the first non-vanishing order. We should get (when the CFL condition is verified):\n",
    "$$\n",
    "    \\frac{u (t_{n+1}, x_j) - u (t_{n}, x_{j})}{\\Delta t} + V\\frac{u (t_{n}, x_{j+1}) - u (t_{n}, x_{j-1})}{2\\Delta x} - \\left(\\frac{V^2\\Delta t}{2}\\right)\\frac{u (t_{n}, x_{j-1}) - 2u (t_{n}, x_{j}) + u (t_{n}, x_{j+1})}{(\\Delta x)^2} = \\frac{V\\Delta x^2}{6}\\left(1-\\frac{(V\\Delta t)^2}{\\Delta x^2}\\right)\\frac{\\partial^3u}{\\partial x^3}(t_n, x_j) + \\mathcal{O}(\\Delta t^3 + \\Delta x^3)\n",
    "$$\n",
    "\n",
    "This shows that the Lax-Wendroff scheme is also a good approximation of the following equation:\n",
    "$$\n",
    "    \\frac{\\partial u}{\\partial t} + V \\frac{\\partial u}{\\partial x} + \\frac{V(\\Delta x)^2}{6}\\left(1 - \\frac{(V\\Delta t)^2}{(\\Delta x)^2}\\right)\\frac{\\partial^3u}{\\partial x^3} = 0\n",
    "$$\n",
    "The third term of this equation corresponds to a dispersive term. This scheme is therefore a **dispersive scheme**.\n",
    "\n",
    "**10)** Choosing values of $V, \\Delta t$ and $\\Delta x$ so that the dispersive term can be big enough, observe this dispersive behaviour of the Lax-Wendroff scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
