lecture-7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lecture 7: Spectral Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## An idea\n",
    "\n",
    "* Suppose we are solving\n",
    "$$\n",
    "\\begin{align*}\n",
    "-\\Delta u &= f \\\\\n",
    "u &= 0\\quad\\text{on }\\Gamma,\n",
    "\\end{align*}\n",
    "$$\n",
    "in $\\Omega = [0,1]^2$ (unit square)\n",
    "\n",
    "* Idea: why not\n",
    "$$\n",
    "v^{(k,\\ell)} := \\sin(k \\pi x) \\sin(\\ell \\pi y)\n",
    "$$\n",
    "as basis functions?\n",
    "\n",
    "* It is a good idea:\n",
    "$$\n",
    "-\\Delta v^{(k,\\ell)} = \\pi^2 (k^2+\\ell^2) v^{(k,\\ell)}\n",
    "$$\n",
    "(in other words, $v^{(k,\\ell)}$ are eigenvalues of $-\\Delta$)\n",
    "\n",
    "* Can use Fast Sine Transform (FST) to obtain the coefficients in\n",
    "$$\n",
    "u^h = \\sum _{k,\\ell} c _{k,\\ell} v^{(k,\\ell)}\n",
    "$$\n",
    "\n",
    "* We then can easily solve:\n",
    "    - Sample $f$ at points $(i h, j h)$, $h=1/N$;\n",
    "    - Find the coefficients $\\hat{f} _{k,\\ell}$ in\n",
    "    $$\n",
    "        f = \\sum _{k,\\ell} \\hat{f} _{k,\\ell} v^{(k,\\ell)}\n",
    "    $$\n",
    "    (i.e., apply the FST)\n",
    "    \n",
    "    - Find the coefficients of the solution\n",
    "    $$\n",
    "        {\\pi^2(k^2+\\ell^2)} \\hat{u}^h _{k,\\ell} = \\hat{f} _{k,\\ell}\n",
    "    $$\n",
    "    - Find the solution (i.e., apply the inverse FST)\n",
    "    $$\n",
    "        u _h = \\sum _{k,\\ell} \\hat{u}^h _{k,\\ell} v^{(k,\\ell)}\n",
    "    $$\n",
    "    \n",
    "* **To think:** Did we discretize the domain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Algorithm:\n",
    "\n",
    "Given $f$ and $N$ (hence the size of the discrete problem is $(N-1)^2$)\n",
    "* Step 1: calculate $f(ih, jh)$\n",
    "* Step 2: find the coefficients $\\hat{f}_{k,l}$ by applying the fast DST\n",
    "* Step 3: calculate $\\hat{u}^h_{k,l}$\n",
    "* Step 4: find $u_h$ by applying the fast inverse DST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Smoothness in x  <->  fast decay in the Fourier domain\n",
    "\n",
    "* Smooth in Fourier domain <-> decays fast in x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Time-dependent example\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "u _t &= \\Delta u \\\\\n",
    "u &= u^0\\quad\\text{at }t=0 \\\\\n",
    "u &= 0\\quad\\text{on }\\Gamma,\n",
    "\\end{align*}\n",
    "$$\n",
    "in $\\Omega = [0,1]^2$ (unit square)\n",
    "\n",
    "* If we assume the same spatial discretization, then\n",
    "$$\n",
    "u _h(t,x,y) = \\sum _{k,\\ell} e^{-\\pi^2 (k^2+\\ell^2) t} \\hat{u}^0 _{k,\\ell} v^{(k,\\ell)}\n",
    "$$\n",
    "\n",
    "* The algorithm is then:\n",
    "    - Find $\\hat{u}^0 _{k,\\ell}$ (apply the FST)\n",
    "    - Calculate $\\hat{u}^h _{k,\\ell}(T) := e^{-\\pi^2 (k^2+\\ell^2) T}\\, \\hat{u}^0 _{k,\\ell}$\n",
    "    - Hence find $u^h(T)$ (apply the inverse FST)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Exercise:\n",
    "\n",
    "$$ u_t = a u_x $$\n",
    "$$ u(t,x+1) = u(t,x)$$\n",
    "$$\\text{(periodic boundary conditions)} $$\n",
    "\n",
    "1) write down the three schemes:\n",
    "\n",
    "* An explicit scheme\n",
    "* An implicit scheme\n",
    "* A spectral method\n",
    "\n",
    "2) Study their stability (write down $\\lambda^{(k)}$ for the time step operator)\n",
    "\n",
    "3) Hint: use that $v^{(k)}(x) = \\exp(i \\pi k x)$ is a (complex-valued) eigenfunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other uses of the spectral methods (1):\n",
    "\n",
    "* Suppose we need to compute\n",
    "$$\n",
    "{\\rm div}(K \\nabla u)\n",
    "$$\n",
    "with non-constant $K$ and with (for a change) periodic boundary conditions (so FFT can be applied) in a unit square.\n",
    "\n",
    "* The algorithm is then:\n",
    "    - Given values of $u$\n",
    "    - Calculate $\\hat{u} = {\\rm FFT}(u)$\n",
    "    - Calculate $i\\pi k \\hat{u} _{k,\\ell} = {\\rm FFT}(\\frac{\\partial}{\\partial x} u)$\n",
    "    - Calculate $\\frac{\\partial}{\\partial x} u = {\\rm iFFT}(i\\pi k \\hat{u} _{k,\\ell})$\n",
    "    - Then multiply by $K$\n",
    "    - and again apply $\\frac{\\partial}{\\partial x}$ using ${\\rm FFT}$ and  ${\\rm iFFT}$.\n",
    "    - Repeat with $\\frac{\\partial}{\\partial y}$.\n",
    "    \n",
    "* Nonlinear problems can be treated in a similar way\n",
    "    - But extra care must be exercised (check out, e.g., what is an _aliasing error_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example (2): Nonlinear Schrodinger equation\n",
    "\n",
    "$$ i u _t = -\\frac{1}{2} u _{xx} + |u| ^2 u $$\n",
    "\n",
    "* Split-step Fourier method (physical splitting + spectral for the derivatives)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other uses of the spectral methods (3):\n",
    "\n",
    "* Can use polynomials instead of trigonometry:\n",
    "$$\n",
    "T _n(x) = \\cos(n \\arccos (x))\n",
    "$$\n",
    "are Chebyshev polynomials!\n",
    "If one generates meshes with nodes $x _k = \\arccos (h k)$, one can use FFT!\n",
    "    - Check out [ChebFun](http://www.chebfun.org/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Other uses of the spectral methods (4):\n",
    "\n",
    "* Can split the domain into several quadrilaterals, and use a spectral discretization in each $\\rightarrow$ **Spectral element method**\n",
    "<img src=\"fig/SEM.gif\">\n",
    "\n",
    "* Closely related to the hp-FEM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Limitations of the Spectral Method:\n",
    "\n",
    "* Efficient only if the solution is smooth (otherwise not more accurate than a FEM of some order)\n",
    "* Can be used only in rectangular domains\n",
    "    - or smooth deformations of rectangular domains\n",
    "    - or their union\n",
    "* Parallelization is not easy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Questions?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}