{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stability of the heat equation FD approximation\n",
    "\n",
    "## Calcuation of the stability condtions\n",
    "\n",
    "The explicit (forward difference in time) formulation of the FD heat equation has a maximum time step that can be used in the calculations. Exceeding this time step will lead to the *instability* of the calculations, and typically, misbehaviour or crashing of the computational code (and, certainly, erroneous results).\n",
    "\n",
    "Consider the simplified heat equation\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{\\partial T}{\\partial t} = \\kappa \\frac{\\partial^2 T}{\\partial x^2}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\kappa = \\frac{k}{\\rho C_p}$. The FD approximation is\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{T_i^{q+1}-T_i^q}{\\Delta t} = \\kappa \\frac{T_{i+1}^{q} - 2T_{i}^{q} + T_{i-1}^{q}}{(\\Delta x)^2}\n",
    "\\end{equation}\n",
    "\n",
    "and, with $r=\\frac{\\kappa\\Delta t}{\\Delta x^2}$, can be written\n",
    "\n",
    "\\begin{equation}\n",
    "  T_i^{q+1} = T_i^q + r (T_{i+1}^{q} - 2T_{i}^{q} + T_{i-1}^{q})\n",
    "\\end{equation}\n",
    "\n",
    "If we assume that the numerical solution of the temperature, $T$ is the sum of the correct temperature values $\\theta$ and the error $\\epsilon$ at each grid point ($T_i = \\theta_i + \\epsilon_i$), the difference equations of temperatures and errors can be separated, leading to\n",
    "\n",
    "\\begin{equation}\n",
    "  \\epsilon_i^{q+1} = \\epsilon_i^q + r (\\epsilon_{i+1}^{q} - 2\\epsilon_{i}^{q} + \\epsilon_{i-1}^{q})\n",
    "\\end{equation}\n",
    "\n",
    "That is, the error *diffuses* from grid point to next grid point, similar to the actual temperature. We want to *suppress* the error, that is, we want\n",
    "\n",
    "\\begin{equation}\n",
    "  |\\epsilon_i^{q+1}| \\leq |\\epsilon_i^q| \\Rightarrow \\Big|\\frac{\\epsilon_i^{q+1}}{\\epsilon_i^q}\\Big| \\leq 1,\n",
    "\\end{equation}\n",
    "\n",
    "since otherwise the error would grow from one time step to the next. \n",
    "\n",
    "Dividing the error diffusion expression by $\\epsilon_i^q$ and applying the inequality\n",
    "\n",
    "\\begin{equation}\n",
    "  \\Big|1 + r (\\epsilon_{i+1}^{q} - 2\\epsilon_{i}^{q} + \\epsilon_{i-1}^{q}) / \\epsilon_i^q\\Big| \\leq 1\n",
    "\\end{equation}\n",
    "\n",
    "Let's consider the worst case scenario, where $\\epsilon_{i+1} = \\epsilon_{i-1} = a$ and $\\epsilon_i = -a$. We get\n",
    "\n",
    "\\begin{equation}\n",
    "  \\Big|1 - 4r\\Big| \\leq 1 \\Rightarrow 0 \\leq r \\leq \\frac{1}{2}\n",
    "\\end{equation}\n",
    "\n",
    "The first inequality is naturally true ($r$ cannot be negative), the second leads to\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{\\kappa \\Delta t}{\\Delta x^2} \\leq \\frac{1}{2}\n",
    "\\end{equation}\n",
    "\n",
    "This inequality can be used to determine a proper time step if grid spacing has been chosen, or proper grid spacing if time step has been chosen. Typically, it is safer to enforce an even stricter condition, e.g., $\\frac{\\kappa \\Delta t}{\\Delta x^2} \\leq \\frac{1}{4}$.\n",
    "\n",
    "**Note**: The stability of the method does not guarantee its correctness!  \n",
    "\n",
    "A physical interpretation exists for the stability criterion: The diffusive time scale $\\sqrt{\\kappa\\Delta t}$ should not exceed the grid spacing, i.e. heat should not diffuse in significant amounts over more than one grid point within one time step."
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}