08_convergence_IVP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Convergence Results for Initial Value Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Convergence for IVPs is a bit different than with BVPs, we want in general\n",
    "$$\n",
    "    \\lim_{\\Delta t \\rightarrow 0} U^N = u(t_f)\n",
    "$$\n",
    "where $t_f$ is the final desired time and $N$ is the number of time steps needed to reach $t_f$ such that\n",
    "$$\n",
    "    N \\Delta t = t_f \\quad \\Rightarrow N = \\frac{t_f}{\\Delta t}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We need to be careful at this juncture however when we are talking about a convergent method.  A method can be convergent for a particular set of equations and particular initial conditions but not others.  Practically speaking we would like convergence results to apply to a reasonably large set of equations and initial conditions.  With these considerations we have the following definition of convergence for IVPs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If applying an $r$-step method to an ODE of the form\n",
    "$$\n",
    "    u'(t) = f(t,u)\n",
    "$$\n",
    "with $f(t,u)$ Lipschitz continuous in $u$, and with any set of starting values satisfying\n",
    "$$\n",
    "    \\lim_{\\Delta t\\rightarrow 0} U^\\nu(\\Delta t) = u_0 \\quad \\text{for} \\quad \\nu = 0, 1, \\ldots, r-1\n",
    "$$\n",
    "(i.e. the bootstrap startup for the multi-step method is consistent with the initial value as $\\Delta t \\rightarrow$), then the method is said to be *convergent* in the sense\n",
    "$$\n",
    "    \\lim_{\\Delta t \\rightarrow 0} U^N = u(t_f).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As we saw previously for a method to be convergent it must be\n",
    " - **consistent** - the local truncation error $\\tau = \\mathcal{O}(\\Delta t^p)$ where $p > 0$ and\n",
    " - **zero-stable** - a similar minimal form of stability implying that the sum total of the errors as $\\Delta t \\rightarrow 0$ is bounded and has the same order as $\\tau$ which we know goes to zero as $\\Delta t \\rightarrow 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## One-Step Method Convergence\n",
    "\n",
    "Consider the simple linear problem\n",
    "$$\n",
    "    \\frac{\\text{d}u}{\\text{d}t} = \\lambda u + g(t) \\quad \\text{with}\\quad u(0) = u_0\n",
    "$$\n",
    "which we know has the solution \n",
    "$$\n",
    "    u(t) = u_0 e^{\\lambda (t - t_0)} + \\int^t_{t_0} e^{\\lambda (t - \\tau)} g(\\tau) d\\tau.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Forward Euler on a Linear Problem\n",
    "\n",
    "Applying Euler's method to this problem leads to \n",
    "$$\\begin{aligned}\n",
    "    U^{n+1} &= U^n + \\Delta t\\lambda U^n \\\\\n",
    "            &= (1 + \\Delta t \\lambda) U^n\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We also know the local truncation error is\n",
    "$$\\begin{aligned}\n",
    "    \\tau^n &= \\left (\\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} \\right ) - \\lambda u(t_n)\\\\\n",
    "    &= \\left (u'(t_n) + \\frac{1}{2} \\Delta t u''(t_n) + \\mathcal{O}(\\Delta t^2) \\right ) - u'(t_n) \\\\\n",
    "    &= \\frac{1}{2} \\Delta t u''(t_n) + \\mathcal{O}(\\Delta t^2)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Noting the original definition of $\\tau^n$ we can rewrite the expression for the local truncation error as\n",
    "\n",
    "$$\n",
    "    u(t_{n+1}) = (1 + \\Delta t \\lambda) u(t_n) + \\Delta t \\tau^n\n",
    "$$\n",
    "\n",
    "which in combination with the application of Euler's method leads to an expression for the global error\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    E^{n+1} = U^{n+1} - u(t^{n+1}) &= (1 + \\Delta t \\lambda) U^n - (1 + \\Delta t \\lambda) u(t_n) - \\Delta t \\tau^n \\\\\n",
    "    &= (1+\\Delta t \\lambda) E^n - \\Delta t \\tau^n \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Expanding this expression out backwards in time to $n=0$ leads to \n",
    "\n",
    "$$\n",
    "    E^n = (1 + \\Delta t \\lambda) E^0 - \\Delta t \\sum^n_{i=1} (1 + \\Delta t \\lambda)^{n-i} \\tau^{i - 1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now see the importance of the term $(1 + \\Delta t \\lambda)$.  We can bound this term by\n",
    "\n",
    "$$\n",
    "    |1 + \\Delta t \\lambda| \\leq e^{\\Delta t \\lambda}\n",
    "$$\n",
    "\n",
    "which then implies the term in the summation can be bounded by\n",
    "\n",
    "$$\n",
    "    |1 + \\Delta t \\lambda|^{n - i} \\leq e^{(n-i) \\Delta t |\\lambda|} \\leq e^{n \\Delta t |\\lambda||} \\leq e^{|\\lambda| t_f}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using this expression in the expression for the global error we find\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    E^n &= (1 + \\Delta t \\lambda) E^0 - \\Delta t \\sum^n_{i=1} (1 + \\Delta t \\lambda)^{n-i} \\tau^{i - 1} \\\\\n",
    "    |E^n| &\\leq e^{|\\lambda| \\Delta t} |E^0| - \\Delta t \\sum^n_{i=1} e^{|\\lambda| t_f} |\\tau^{i - 1}| \\\\\n",
    "          &\\leq e^{|\\lambda| t_f} \\left(|E^0| - \\Delta t \\sum^n_{i=1} |\\tau^{i - 1}|\\right) \\\\\n",
    "          &\\leq e^{|\\lambda| t_f} \\left(|E^0| - n \\Delta t \\max_{1 \\leq i \\leq n} |\\tau^{i - 1}|\\right)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In other words the global error is bounded by the original global error and the maximum one-step error made multiplied by the number of time steps taken.  If $N = \\frac{t_f}{\\Delta t}$ as before and taking into account the local truncation error we can simplify this expression further to\n",
    "\n",
    "$$\n",
    "    |E^n| \\leq e^{|\\lambda| t_f} \\left[|E^0| + t_f \\left(\\frac{1}{2} \\Delta t |u''| + \\mathcal{O}(\\Delta t^2)\\right ) \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we assume that we have used the correct initial condition $u_0$ then $E_0 \\rightarrow 0$ as $\\Delta t \\rightarrow 0$ and we see that the method is truly convergent as \n",
    "\n",
    "$$\n",
    "    |E^n| \\leq e^{|\\lambda| t_f} t_f \\left(\\frac{1}{2} \\Delta t |u''| + \\mathcal{O}(\\Delta t^2)\\right ) = \\mathcal{O}(\\Delta t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Relation to Stability for BVPs\n",
    "\n",
    "We can see the relationship between the previous version of stability and the one outlined above.  Try writing the forward Euler method as a linear system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Forward Euler:\n",
    "$$\n",
    "    A = \\frac{1}{\\Delta t} \\begin{bmatrix}\n",
    "        1 \\\\\n",
    "        -(1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        & -(1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        & & -(1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        & & & \\ddots & \\ddots \\\\\n",
    "        & & & & -(1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        & & & & & -(1 + \\Delta t \\lambda) & 1\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "    U = \\begin{bmatrix} U^1 \\\\ U^2 \\\\ \\vdots \\\\ U^N \\end{bmatrix} ~~~~ \n",
    "    F = \\begin{bmatrix} (1 / \\Delta t + \\lambda) U^0 + g(t_0) \\\\ g(t_1) \\\\ \\vdots \\\\ g(t_{N-1}) \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Following our previous stability result and taking $\\hat{U~}$ to be the vector obtained from the true solution ($\\hat{U~}^i = u(t_i)$) we then have\n",
    "$$\n",
    "    A U = F  ~~~~~~ A \\hat{U~} = F + \\tau\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "    A (\\hat{U~} - U) = \\tau.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Noting that $\\hat{U~} - U = E$ we can then invert that matrix $A$ to find the relationship between the truncation error $\\tau$ and the global error $E$.  As before then we require that $A^{-1}$ is invertible (which is trivial in this case) and that $||A^{-1}|| < C$ in some norm.  We can see this as\n",
    "$$\n",
    "    A^{-1} = \\Delta t \\begin{bmatrix}\n",
    "        1 \\\\\n",
    "        (1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        (1 + \\Delta t \\lambda)^2 & (1 + \\Delta t \\lambda) & 1 \\\\\n",
    "        (1 + \\Delta t \\lambda)^3 & (1 + \\Delta t \\lambda)^2 & (1 + \\Delta t \\lambda) & 1\\\\\n",
    "        \\vdots & & & \\ddots \\\\\n",
    "        (1 + \\Delta t \\lambda)^{N-1} & (1 + \\Delta t \\lambda)^{N-2} (1 + \\Delta t \\lambda)^{N-3} & \\cdots & (1 + \\Delta t \\lambda) & 1 \n",
    "   \\end{bmatrix}\n",
    "$$\n",
    "whose infinity norm is\n",
    "$$\n",
    "    ||A^{-1}||_\\infty = \\Delta t \\sum^N_{m=1} | (1 + \\Delta t \\lambda)^{N-M} |\n",
    "$$\n",
    "and therefore\n",
    "$$\n",
    "    ||A^{-1}||_\\infty \\leq \\Delta t N e^{|\\lambda| T} = T e^{|\\lambda| T}.\n",
    "$$\n",
    "As $\\Delta t \\rightarrow 0$ this is bounded for **fixed T**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### General One-Step Method Convergence\n",
    "\n",
    "Consider the general one step method denoted by\n",
    "$$\n",
    "    U^{n+1} = U^n + \\Delta t \\Phi(U^n, t_n, \\Delta t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Assuming $\\Phi$ is continuous in $t$ and $\\Delta t$ and Lipschitz continous in $u$ with Lipschitz contsant $L$ (related to the Lipschitz constant of $f$).  If the one-step method is consistent\n",
    "$$\n",
    "    \\Phi(u,t,0) = f(u,t)\n",
    "$$\n",
    "for all $u$, $t$, and $\\Delta t$ and the local truncation error is\n",
    "$$\n",
    "    \\tau^n =\\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - \\phi(u(t_n), t_n, \\Delta t)\n",
    "$$\n",
    "then the one-step method is convergent.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the general approach we used for forward Euler we know that the true solution and $\\tau$ are realted through\n",
    "$$\n",
    "    u(t_{n+1}) = u(t_n) + \\Delta t \\Phi(u(t_n), t_n, \\Delta t) + \\Delta t \\tau^n\n",
    "$$\n",
    "which subtracted from the approximate solution\n",
    "$$\n",
    "    U^{n+1} = U^n + \\Delta t \\Phi(U^n, t_n, \\Delta t)\n",
    "$$\n",
    "leads to\n",
    "$$\n",
    "    E^{n+1} = E^n + \\Delta t (\\Phi(U^n, t_n, \\Delta t) - \\Phi(u(t_n), t_n, \\Delta t)) - \\Delta t \\tau^n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the Lipschitz continuity of $\\Phi$ we then have\n",
    "$$\n",
    "    |E^{n+1}| \\leq |E^n| + \\Delta t L |E^n| + \\Delta t |\\tau^n|.\n",
    "$$\n",
    "which has the same form as we saw in the proof for forward Euler."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Zero-Stability for Linear Multistep Methods\n",
    "\n",
    "We can also make general statements for linear multistep methods although it is important to note that we have additional requirements for linear multistep methods so that they are convergent.  As an example consider the method\n",
    "$$\n",
    "    U^{n+2} - 3 U^{n+1} + 2 U^n = - \\Delta t f(U^n)\n",
    "$$\n",
    "so that $\\alpha_0 = 2$, $\\alpha_1 = -3$, and $\\alpha_2 = 1$ and $\\beta_0 = -1$ with the rest equal to zero.  Note that these coefficients satisfy our conditions for being consistent with a truncation error\n",
    "$$\n",
    "    \\tau^n = \\frac{1}{\\Delta t} (u(t_{n+2}) - 3 u(t_{n+1}) + 2 u(t_n) + \\Delta t u'(t_n)) = \\frac{5}{2} \\Delta t u''(t_n) + \\mathcal{O}(\\Delta t^2).\n",
    "$$\n",
    "It turns out that although this method is consistent the global error does not converge in general!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\n",
    "    U^{n+2} - 3 U^{n+1} + 2 U^n = - \\Delta t f(U^n)\n",
    "$$\n",
    "Consider the above method with the trivial ODE\n",
    "$$\n",
    "    u'(t) = 0 \\quad u(0) = 0\n",
    "$$\n",
    "so that we are left with the method\n",
    "$$\n",
    "    U^{n+2} - 3 U^{n+1} + 2 U^n = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we have exact values for $U^0$ and $U^1$ then this method would lead to $U^n = 0$. In general however we only have an approximation to $U^1$ so what happens then?  We can solve the linear difference equation in terms of $U^0$ and $U^1$ to find\n",
    "$$\n",
    "    U^n = 2 U^0 - U^1 + 2^n (U^1 - U^0).\n",
    "$$\n",
    "If we assume a error on the order of $\\mathcal{O}(\\Delta t)$ for $U^1$ this quickly leads to large values even for small $n$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Characteristic Polynomials and Linear Difference Equations\n",
    "\n",
    "As an short aside, say we wanted to solve\n",
    "\n",
    "$$\\sum^r_{j=0} \\alpha_j U^{n+j} = 0$$\n",
    "\n",
    "given initial conditions $U^0, U^1, \\ldots, U^{r-1}$ which has a solution in the general form $U^n = \\xi^n$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plugging this into the equation we have\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j \\xi^{n+j} = 0\n",
    "$$\n",
    "which simplifies to\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j \\xi^j = 0 \n",
    "$$\n",
    "\n",
    "by dividing by $\\xi^n$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $\\xi$ then is a root of the polynomial\n",
    "$$\n",
    "    \\rho(\\xi) = \\sum^r_{j=0} \\alpha_j \\xi^j\n",
    "$$\n",
    "then $\\xi$ solves the equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note that since these are linear methods that a linear combination of solutions is also a solution so the general form of a solution has the form\n",
    "$$\n",
    "    U^n = c_1 \\xi_1^n + c_2 \\xi_2^n + \\cdots + c_r \\xi^n_r.\n",
    "$$\n",
    "Given initial values for $U^0, U^1, \\ldots$ we can uniquely determine the $c_j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### General Zero-Stability Result for LMM\n",
    "\n",
    "An $r$-step LMM is *zero-stable* if the roots of the characteristic polynomial $\\rho(\\xi)$ satisfy\n",
    "$$\n",
    "    |\\xi_j| \\leq 1 \\quad \\quad \\text{for} \\quad j=1,2,3,\\ldots,r\n",
    "$$\n",
    "if $\\xi_j$ is not repeated and $|\\xi_j| < 1$ for repeated roots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Consider the linear multistep method\n",
    "$$\n",
    "    U^{n+2} - 2 U^{n+1} + U^n = \\frac{\\Delta t}{2} (f(U^{n+2}) - f(U^n)).\n",
    "$$\n",
    "Applying this to the ODE $u'(t) = 0$ leads to the difference equation\n",
    "$$\n",
    "    U^{n+2} - 2 U^{n+1} + U^n = 0\n",
    "$$\n",
    "whose characteristic polynomial is\n",
    "$$\n",
    "    \\rho(\\xi) = \\xi^2 - 2 \\xi + 1 = (\\xi - 1)^2\n",
    "$$\n",
    "leading to the general solution\n",
    "$$\n",
    "    U^n = c_1 + c_2 n.\n",
    "$$\n",
    "Here we see that given a $U^0$ and $U^1$ that the solution will still grow linearly with $n$ which will again lead to a divergent solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Consider the linear multistep method\n",
    "$$\n",
    "    U^{n+3} - 2 U^{n+2} + \\frac{5}{4} U^{n+1} - \\frac{1}{4} U^n = \\frac{\\Delta t}{4} f(U^n).\n",
    "$$\n",
    "Apply this to the ODE $u'(t) = 0$ and determine whether this method is zero-stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Applied to the ODE $u'(t) = 0$ we have the linear difference equation\n",
    "$$\n",
    "    U^{n+3} - 2 U^{n+2} + \\frac{5}{4} U^{n+1} - \\frac{1}{4} U^n = 0\n",
    "$$\n",
    "leading to \n",
    "$$\n",
    "    \\rho(\\xi) = \\xi^3 - 2 \\xi^2 + \\frac{5}{4} \\xi - \\frac{1}{4} = 0\n",
    "$$\n",
    "whose solutions are\n",
    "$$\n",
    "    \\xi_1 = 1, \\xi_2 = \\xi_3 = 1 / 2\n",
    "$$\n",
    "with the general solution\n",
    "$$\n",
    "    U^n = c_1 + c_2 \\frac{1}{2^n} + c_3 n \\frac{1}{2^n}\n",
    "$$\n",
    "which does converge due to the factor of $1 / 2^n$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: Adams Methods\n",
    "\n",
    "The general form for all Adams methods take the form\n",
    "$$\n",
    "    U^{n+r} = U^{n+r-1} + \\Delta t \\sum^r_{j=0} \\beta_j f(U^{n+j})\n",
    "$$\n",
    "which has the characteristic polynomial (for the ODE $u'(t)=0$)\n",
    "$$\n",
    "    \\rho(\\xi) = \\xi^r - \\xi^{r-1} = (\\xi - 1) \\xi^{r-1}\n",
    "$$\n",
    "leading to the roots $\\xi_1 = 1$ and $\\xi_2 = \\xi_3 = \\cdots = \\xi_r = 0$ which satisfy the general zero-stability result and therefore all Adams methods are convergent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Absolute Stability\n",
    "Although zero-stability guarantees stability it is much more difficult to work with in general as the limit $\\Delta t \\rightarrow 0$ can be difficult to compute.  Instead we often consider a finite $\\Delta t$ and examine if the method is stable for this particular choice of $\\Delta t$.  This has the practical upside that it will also tell us what particular $\\Delta t$ will ensure that our method is indeed stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "Consider the problem\n",
    "\n",
    "$$u'(t) = \\lambda (u - \\cos t) - \\sin t \\quad \\text{with} \\quad u(0) = 1$$\n",
    "\n",
    "whose exact solution is \n",
    "\n",
    "$$u(t) = \\cos t.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can compute an estimate for what $\\Delta t$ we need to use by examining the truncation error \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\tau &= \\frac{1}{2} \\Delta t u''(t) + \\mathcal{O}(\\Delta t^2) \\\\\n",
    "      &= -\\frac{1}{2} \\Delta t \\cos t + \\mathcal{O}(\\Delta t^2)\n",
    "\\end{aligned}$$\n",
    "\n",
    "and therefore\n",
    "\n",
    "$$|E^n| \\leq \\Delta t \\max_{0 \\leq t \\leq t_f} |\\cos t| = \\Delta t.$$\n",
    "\n",
    "If we want a solution where $|E^n| < 10^{-3}$ then $\\Delta t \\approx 10^{-3}$.  Turning to the application of Euler's method lets apply this to the case where $\\lambda = -10$ and $\\lambda = -2100$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler\n",
    "f = lambda t, lam, u: lam * (u - numpy.cos(t)) - numpy.sin(t)\n",
    "u_exact = lambda t: numpy.cos(t)\n",
    "\n",
    "t_f = 2.0\n",
    "num_steps = [2**n for n in range(4, 10)]\n",
    "# num_steps = [2**n for n in range(15,20)]\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_10 = numpy.empty(len(num_steps))\n",
    "error_2100 = numpy.empty(len(num_steps))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    \n",
    "    # Compute Euler solution\n",
    "    U = numpy.empty(t.shape)\n",
    "    U[0] = 1.0\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U[n+1] = U[n] + delta_t[i] * f(t_n, -10.0, U[n])\n",
    "    error_10[i] = numpy.abs(U[-1] - u_exact(t_f)) / numpy.abs(u_exact(t_f))\n",
    "    \n",
    "    U = numpy.empty(t.shape)\n",
    "    U[0] = 1.0\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U[n+1] = U[n] + delta_t[i] * f(t_n, -2100.0, U[n])\n",
    "    error_2100[i] = numpy.abs(U[-1] - u_exact(t_f)) / numpy.abs(u_exact(t_f))\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "\n",
    "axes.loglog(delta_t, error_10, 'bo', label='Forward Euler')\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_10[1], 1.0) * delta_t**1.0, 'r--', label=\"1st Order\")\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "\n",
    "axes.loglog(delta_t, error_2100, 'bo', label='Forward Euler')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_2100[1], 1.0) * delta_t**1.0, 'r--', label=\"1st Order\")\n",
    "\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So what went wrong with $\\lambda = -2100$?  The global error should go as\n",
    "\n",
    "$$E^{n+1} = (1 + \\Delta t \\lambda) E^n - \\Delta t T^n$$\n",
    "\n",
    "If $\\Delta t \\approx 10^{-3}$ then for the case $\\lambda = -10$ the previous global error is multiplied by\n",
    "\n",
    "$$1 + 10^{-3} \\cdot -10 = 0.99$$\n",
    "\n",
    "which means the contribution from $E^n$ will slowly decrease as we take more time steps.  For the other case we have\n",
    "\n",
    "$$1 + 10^{-3} \\cdot -2100 = -1.1$$\n",
    "\n",
    "which means that for this $\\Delta t$ the error made in previous time steps will grow!  For this not to happen we would have to have $\\Delta t < 1 / 2100$ which would lead to convergence again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute Stability of the Forward Euler Method\n",
    "\n",
    "Consider again the simple test problem $u'(t) = \\lambda u$.  We know from before that applying Euler's method to this problem leads to an update of the form\n",
    "\n",
    "$$U_{n+1} = (1 + \\Delta t \\lambda) U_n.$$\n",
    "\n",
    "As may have been clear from the last example, we know that if\n",
    "\n",
    "$$|1 + \\Delta t \\lambda| \\leq 1$$\n",
    "\n",
    "that the method will be stable, this is called **absolute stability**.  Note that the product of $\\Delta t \\lambda$ is what matters here and often we consider a **region of absolute stability** on the complex plain defined by the equation outlined where now $z = \\Delta t \\lambda$.  This allows the values of $\\lambda$ to be complex which can be an important case to consider, especially for systems of equations where the $\\lambda$s are identified as the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), 'b')\n",
    "axes.plot([-3, 3],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-3, 3],'k--')\n",
    "axes.set_xlim((-3, 3.0))\n",
    "axes.set_ylim((-3,3))\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Forward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### General Stability Regions for Linear Multistep Methods\n",
    "\n",
    "Going back to linear multistep methods and applying them in general to our test problem we have\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j U_{n+j} = \\Delta t \\sum^r_{j=0} \\beta_j \\lambda U_{n+j}\n",
    "$$\n",
    "which can be written as \n",
    "$$\n",
    "    \\sum^r_{j=0} (\\alpha_j - \\beta_j \\Delta t \\lambda) U_{n+j} = 0\n",
    "$$\n",
    "or using our notation of $z = \\Delta t \\lambda$ we have\n",
    "$$\n",
    "    \\sum^r_{j=0} (\\alpha_j - \\beta_j z) U_{n+j} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This has a similar form to the linear difference equations considered above!  Letting\n",
    "$$\n",
    "    \\rho(\\xi) = \\sum^r_{j=0} \\alpha_j \\xi^j\n",
    "$$\n",
    "and \n",
    "$$\n",
    "    \\sigma(\\xi) = \\sum^r_{j=0} \\beta_j \\xi^j\n",
    "$$\n",
    "we can write the expression above as\n",
    "$$\n",
    "    \\pi(\\xi, z) = \\rho(\\xi) - z \\sigma(\\xi)\n",
    "$$\n",
    "called the **stability polynomial** of the the linear multi-step method.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out that if the roots $\\xi_i$ of this polynomial satisfy\n",
    "$$\n",
    "    |\\xi_i| \\leq 1\n",
    "$$\n",
    "then the multi-step method is absolutely-stable.  We then define the region of absolute stability as the values for $z$ for which this is true.  This approach can also be applied to one-step methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Forward Euler's Method\n",
    "\n",
    "Examining forward Euler's method we have\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    0 &= U_{n+1} - U_n - \\Delta t \\lambda U_n \\\\\n",
    "    &= U_{n+1} - U_n (1 + \\Delta t \\lambda)\\\\\n",
    "    &= \\xi - 1 (1 + z)\\\\\n",
    "    &=\\pi(\\xi, z)\n",
    "\\end{aligned}$$\n",
    "\n",
    "whose root is $\\xi = 1 + z$ and we have re-derived the stability region we had found before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute Stability of the backward Euler Method\n",
    "\n",
    "The backward version of Euler's method is defined as\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t f(t_{n+1}, U_{n+1}).\n",
    "$$\n",
    "Check to see if backward Euler is absolute stable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we again consider the test problem from before we find that\n",
    "$$\\begin{aligned}\n",
    "    0 &= U_{n+1} (1 - \\Delta t \\lambda) - U_n \\\\\n",
    "    &= \\xi (1 - z) - 1\n",
    "\\end{aligned}$$\n",
    "which has the root $\\xi = \\frac{1}{1 - z}$.  We then have\n",
    "$$\n",
    "    \\left|\\frac{1}{1-z}\\right| \\leq 1 \\leftrightarrow |1 - z| \\geq 1\n",
    "$$\n",
    "so in fact the stability region encompasses the entire complex plane except for a circle centered at $(1, 0)$ of radius 1 implying that the backward Euler method is in fact stable for any choice of $\\Delta t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Application to Stiff ODEs\n",
    "\n",
    "Consider again the ODE we examined before\n",
    "\n",
    "$$u'(t) = \\lambda (u - \\cos t) - \\sin t$$\n",
    "\n",
    "except this time with general initial condition $u(t_0) = \\eta$.  What happens to solutions that are slightly different from $\\eta = 1$ or $t_0 = 0$?  The general solution of the ODE is\n",
    "\n",
    "$$u(t) = e^{\\lambda (t - t_0)} (\\eta - \\cos t_0)) + \\cos t$$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot \"hairy\" solutions to the ODE\n",
    "u = lambda t_0, eta, lam, t: numpy.exp(lam * (t - t_0)) * (eta - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "for lam in [-1, -10]:\n",
    "    fig = plt.figure()\n",
    "    axes = fig.add_subplot(1, 1, 1)\n",
    "    for eta in numpy.linspace(-1, 1, 10):\n",
    "        for t_0 in numpy.linspace(0.0, 9.0, 10):\n",
    "            t = numpy.linspace(t_0,10.0,100)\n",
    "            axes.plot(t, u(t_0, eta, lam, t),'b')\n",
    "    t = numpy.linspace(0.0,10.0,100)\n",
    "    axes.plot(t, numpy.cos(t), 'r', linewidth=5)\n",
    "        \n",
    "    axes.set_title(\"Perturbed Solutions $\\lambda = %s$\" % lam)\n",
    "    axes.set_xlabel('$t$')\n",
    "    axes.set_ylabel('$u(t)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot \"inverse hairy\" solutions to the ODE\n",
    "u = lambda t_0, eta, lam, t: numpy.exp(lam * (t - t_0)) * (eta - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "fig = plt.figure()\n",
    "num_steps = 10\n",
    "error = numpy.ones(num_steps) * 1.0\n",
    "t_hat = numpy.linspace(0.0, 10.0, num_steps + 1)\n",
    "t_whole = numpy.linspace(0.0, 10.0, 1000)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "eta = 1.0\n",
    "lam = 0.1\n",
    "\n",
    "for n in range(1,num_steps):\n",
    "    t = numpy.linspace(t_hat[n-1], t_hat[n], 100)\n",
    "    U = u(t_hat[n-1], eta, lam, t)\n",
    "    axes.plot(t, U, 'b')\n",
    "    axes.plot(t_whole, u(t_hat[n-1], eta, lam, t_whole),'b--')\n",
    "    axes.plot([t[-1], t[-1]], (U[-1], U[-1] + -1.0**n * error[n]), 'r')\n",
    "    eta = U[-1] + -1.0**n * error[n]\n",
    "\n",
    "t = numpy.linspace(0.0, 10.0, 100)\n",
    "axes.plot(t, numpy.cos(t), 'g')\n",
    "\n",
    "axes.set_title(\"Perturbed Solutions $\\lambda = %s$\" % lam)\n",
    "axes.set_xlabel('$t$')\n",
    "axes.set_ylabel('$u(t)$')\n",
    "axes.set_ylim((-10,10))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Chemical systems\n",
    "\n",
    "Consider the transition of a chemical $A$ to a chemical $C$ through the process\n",
    "\n",
    "$$A \\overset{K_1}{\\rightarrow} B \\overset{K_2}{\\rightarrow} C.$$\n",
    "\n",
    "If we let\n",
    "\n",
    "$$\\vec{u} = \\begin{bmatrix} [A] \\\\ [B] \\\\ [C] \\end{bmatrix}$$\n",
    "\n",
    "then we can model this simple chemical reaction with the system of ODEs\n",
    "\n",
    "$$\\frac{\\text{d} \\vec{u}}{\\text{d} t} = \n",
    "\\begin{bmatrix}\n",
    "    -K_1 & 0 & 0 \\\\\n",
    "    K_1 & -K_2 & 0 \\\\\n",
    "    0 & K_2 & 0\n",
    "\\end{bmatrix} \\vec{u}$$\n",
    "\n",
    "The solution of this system is of the form\n",
    "\n",
    "$$u_j(t) = c_{j1} e^{-K_1 t} + c_{j2}e^{-K_2 t} + c_{j3}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Solve the chemical systems example\n",
    "# Problem parameters\n",
    "K_1 = 3\n",
    "K_2 = 1\n",
    "# K_1 = 30.0\n",
    "# K_2 = 1.0\n",
    "\n",
    "A = numpy.array([[-K_1, 0, 0], [K_1, -K_2, 0], [0, K_2, 0]])\n",
    "f = lambda u: numpy.dot(A, u)\n",
    "\n",
    "t = numpy.linspace(0.0, 8.0, 128)\n",
    "delta_t = t[1] - t[0]\n",
    "\n",
    "U = numpy.empty((t.shape[0], 3))\n",
    "U[0, :] = [2.5, 5.0, 2.0]\n",
    "\n",
    "for n in range(t.shape[0] - 1):\n",
    "    U[n+1, :] = U[n, :] + delta_t * f(U[n, :])\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t, U)\n",
    "\n",
    "axes.set_title(\"Chemical System\")\n",
    "axes.set_xlabel(\"$t$\")\n",
    "axes.set_title(\"$[A], [B], [C]$\")\n",
    "axes.set_ylim((0.0, 10.))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### What is stiffness?\n",
    "\n",
    "In general a **stiff** ODE is one where $u'(t) \\ll f'(t, u)$.  For systems of ODEs the **stiffness ratio**\n",
    "\n",
    "$$\\frac{\\max_p |\\lambda_p|}{\\min_p |\\lambda_p|}$$\n",
    "\n",
    "can be used to characterize the stiffness of the system.  In our last example this ratio was $K_1 / K_2$ if $K_1 > K_2$.  As we increased this ratio we observed that the numerical method became unstable only a reduction in $\\Delta t$ lead to stable solution again.  For explicit time step methods this is problematic as the reduction of the time step for only one of the species leads to very expensive evaluations.  For example, forward Euler has the stability criteria\n",
    "\n",
    "$$|1 + \\Delta t \\lambda| < 1$$\n",
    "\n",
    "where $\\lambda$ will have to be the maximum eigenvalue of the system.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "K_1 = 3.0\n",
    "K_2 = 1.0\n",
    "delta_t = 1.0\n",
    "eigenvalues = [-K_1, -K_2]\n",
    "\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), color='lightgray')\n",
    "for lam in eigenvalues:\n",
    "    print(lam * delta_t)\n",
    "    axes.plot(lam * delta_t, 0.0, 'ko')\n",
    "axes.plot([-3, 3],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-3, 3],'k--')\n",
    "# axes.set_xlim((-3, 1))\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Forward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A-Stability\n",
    "What if we could expand the absolute stability region to encompass more of the left-half plane or even better, all of it.  A method that has this property is called **A-stable**.  We have already seen one example of this with backward Euler which has a stability region of\n",
    "\n",
    "$$|1 - z| \\geq 1$$\n",
    "\n",
    "which covers the full left-half plane.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out that for linear multi-step methods a theorem by Dahlquist proves that there are no LMMs that satisfies the A-stability criterion beyond second order (trapezoidal rule).  There are higher-order Runge-Kutta methods do however."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Perhaps this is too restrictive though.  Often large eigenvalues for systems (for instance coming from a PDE discretization for the heat equation) lie completely on the real line.  If the stability region can encompass as much of the real line as possible while leaving out the rest of the left-half plane we can possibly get a more efficient method.  There are a number of methods that can be constructed that have this property but are higher-order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Backward Euler\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "K_1 = 3.0\n",
    "K_2 = 1.0\n",
    "delta_t = 1.0\n",
    "eigenvalues = [-K_1, -K_2]\n",
    "\n",
    "axes.set_facecolor('lightgray')\n",
    "axes.fill(numpy.cos(t) + 1.0, numpy.sin(t), 'w')\n",
    "for lam in eigenvalues:\n",
    "    print(lam * delta_t)\n",
    "    axes.plot(lam * delta_t, 0.0, 'ko')\n",
    "axes.plot([-3, 3],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-3, 3],'k--')\n",
    "# axes.set_xlim((-3, 1))\n",
    "axes.set_ylim((-2,2))\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Backward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### L-Stability\n",
    "\n",
    "It turns out not all A-stable methods are alike.  Consider the backward Euler method and the trapezoidal method.  The stability polynomial for the trapezoidal method is\n",
    "$$\\begin{aligned}\n",
    "    0 &= U_{n+1} - U_n - \\Delta t \\frac{1}{2} (\\lambda U_n + \\lambda U_{n+1}) \\\\\n",
    "      &= U_{n+1}\\left(1 - \\frac{1}{2} \\Delta t \\lambda \\right ) - U_n \\left(1 + \\frac{1}{2}\\Delta t \\lambda \\right) \\\\\n",
    "      &= \\left(\\xi - \\frac{1 + \\frac{1}{2}z}{1 - \\frac{1}{2} z}\\right) \\left(1 - \\frac{1}{2} z \\right )\\\\\n",
    "\\end{aligned}$$\n",
    "which shows that it is A-stable.  Lets apply both these methods to a problem we have seen before and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler\n",
    "f = lambda t, lam, u: lam * (u - numpy.cos(t)) - numpy.sin(t)\n",
    "u_exact = lambda t_0, eta, lam, t: numpy.exp(lam * (t - t_0)) * (eta - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "t_0 = 0.0\n",
    "eta = 1.5\n",
    "lam = -1e6\n",
    "\n",
    "num_steps = [10, 20, 40, 50]\n",
    "\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_euler = numpy.empty(len(num_steps))\n",
    "error_trap = numpy.empty(len(num_steps))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    u = u_exact(t_0, eta, lam, t_f)\n",
    "    \n",
    "    # Compute Euler solution\n",
    "    U_euler = numpy.empty(t.shape)\n",
    "    U_euler[0] = eta\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U_euler[n+1] = (U_euler[n] - lam * delta_t[i] * numpy.cos(t_n) - delta_t[i] * numpy.sin(t_n)) / (1.0 - lam * delta_t[i])\n",
    "    error_euler[i] = numpy.abs(U_euler[-1] - u) / numpy.abs(u)\n",
    "    \n",
    "    # Compute using trapezoidal\n",
    "    U_trap = numpy.empty(t.shape)\n",
    "    U_trap[0] = eta\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U_trap[n+1] = (U_trap[n] + delta_t[i] * 0.5 * f(t_n, lam, U_trap[n]) - 0.5 * lam * delta_t[i] * numpy.cos(t_n) - 0.5 * delta_t[i] * numpy.sin(t_n)) / (1.0 - 0.5 * lam * delta_t[i])\n",
    "    error_trap[i] = numpy.abs(U_trap[-1] - u)\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t, U_euler, 'ro-')\n",
    "axes.plot(t, u_exact(t_0, eta, lam, t),'k')\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.loglog(delta_t, error_euler, 'bo')\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_euler[1], 1.0) * delta_t**1.0, 'r--', label=\"1st Order\")\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_euler[1], 2.0) * delta_t**2.0, 'b--', label=\"2nd Order\")\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors for Backwards Euler\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "# Plots for trapezoid\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t, U_trap, 'ro-')\n",
    "axes.plot(t, u_exact(t_0, eta, lam, t),'k')\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.loglog(delta_t, error_trap, 'bo', label='Forward Euler')\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_trap[1], 1.0) * delta_t**1.0, 'r--', label=\"1st Order\")\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_trap[1], 2.0) * delta_t**2.0, 'b--', label=\"2nd Order\")\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors for Trapezoidal Rule\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It turns out that if we look at a one-step method and define the following ratio\n",
    "\n",
    "$$U_{n+1} = R(z) U_n$$\n",
    "\n",
    "we can define another form of stability, called **L-stability**, where we require that the method is A-stable and that\n",
    "\n",
    "$$\\lim_{z \\rightarrow \\infty} |R(z)| = 0.$$\n",
    "\n",
    "Backwards Euler is L-stable while the trapezoidal method is not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Backward Differencing Formulas\n",
    "\n",
    "A class of LMM methods that are useful for stiff ODE problems are the backward difference formula (BDF) methods which have the form\n",
    "\n",
    "$$\\alpha_0 U_n + \\alpha_1 U_{n+1} + \\cdots + \\alpha_r U_{n+r} = \\Delta \\beta_r f(U_{n+r})$$\n",
    "\n",
    "These methods can be derived directly from backwards finite differences from the point $U_{n+r}$ and the rest of the points back in time.  One can then derive r-step methods that are rth-order accurate this way.  Some of the methods are \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    r = 1:& & U_{n+1} - U_n = \\Delta t f(U_{n+1}) \\\\\n",
    "    r = 2:& &3 U_{n+2} - 4 U_{n+1} + U_n = 2 \\Delta t f(U_{n+1}) \\\\\n",
    "    r = 3:& &11U_{n+3} - 18U_{n+2} + 9U_{n+1} - 2 U_n = 6 \\Delta t f(U_{n+3}) \\\\\n",
    "    r = 4:& &25 U_{n+4} - 48 U_{n+3} +36 U_{n+2} -16 U_{n+1} +3 U_n = 12 \\Delta t f(U_{n+4})\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plotting Stability Regions\n",
    "\n",
    "If we think of the roots of the stability polynomial $\\xi_j$ as complex numbers and write them in exponential form \n",
    "$$\n",
    "    \\xi_j = |\\xi_j| e^{i \\theta}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Here $|\\xi_j|$ is the modulus (or magnitude) or the complex number and is defined as $|\\xi_j| = x^2 + y^2$ where $\\xi_j = x + i j$.  If the $\\xi_j$s are on the boundary of the absolute stability region then we know that $|\\xi_j| = 1$.  Using this in conjunction with the stability polynomial then leads to \n",
    "$$\n",
    "    \\rho(e^{i\\theta}) - z \\sigma(e^{i\\theta}) = 0\n",
    "$$\n",
    "which solving for $z$ leads to\n",
    "$$\n",
    "    z(\\theta) = \\frac{\\rho(e^{i\\theta})}{\\sigma(e^{i\\theta})}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As an example consider the Adams-Bashforth 2-stage method.  The stability polynomial can be found as\n",
    "$$\\begin{aligned}\n",
    "    U_{n+2} &= U_{n+1} + \\frac{\\Delta t}{2} (-f(U_n) + 3 f(U_{n+1})) \\\\\n",
    "    \\pi(\\xi, z) &= U_{n+2} - U_{n+1} - \\frac{\\Delta t}{2} (-f(U_n) + 3 f(U_{n+1})) = 0 \\\\\n",
    "    &= U_{n+2} - U_{n+1} - \\frac{1}{2} (\\Delta t \\lambda U_n - 3 \\Delta t \\lambda U_{n+1}) \\\\\n",
    "    &= 2 \\xi^2 - 2 \\xi + 3 z\\xi - z \\\\\n",
    "    &= \\rho(\\xi, z) + z \\sigma(\\xi, z)\n",
    "\\end{aligned}$$\n",
    "where\n",
    "$$\n",
    "    \\rho(\\xi, z) = 2 ( \\xi - 1) \\xi ~~~ \\text{and} ~~~ \\sigma(\\xi, z) = 3 \\xi - 1\n",
    "$$\n",
    "so that\n",
    "$$\n",
    "    z(\\theta) = \\frac{2 (\\xi - 1) \\xi}{3 \\xi - 1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This does not necessarily ensure that given a $\\theta$ that $z(\\theta)$ will lie on the absolute stability region's boundary.  This can occur when $\\xi_j = 1$ but to the left and right of the curve $\\xi_j > 1$ and so therefore does not mark the boundary of the region.  To determine whether a particular region outlined by this curve is inside or outside of the stability region we can evaluate all the roots of $\\pi(\\xi, z)$ at some $z$ inside of the region in question and see if then $\\forall j, \\xi_j < 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For one-step methods this becomes easier, if we look at the ratio $R(z)$ we defined earlier as\n",
    "$$\n",
    "    U_{n+1} = R(z) U_n\n",
    "$$\n",
    "in the case of the pth order Taylor series method applied to $u'(t) = \\lambda u$ we get\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    U_{n+1} &= U_n + \\Delta t \\lambda U_n + \\frac{1}{2}\\Delta t^2 \\lambda^2 U_n + \\cdots + \\frac{1}{p!}\\Delta t^p \\lambda^p U_n \\\\\n",
    "    &=\\left(1 + z + \\frac{1}{2} z^2 + \\frac{1}{6} z^3 + \\cdots +\\frac{1}{p!}z^p\\right) U_n \\Rightarrow \\\\\n",
    "    R(z) &= 1 + z + \\frac{1}{2} z^2 + \\frac{1}{6} z^3 + \\cdots +\\frac{1}{p!}z^p.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Setting $R(z) = e^{i\\theta}$ could lead to a way for solving for the boundary but (where $|R(z)| = 1$) but this is very difficult to do in general.  Instead if we plot the contours of $R(z)$ in the complex plain we can pick out the $R(z)=1$ contour and plot that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "theta = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "# ==================================\n",
    "#  Forward euler\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2.0)\n",
    "fig.set_figheight(fig.get_figheight() * 2.0)\n",
    "axes = fig.add_subplot(2, 2, 1)\n",
    "x = numpy.linspace(-2, 2, 100)\n",
    "y = numpy.linspace(-2, 2, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "Z = X + 1j * Y\n",
    "\n",
    "# Evaluate which regions are actually in stability region\n",
    "print(\"Forward Euler's method - Stability region tests\")\n",
    "z = -1 + 1j * 0\n",
    "print(\"  Inside of circle: \", numpy.abs(1.0 + z))\n",
    "z = -3 + 1j * 0\n",
    "print(\"  Outside of circle: \", numpy.abs(1.0 + z))\n",
    "\n",
    "axes.contour(X, Y, numpy.abs(1.0 + Z), levels=[1.0])\n",
    "axes.plot(x, numpy.zeros(x.shape),'k')\n",
    "axes.plot(numpy.zeros(y.shape), y,'k')\n",
    "axes.set_aspect('equal')\n",
    "axes.set_title(\"Forward Euler\")\n",
    "\n",
    "# ==================================\n",
    "#  Backwards Euler\n",
    "axes = fig.add_subplot(2, 2, 2)\n",
    "x = numpy.linspace(-2, 2, 100)\n",
    "y = numpy.linspace(-2, 2, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "Z = X + 1j * Y\n",
    "\n",
    "# Evaluate which regions are actually in stability region\n",
    "print(\"Backward Euler's method - Stability region tests\")\n",
    "z = 1 + 1j * 0\n",
    "print(\"  Inside of circle: \", numpy.abs(1.0 + z))\n",
    "z = -1 + 1j * 0\n",
    "print(\"  Outside of circle: \", numpy.abs(1.0 + z))\n",
    "\n",
    "axes.contour(X, Y, numpy.abs(1.0 / (1.0 - Z)), levels=[1.0])\n",
    "axes.plot(x, numpy.zeros(x.shape),'k')\n",
    "axes.plot(numpy.zeros(y.shape), y,'k')\n",
    "axes.set_aspect('equal')\n",
    "axes.set_title(\"Backwards Euler\")\n",
    "\n",
    "# ==================================\n",
    "#  Taylor series method of order 4\n",
    "axes = fig.add_subplot(2, 2, 3)\n",
    "x = numpy.linspace(-5, 5, 100)\n",
    "y = numpy.linspace(-5, 5, 100)\n",
    "X, Y = numpy.meshgrid(x, y)\n",
    "Z = X + 1j * Y\n",
    "\n",
    "# Evaluate which regions are actually in stability region\n",
    "print(\"Taylor series method of order 4 - Stability region tests\")\n",
    "z = -1 + 1j * 0\n",
    "print(\"  Inside of strange region: \", numpy.abs(1.0 + z))\n",
    "z = -4 + 1j * 0\n",
    "print(\"  Outside of strange region: \", numpy.abs(1.0 + z))\n",
    "\n",
    "axes.contour(X, Y, numpy.abs(1 + Z + 0.5 * Z**2 + 1.0/6.0 * Z**3 + 1.0 / 24.0 * Z**4), levels=[1.0])\n",
    "axes.plot(x, numpy.zeros(x.shape),'k')\n",
    "axes.plot(numpy.zeros(y.shape), y,'k')\n",
    "axes.set_aspect('equal')\n",
    "axes.set_title(\"4th Order Taylor Series\")\n",
    "\n",
    "# ==================================\n",
    "# 2-step Adams-Bashforth\n",
    "theta = numpy.linspace(0.0, 2.0 * numpy.pi, 1000)\n",
    "xi = numpy.exp(1j * theta)\n",
    "\n",
    "rho_2AB = lambda xi: (xi - 1.0) * xi\n",
    "sigma_2AB = lambda xi: (3.0 * xi - 1.0) / 2.0\n",
    "z_2AB = rho_2AB(xi) / sigma_2AB(xi)\n",
    "z = rho_2AB(xi) / sigma_2AB(xi)\n",
    "\n",
    "axes = fig.add_subplot(2, 2, 4)\n",
    "axes.plot(z_2AB.real, z_2AB.imag)\n",
    "axes.plot(x, numpy.zeros(x.shape),'k')\n",
    "axes.plot(numpy.zeros(y.shape), y,'k')\n",
    "axes.set_title(\"2-step Adams-Bashforth\")\n",
    "axes.set_aspect('equal')\n",
    "axes.set_xlim([-2, 3])\n",
    "axes.set_ylim([-2, 2])\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "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.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}