Tutorial-How_NRPy_Computes_Finite_Difference_Coeffs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# How NRPy+ Computes Finite Difference Coefficients\n",
    "\n",
    "## Author: Zach Etienne\n",
    "\n",
    "### NRPy+ Source Code for this module: [finite_difference.py](../edit/finite_difference.py)\n",
    "\n",
    "## Introduction to Finite Differencing: Basic Definitions\n",
    "\n",
    "Finite difference derivatives are written in the form\n",
    "$$u^{(n)}(x_i)_{\\text{FD}} = \\sum_{j} u_j a_j,$$\n",
    "where the $a_j$'s are known as *finite difference coefficients*. For a given finite difference representation, the set of $a_j$'s are unique.\n",
    "\n",
    "There are many ways to compute finite difference coefficients $a_j$, though perhaps the most popular involves Taylor series expansions about sampled points near the point where we wish to evaluate the derivative. Below we derive the basic algorithm for a specific case, which is trivially generalizable. The generalized version is implemented in NRPy+'s **compute_fdcoeffs_fdstencl()** function, in the finite_difference module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Centered Finite Difference Representation of $u'(x_0) = u'_0$ Accurate to Fourth-order in $\\Delta x$\n",
    "\n",
    "As an illustration, let's first derive for uniform grids the centered, first-order, finite-difference coefficients accurate to fourth-order in $\\Delta x$. The fourth-order-accurate, uniformly sampled, centered finite-difference derivative $u'(x_0)$ is equivalent to the derivative of the unique polynomial that passes through $u(x)$ at sampled points $\\left\\{ x_{-2},x_{-1},x_{0},x_{1},x_{2} \\right\\}$, where $x_i=x_0 + i \\Delta x$.\n",
    "\n",
    "The Taylor series expansion of a function $u(x)$ about a point $x_0$ is given by\n",
    "\n",
    "$$u(x) = \\sum_{n=0}^\\infty \\frac{u^{(n)}(x_0)}{n!} (x-x_0)^n,$$\n",
    "\n",
    "where $u^{(n)}(x_0)$ is the $n$th derivative of $u$ evaluated at point $x_0$. Based on this, we can immediately write the Taylor expansion of $f$ at a point $x=x_0+j\\Delta x$. In this case,\n",
    "\n",
    "\\begin{align}\n",
    "u(x_0+j\\Delta x) &= \\sum_{n=0}^\\infty \\frac{u^{(n)}(x_0)}{n!} (j\\Delta x)^n,\\text{ or equivalently:} \\\\\n",
    "u_j &= \\sum_{n=0}^\\infty \\frac{u^{(n)}_0}{n!} (j\\Delta x)^n.\n",
    "\\end{align}\n",
    "\n",
    "Our goal is to compute $u^{(1)}(x_0)=u'_0$ at some point $x_0$, with a dominant error term proportional to $(\\Delta x)^4$. We accomplish this by Taylor expanding $u(x_j)$ about $x=x_0$ for $j\\in \\left\\{-2,-1,0,1,2\\right\\}$, each up to the $n=5$ term:\n",
    "\n",
    "\\begin{align}\n",
    "u_{-2} &= u_0 - (2 \\Delta x) u'_0 + \\frac{(2 \\Delta x)^2}{2} u''_0 - \\frac{(2 \\Delta x)^3}{3!} u'''_0 + \\frac{(2 \\Delta x)^4}{4!} u^{(4)}_0 +\\mathcal{O}\\left((\\Delta x)^5\\right) \\\\\n",
    "u_{-1} &= u_0 - (\\Delta x) u'_0 + \\frac{(\\Delta x)^2}{2} u''_0 - \\frac{(\\Delta x)^3}{3!} u'''_0 + \\frac{(\\Delta x)^4}{4!} u^{(4)}_0 +\\mathcal{O}\\left((\\Delta x)^5\\right)\\\\\n",
    "u_{0} &= u_0 \\\\\n",
    "u_{1} &= u_0 + (\\Delta x) u'_0 + \\frac{(\\Delta x)^2}{2} u''_0 + \\frac{(\\Delta x)^3}{3!} u'''_0 + \\frac{(\\Delta x)^4}{4!} u^{(4)}_0 +\\mathcal{O}\\left((\\Delta x)^5\\right)\\\\\n",
    "u_{2} &= u_0 + (2 \\Delta x) u'_0 + \\frac{(2 \\Delta x)^2}{2} u''_0 + \\frac{(2 \\Delta x)^3}{3!} u'''_0 + \\frac{(2 \\Delta x)^4}{4!} u^{(4)}_0 +\\mathcal{O}\\left((\\Delta x)^5\\right)\\\\\n",
    "\\end{align}\n",
    "\n",
    "Let's combine the above equations to find coefficients $a_j$ such that $(a_{-2} u_{-2} + a_{-1} u_{-1}...)/(\\Delta x) = u'_0 + \\mathcal{O}\\left((\\Delta x)^4\\right)$:\n",
    "\n",
    "\\begin{align}\n",
    "& (a_{-2} u_{-2} + a_{-1} u_{-1} + a_0 u_0 + a_{1} u_{1} +a_{2} u_{2})/(\\Delta x) \\\\\n",
    "= & \\left( u_0 - (2 \\Delta x) u'_0 + \\frac{(2 \\Delta x)^2}{2} u''_0 -\\frac{(2 \\Delta x)^3}{3!} u'''_0+\\frac{(2 \\Delta x)^4}{4!} u^{(4)}_0 \\right) a_{-2} \\\\\n",
    "& + \\left( u_0 - (\\Delta x) u'_0 + \\frac{(\\Delta x)^2}{2} u''_0 - \\frac{(\\Delta x)^3}{3!} u'''_0+\\frac{(\\Delta x)^4}{4!} u^{(4)}_0 \\right) a_{-1} \\\\\n",
    "& + \\left( u_0 \\right) a_{0} \\\\\n",
    "& + \\left( u_0 + (\\Delta x) u'_0 + \\frac{(\\Delta x)^2}{2} u''_0 + \\frac{(\\Delta x)^3}{3!} u'''_0+\\frac{(\\Delta x)^4}{4!} u^{(4)}_0 \\right) a_{1} \\\\\n",
    "& + \\left( u_0 + (2 \\Delta x) u'_0 + \\frac{(2 \\Delta x)^2}{2} u''_0 + \\frac{(2 \\Delta x)^3}{3!} u'''_0+\\frac{(2 \\Delta x)^4}{4!} u^{(4)}_0 \\right) a_{2}\n",
    "\\end{align}\n",
    "\n",
    "First notice that each time we take a derivative in the Taylor\n",
    "expansion, we multiply by a $\\Delta x$. Notice that this helps to keep\n",
    "the units consistent (e.g., if $x$ were in units of meters). Let's\n",
    "just **absorb those $\\Delta x$'s into the derivatives (we will extract them again later)** and\n",
    "rearrange terms:\n",
    "\n",
    "\\begin{align}\n",
    "& a_{-2} u_{-2} + a_{-1} u_{-1} + a_0 u_0 + a_{1} u_{1} + a_{2} u_{2} \\\\\n",
    "& = \\left( a_{-2} + a_{-1} + a_0 + a_{1} + a_{2} \\right) \\times u_0 \\\\\n",
    "& + \\left( -2 a_{-2} - a_{-1} + a_{1} + 2 a_{2} \\right) \\times u'_0 \\\\\n",
    "& + \\left( 2^2 a_{-2} + a_{-1} + a_{1} + 2^2 a_{2} \\right)/2! \\times u''_0 \\\\\n",
    "& + \\left( -2^3 a_{-2} - a_{-1} + a_{1} + 2^3 a_{2} \\right)/3! \\times u'''_0 \\\\\n",
    "= & u'_0\n",
    "\\end{align}\n",
    "\n",
    "In order for the above to hold true for any nonzero values of\n",
    "$\\left\\{ u_0,u'_0,u''_0,u'''_0,u^{(4)}_0\\right\\}$, the following set\n",
    "of equations must also hold:\n",
    "\\begin{align}\n",
    "0 &= a_{-2} + a_{-1} + a_0 + a_{1} + a_{2}\\\\\n",
    "1 &= -2 a_{-2} - a_{-1} + a_{1} + 2 a_{2}\\\\\n",
    "0 \\times 2! &= 2^2 a_{-2} + a_{-1} + a_{1} + 2^2 a_{2}\\\\\n",
    "0 \\times 3! &= -2^3 a_{-2} - a_{-1} + a_{1} + 2^3 a_{2} \\\\\n",
    "0 \\times 4! &= 2^4 a_{-2} + a_{-1} + a_{1} + 2^3 a_{2}.\n",
    "\\end{align}\n",
    "\n",
    "Now we write this expression in matrix form (note that $0!=1$):\n",
    "\\begin{equation}\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "0\\times 0! \\\\\n",
    "1\\times 1! \\\\\n",
    "0\\times 2! \\\\\n",
    "0\\times 3! \\\\\n",
    "0\\times 4! \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "=\n",
    "\\left(\n",
    "\\begin{array}{ccccc}\n",
    " 1 &  1 & 1 & 1 & 1 \\\\\n",
    "(-2)^1 &(-1)^1 & 0 & 1 & 2 \\\\\n",
    "(-2)^2 &(-1)^2 & 0 & 1 & 2^2 \\\\\n",
    "(-2)^3 &(-1)^3 & 0 & 1 & 2^3 \\\\\n",
    "(-2)^4 &(-1)^4 & 0 & 1 & 2^4 \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "a_{-2} \\\\\n",
    "a_{-1} \\\\\n",
    "a_{0} \\\\\n",
    "a_{1} \\\\\n",
    "a_{2} \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\end{equation}\n",
    "\n",
    "So we have reduced the computation of finite difference coefficients\n",
    "to the inversion of an $N\\times N$ matrix equation. Notice that the\n",
    "elements of the matrix will vary from the one given above if the grid\n",
    "spacing is not constant, but are otherwise invariant to $\\Delta x$.\n",
    "\n",
    "The inverted matrix reads\n",
    "\\begin{equation}\n",
    "\\left(\n",
    "\\begin{array}{ccccc}\n",
    "0 & 1/12 & -1/24 & -1/12 & 1/24 \\\\\n",
    "0 & -2/3 & 2/3 & 1/6 & -1/6 \\\\\n",
    "1 & 0 & -5/4 & 0 & 1/4 \\\\\n",
    "0 & 2/3 & 2/3 & -1/6 & -1/6 \\\\\n",
    "0 & -1/12 & -1/24 & 1/12 & 1/24 \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\label{fourthorder_inv_matrix}\n",
    "\\end{equation}\n",
    "\n",
    "The coefficients for the $M$th derivative can be immediately read by\n",
    "multiplying the $(M+1)$st column by $M!/(\\Delta x)^M$. For example, the zeroth derivative at point $x_0$ is given by \n",
    "\n",
    "$$\\frac{0!}{(\\Delta x)^0} \\times (0 u_{-2} + 0 u_{-1} + u_0 + 0 u_{1} + 0 u_{2}) = u_0,$$\n",
    "\n",
    "which is exact. The first derivative finite difference approximation at point $x_0$ is given by\n",
    "\n",
    "$$\\frac{1!}{(\\Delta x)^1} \\times \\left(\\frac{1}{12}( u_{-2} - u_{2}) + \\frac{2}{3}( -u_{-1} + u_{1})\\right) \\approx (\\partial_x u)_0,$$\n",
    "\n",
    "and the second derivative finite difference approximation at point $x_0$ is given by \n",
    "\n",
    "$$\\frac{2!}{(\\Delta x)^2} \\times \\left(-\\frac{1}{24}(u_{-2} + u_{2}) + \\frac{2}{3}(u_{-1} + u_{1}) - \\frac{5}{4} u_0 \\right) \\approx (\\partial_x^2 u)_0,$$\n",
    "\n",
    "In short, this matrix yields the finite difference derivative coefficients with the lowest possible error given a stencil size of 5 gridpoints. It can be shown by analyzing cancellations in higher order terms of the Taylor series expansions that the first and second derivative coefficients are correct to $(\\Delta x)^4$ and third and fourth derivatives are correct to $(\\Delta x)^2$. \n",
    "\n",
    "### Exercise 1: Find the exact expressions for the dominant error term on all derivatives that can be computed from this matrix (zeroth through fourth derivatives).\n",
    "\n",
    "### Exercise 2: Construct the matrix whose inverse yields the 5-point stencil *upwinded* derivative coefficients (i.e., the stencil includes points $\\{u_{-4},u_{-3},u_{-2},u_{-1},u_{0}\\}$).\n",
    "\n",
    "NRPy+ implements this simple matrix inversion strategy to evaluate finite difference coefficients."
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}