{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TD4: 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 compute 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": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"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$ in 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 a 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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}