09_ODE_ivp_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.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>\n",
    "\n",
    "Note:  The presentation below largely follows part II in \"Finite Difference Methods for Ordinary and Partial Differential Equations\" by LeVeque (SIAM, 2007)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Solution to ODE Initial Value Problems - Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Overview -- so far\n",
    "\n",
    "So far we have discussed 3 basic techniques for integration of ODE IVP's\n",
    "\n",
    "* **Single-Step Multi-Stage** schemes (explicit and implicit)\n",
    "* **Taylor's Series** Methods\n",
    "* **Linear Multi-Step** Schemes (explicit and implicit)\n",
    "\n",
    "as well as \n",
    "* **truncation error** of each method (and its relation to step-error)\n",
    "* **adaptive stepping** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Overview -- to do\n",
    "\n",
    "This notebook will continue the analysis of these methods to discuss the other important issues of\n",
    "* **Convergence** (Global Error) -- does the discrete solution converge to the continuous analytic solution?\n",
    "* **Stability** -- Are there limits to how large a time-step we can take?\n",
    "* **Stiff Equations**:  consequences of stability\n",
    "* **Implicit Methods**: Solution of non-linear systems of equation using **Newton's Method**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Convergence\n",
    "\n",
    "We can think of an ODE method as providing a sequence of approximations $U_N$ where here $N$ is the number of time steps needed to reach the final time of interest $t_f$, in effect we are increasing the resolution of the method.  We then say that a method is convergent if this sequence converges to the true solution at the same time\n",
    "$$\n",
    "    \\lim_{N\\rightarrow \\infty} U_N = u(t_f).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can also define this in a more familiar way in terms of $\\Delta t$ such that\n",
    "$$\n",
    "    N \\Delta t = t_f ~~~~ \\Rightarrow ~~~~ N = \\frac{t_f}{\\Delta t}\n",
    "$$\n",
    "so that our definition of convergence becomes\n",
    "$$\n",
    "    \\lim_{\\Delta t \\rightarrow 0} U_{\\Delta t} = u(t_f).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In general for a method to be convergent it must be\n",
    " - **consistent** which as before meant that the local truncation error $T = \\mathcal{O}(\\Delta t^p)$ where $p > 0$,\n",
    " - **zero-stable** which implies that the sum total of the errors as $\\Delta t \\rightarrow 0$ is bounded and has the same order as $T$ (the truncation error) which we know goes to zero as $\\Delta t \\rightarrow 0$.\n",
    " \n",
    " \n",
    "put another way,  a method is convergent if the **global error**\n",
    "\n",
    "$$\n",
    "    E = | U_N(t_f) - u(t_f)| \\rightarrow 0\\quad\\mathrm{as}\\quad N\\rightarrow\\infty\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Forward Euler on a Linear Problem\n",
    "\n",
    "Consider the simple linear problem\n",
    "\n",
    "$$\n",
    "    \\frac{\\text{d}u}{\\text{d}t} = \\lambda u \\quad \\text{with} \\quad u(0) = u_0\n",
    "$$\n",
    "\n",
    "which we know has the solution $u(t) = u_0 e^{\\lambda t}$.  Applying Euler's method to this problem leads to \n",
    "\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 defined by\n",
    "$$\\begin{aligned}\n",
    "    T_n &= \\frac{1}{\\Delta t} \\left[ U_{n+1} - u_{n+1} \\right ] \\\\\n",
    "        &= \\frac{1}{\\Delta t} \\left[ (1 + \\lambda \\Delta t) u_n - u_{n+1} \\right ]\n",
    "\\end{aligned}$$\n",
    "which can be rearranged to find\n",
    "\n",
    "$$\n",
    "    u_{n+1} = (1 + \\Delta t \\lambda) u_n - \\Delta t T_n.\n",
    "$$\n",
    "\n",
    "Note that all values here are in terms of the exact solution where as the application of Euler's method is in terms of the approximate solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Defining the global error as\n",
    "$$\n",
    "    E_{n} = u_n - U_n\n",
    "$$\n",
    "\n",
    "we can subtract the last two expressions to find\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    u_{n+1} - U_{n+1} &= \\underbrace{(1 + \\Delta t \\lambda) u_n - \\Delta t T_n}_{u_{n+1}} - \\underbrace{(1 + \\Delta t \\lambda) U_n}_{U_{n+1}} \\\\\n",
    "    E_{n+1} &= (1 + \\Delta t \\lambda) (u_n - U_n) - \\Delta t T_n \\\\\n",
    "            &= (1 + \\Delta t \\lambda) E_n - \\Delta t T_n,\n",
    "\\end{aligned}$$\n",
    "\n",
    "a recursive definition for the global error that connects it to the local truncation error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Working backwards from the global error $E_n$ at maximum time step $n = t_f/\\Delta t$,  the first few iterates of this look like\n",
    "\n",
    "$$\\begin{align}\n",
    "    E_n &= (1 + \\lambda \\Delta t) E_{n - 1} - \\Delta t T_{n-1} \\\\\n",
    "    &= (1 + \\lambda \\Delta t) \\left [(1 + \\lambda \\Delta t) E_{n - 2} - \\Delta t T_{n-2} \\right ] - \\Delta t T_{n-1} \\\\\n",
    "    &= (1 + \\lambda \\Delta t) \\left [(1 + \\lambda \\Delta t) \\left \\{(1 + \\lambda \\Delta t) E_{n - 3} - \\Delta t T_{n-3} \\right \\} - \\Delta t T_{n-2} \\right ] - \\Delta t T_{n-1}\n",
    "\\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which at this point looks like\n",
    "$$\n",
    "    E_n = (1 + \\lambda \\Delta t)^3E_{n-3} - \\Delta t\\sum_{i=1}^3 (1+\\Delta t\\lambda)^{i-1}T_{n-i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So taking it back $n$ steps ($3\\rightarrow n$) gives\n",
    "$$\n",
    "    E_n = (1 + \\lambda \\Delta t)^nE_0 - \\Delta t\\sum_{i=1}^n (1+\\Delta t\\lambda)^{i-1}T_{n-i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Expanding this expression out backwards in time to $n=0$ leads to \n",
    "$$\n",
    "    E_n = (1 + \\Delta t \\lambda)^n E_0 - \\Delta t \\sum^n_{i=1} (1 + \\Delta t \\lambda)^{i-1} T_{n - i}.\n",
    "$$\n",
    "\n",
    "Noting that $(1 + \\Delta t \\lambda)$ are the first two terms in the Taylor series for $e^{\\lambda\\Delta t}$ we will bound this term by\n",
    "\n",
    "$$\n",
    "    e^{\\Delta t \\lambda} \\geq |1 + \\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 - 1} \\leq e^{(n-1) \\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",
    "$$\\begin{aligned}\n",
    "    E_n &= (1 + \\Delta t \\lambda)^n E_0 - \\Delta t \\sum^n_{i=1} (1 + \\Delta t \\lambda)^{i-1} T_{n-i} \\\\\n",
    "    |E_n| &\\leq e^{|\\lambda| \\Delta t n} |E_0| + \\Delta t \\sum^n_{i=1} e^{|\\lambda| t_f} |T_{n-i}| \\\\\n",
    "          &\\leq e^{|\\lambda| t_f} \\left(|E_0| + \\Delta t \\sum^n_{i=1} |T_{n - i}|\\right) \\\\\n",
    "          &\\leq e^{|\\lambda| t_f} \\left(|E_0| + n \\Delta t \\max_{1 \\leq i \\leq n} |T_{n - i}|\\right)\n",
    "\\end{aligned}$$\n",
    "or"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "E_n &\\leq e^{|\\lambda| t_f} \\left(|E_0| + t_f \\max_{1 \\leq i \\leq n} |T_{n - i}|\\right)\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "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 truncation 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",
    "$$\n",
    "\n",
    "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 global error $E_n$ is bounded by $\\mathcal{O}(\\Delta t)$, the same as the local truncation error:\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": [
    "More generally,  it can be shown that for all single-step multi-stage schemes (of which Euler is the simplest), the Global error has the same order as the Truncation error, all of which scale as \n",
    "\n",
    "$$\n",
    "    \\Delta t^p,\\quad p>0\n",
    "$$\n",
    "\n",
    "and therefore all of these schemes converge."
   ]
  },
  {
   "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 (and in general isn't ideal as it requires an infinite amount of computation).  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": "skip"
    }
   },
   "source": [
    "### Example\n",
    "Consider the problem\n",
    "$$\n",
    "    u'(t) = \\lambda (u - \\cos t) - \\sin t ~~~~ \\text{with} ~~~~ u(0) = 1\n",
    "$$\n",
    "\n",
    "whose exact solution is \n",
    "$$\n",
    "    u(t) = \\cos t.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can compute an estimate for what $\\Delta t$ we need to use by examining the truncation error for Euler's method\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    T &= \\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",
    "\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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Implement vectorized Forward Euler  \n",
    "def euler(f, t_span, u0, N):\n",
    "    \"\"\" simple implementation of constant step-size forward euler method for vector valued f\n",
    "        This doc string should have so much more in it\n",
    "    \"\"\"\n",
    "    t = numpy.linspace(t_span[0], t_span[1],N)\n",
    "    if numpy.isscalar(u0):\n",
    "        u = numpy.empty(t.shape)\n",
    "    else:\n",
    "        u = numpy.empty((len(t),len(u0)))\n",
    "    u[0] = u0\n",
    "    delta_t = t[1] - t[0]\n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        K1 = delta_t * f(t_n, u[n])\n",
    "        u[n + 1] = u[n] + K1  \n",
    "    return t, u.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler for different values of lambda\n",
    "f = lambda t, u, lam: lam * (u - numpy.cos(t)) - numpy.sin(t)\n",
    "u_exact = lambda t: numpy.cos(t)\n",
    "u_0 = u_exact(0.)\n",
    "\n",
    "t_span = [0., 2.]\n",
    "num_steps = [2**n for n in range(4, 9)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#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",
    "t_f = t_span[1]\n",
    "u_f = u_exact(t_f)\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    # Compute Euler solution\n",
    "    ff = lambda t, u: f(t, u, -10.)\n",
    "    t, U_10 = euler(ff, t_span, u_0, N) \n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    error_10[i] = numpy.abs(U_10[-1] - u_f) / numpy.abs(u_f)\n",
    "    ff = lambda t, u: f(t, u, -2100.)    \n",
    "    t, U_2100 = euler(ff, t_span, u_0, N) \n",
    "    error_2100[i] = numpy.abs(U_2100[-1] - u_f) / numpy.abs(u_f)\n",
    "\n",
    "print('Error: lambda =   10: {}'.format(error_10))\n",
    "print('Error: lambda = 2100: {}'.format(error_2100))    \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(\"Error - $\\lambda = 10$\")\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(\"Error - $\\lambda = 2100$\")\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": "skip"
    }
   },
   "source": [
    "So what went wrong with $\\lambda = -2100$, the global error should go as\n",
    "\n",
    "$$\n",
    "    E_n = \\mathcal{O}(\\Delta t)?\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If $\\Delta t \\approx 10^{-3}$ then for the case $\\lambda = -10$ the previous **global error** is multiplied by $1+\\lambda\\Delta t$\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": [
    "### A simple example\n",
    "Consider our simplest test problem\n",
    "$$\n",
    "    u'(t) = \\lambda u, \\quad u(0) = u_0\n",
    "$$\n",
    "\n",
    "whose exact solution is \n",
    "$$\n",
    "    u(t) = u_0 e^{\\lambda t}.\n",
    "$$\n",
    "which for real $\\lambda < 0$ should decay to zero as $t\\rightarrow\\infty$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Question?\n",
    "\n",
    "What is the **largest** step-size $\\Delta t$ for an Euler method, such that the discrete solution $U_N\\rightarrow 0$ as $N\\rightarrow\\infty$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Euler's method on this simplest problem can be written\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    U_{n+1} &= U_n + \\Delta t f(U_n) = U_n + \\Delta t \\lambda U_n\\\\\n",
    "            &= (1 + \\lambda\\Delta t)U_n\\\\\n",
    "\\end{align}\n",
    "$$            \n",
    "\n",
    "i.e. at each step the last solution is multiplied by the factor $(1 + \\lambda\\Delta t)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "After $N$ steps\n",
    "\n",
    "$$\n",
    "    U_N = (1 + \\lambda\\Delta t)^N U_0\n",
    "$$ \n",
    "\n",
    "Assuming $\\lambda < 0$ what is the largest step size $\\Delta t$ such that $U_N\\rightarrow 0$ as $N\\rightarrow\\infty$?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Hopefully it's clear that for this scheme to be stable requires\n",
    "\n",
    "$$\n",
    "    |1 + \\lambda\\Delta t| < 1\n",
    "$$\n",
    "\n",
    "The key parameter is actually the combination of $\\lambda$ and the time step $z=\\lambda\\Delta t$, so stability requires that\n",
    "\n",
    "$$\n",
    "    | 1 + z | < 1\n",
    "$$\n",
    "or  𝑧∈[?,?]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6))\n",
    "z = numpy.linspace(-3,2,200)\n",
    "plt.plot(z,numpy.abs(1 + z))\n",
    "plt.plot(z,numpy.ones(z.shape),'k--')\n",
    "plt.plot([0.,0],[0.,3.])\n",
    "plt.xlabel('$z$',fontsize=16)\n",
    "plt.ylabel('$|1 + z|$',fontsize=16)\n",
    "plt.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Given $z\\in[-2,0]$,  the  time step is simply\n",
    "\n",
    "$$\n",
    "    \\Delta t = \\frac{z}{\\lambda}\n",
    "$$\n",
    "so the maximum time step we can take for $\\lambda<0$ is\n",
    "\n",
    "$$\n",
    "    \\Delta t = -\\frac{2}{\\lambda} > 0\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler for different values of time steps and $\\lambda = -1$\n",
    "lam = -1.\n",
    "f = lambda t, u :  lam * u\n",
    "u_exact = lambda t: numpy.exp(lam * t)\n",
    "u_0 = u_exact(0.)\n",
    "\n",
    "t_span = [0., 10.]\n",
    "\n",
    "N = 21\n",
    "N = 11\n",
    "N = 9\n",
    "N = 6\n",
    "\n",
    "N = 5\n",
    "abs_z = numpy.abs(lam*t_span[-1]/(N-1))\n",
    "print('Abs(z)={}'.format(abs_z))\n",
    "t, u = euler(f, t_span, u_0, N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "tt = numpy.linspace(0.,t_span[1],100)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(tt, u_exact(tt),label='$u_{{exact}}$')\n",
    "axes.plot(t, u, 'ro-', label='$u_{{euler}}$')\n",
    "axes.set_xlabel('t')\n",
    "axes.set_title('$\\Delta t = {}$, $|z|={}$'.format(t[1]-t[0],abs_z), fontsize=18)\n",
    "axes.grid()\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute Stability of the Forward Euler Method for complex $z$\n",
    "\n",
    "The region $|1 + z| < 1$  is the **region of absolute stability**.  For more general problems $\\lambda$ may be complex (stemming from eigenvalues of real problems),  so we generally consider the region of stability on the complex plane.\n",
    "\n",
    "$$\n",
    "    | 1 + z | < 1\n",
    "$$\n",
    "for $z\\in\\mathbb{C}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure(figsize=(8,6))\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), color=(255.0/255.0,145.0/255.0,0/255.0,1.0))\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_xlabel('Re(z)')\n",
    "axes.set_ylabel('Im(z)')\n",
    "\n",
    "\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Forward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Using the stability plots to choose a stable $\\Delta t$\n",
    "\n",
    "If you know the spectrum of eigenvalues $\\lambda$ for a given problem,  you can plot them on the stability diagram and then find the value of $\\Delta t$ such that all $z_i = \\lambda_i\\Delta t$ are within the stability region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "E.g. consider Euler's method on a single eigenvalue \n",
    "$$\n",
    "    \\lambda = -2 + 2i\n",
    "$$\n",
    "with modulus $|\\lambda| = 2\\sqrt{2}$ which  plots outside of the stability region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "lam = (-2,2)\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), color=(255.0/255.0,145.0/255.0,0/255.0,1.0))\n",
    "axes.plot(lam[0],lam[1],'rx',markersize=8)\n",
    "axes.plot((lam[0], 0.), (lam[1], 0.), 'r--')\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_xlabel('Re(z)')\n",
    "axes.set_ylabel('Im(z)')\n",
    "\n",
    "\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Forward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "However $z=\\lambda\\Delta t$ so \n",
    "$$|z|=\\Delta t|\\lambda|$$\n",
    "\n",
    "Thus absolute stability requires that we find a time-step such that $|1 + z| < 1$. If we just look for the intersection with the boundaries of the stability region\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    | 1 +z| &= | 1 +\\Delta t(-2+2i) | = 1\\\\\n",
    "    & = \\sqrt{ (1-2\\Delta t)^2 + 4\\Delta t^2} = 1 \\\\\n",
    "\\end{align}\n",
    "$$ \n",
    "\n",
    "implies (after a little work) that \n",
    "\n",
    "$$\n",
    "    \\Delta t(2\\Delta t - 1) = 0\n",
    "$$\n",
    "or $\\Delta t \\in (0, 1/2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "lam = (-2,2)\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), color=(255.0/255.0,145.0/255.0,0/255.0,1.0))\n",
    "axes.plot(lam[0],lam[1],'rx',markersize=8)\n",
    "axes.text(lam[0],1.1*lam[1], '$\\lambda$',fontsize=18)\n",
    "axes.plot(0.5*lam[0],0.5*lam[1],'kx',markersize=8)\n",
    "axes.plot(0.5*lam[0],0.5*lam[1],'kx',markersize=8)\n",
    "axes.text(0.5*lam[0],0.6*lam[1], '$\\Delta t = 1/2$', fontsize=18)\n",
    "axes.plot((lam[0], 0.), (lam[1], 0.), 'r--')\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_xlabel('Re(z)')\n",
    "axes.set_ylabel('Im(z)')\n",
    "\n",
    "\n",
    "\n",
    "axes.set_title(\"Absolute Stability Region for Forward Euler\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Why complex eigenvalues?  And why should we care?\n",
    "\n",
    "#### A simple example:  The simple harmonic oscillator\n",
    "\n",
    "The linear harmonic oscillator can be written as a 2nd-order ODE for position $x$ of a mass on a spring\n",
    "$$\n",
    "    \\frac{d^2 x}{dt^2} + x = 0, \\quad x(0)=x_0,\\,  \\frac{dx}{dt} = 0.\n",
    "$$\n",
    "Which has solutions \n",
    "$$\n",
    "    x(t) = x_0\\cos(t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Alternatively, we can write this as a 2 dimensional system of first order equations by  defining the velocity $v= dx/dt$ and\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\frac{dx}{dt} &= v \\\\\n",
    "    \\frac{dv}{dt} &= -x \\\\\n",
    "\\end{align}\n",
    "$$ \n",
    "or more cleanly as \n",
    "$$\n",
    "    \\frac{d\\mathbf{u}}{dt} = A\\mathbf{u},\\quad\\mathbf{u}(0)=\\mathbf{u_0}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\mathbf{u}=\\begin{bmatrix} x\\\\ v\\\\\\end{bmatrix}\\quad A= \\begin{bmatrix} 0 & 1 \\\\ -1 & 0 \\\\ \\end{bmatrix}\\quad\\mathbf{u}_0=\\begin{bmatrix} x_0\\\\ 0\\\\\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All diagonalizable linear dynamical systems have solutions of the form \n",
    "\n",
    "$$\n",
    "    \\mathbf{u}(t) = \\sum_{i=1}^n c_i e^{\\lambda_i t}\\mathbf{s}_i\n",
    "$$ \n",
    "\n",
    "where $\\lambda_i$, $\\mathbf{s}_i$ are corresponding eigenvalues and eigenvectors of $A$ and each eigenvector satisfies\n",
    "\n",
    "$$\n",
    "    \\frac{d\\mathbf{s}_i}{dt} = \\lambda_i\\mathbf{s}_i\n",
    "$$\n",
    "\n",
    "which is exactly our model problem.\n",
    "\n",
    "In the case where \n",
    "$$\n",
    "A= \\begin{bmatrix} 0 & 1 \\\\ -1 & 0 \\\\ \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The eigenvalues are _ _ _ _ _ _ _ _ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Even with complex Eigenvalues and Eigenvectors, the solutions are real and can be shown to be perfect circles on the phase plane \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    x(t) &= \\cos(t)\\\\\n",
    "    v(t) &=-\\sin(t)\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "However,  the plotting the eigenvalues on the stability diagram show's that Euler's method will be completely unstable for this problem for all time steps $\\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "x = numpy.linspace(-2, 2, 11)\n",
    "v = numpy.linspace(-2, 2, 11)\n",
    "X,V = numpy.meshgrid(x,v)\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.quiver(X, V, V, -X)\n",
    "axes.plot(1.5*numpy.cos(t), -1.5*numpy.sin(t))\n",
    "axes.plot([-2, 2],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-2, 2],'k--')\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('v',fontsize=16)\n",
    "\n",
    "axes.set_title('Phase Plane: simple Harmonic oscillator',fontsize=16)\n",
    "axes.set_aspect('equal')\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), color=(255.0/255.0,145.0/255.0,0/255.0,1.0))\n",
    "axes.plot([-3, 3],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-3, 3],'k--')\n",
    "axes.plot([0., 0.0], [-1., 1.],'rx',markersize=10)\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\",fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda t, u : numpy.array([ u[1], -u[0]])\n",
    "t_span = [0., 4*numpy.pi]\n",
    "u_0 = numpy.array([1., 0.])\n",
    "N = 100\n",
    "t, u = euler(f,t_span,u_0, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.quiver(X, V, V, -X)\n",
    "axes.plot(1.*numpy.cos(t), -numpy.sin(t))\n",
    "axes.plot(u[0],u[1],'r')\n",
    "axes.plot([-2, 2],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-2, 2],'k--')\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('v',fontsize=16)\n",
    "\n",
    "axes.set_title(\"Euler's method: simple Harmonic oscillator\", fontsize=16)\n",
    "axes.set_aspect('equal')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Absolute Stability of Single-Step Multi-Stage schemes -- The R method\n",
    "\n",
    "This approach can be applied to **all** Single-Step Multi-Stage schemes including all the Runge-Kutta and Taylor Series methods.  The recipe is straightforward\n",
    "\n",
    "*  Apply the stepping scheme for one step of the model problem $u' = \\lambda u$, $u(0)=u_0$ assuming complex $\\lambda$\n",
    "*  This will result in a discrete approximation\n",
    "\n",
    "$$\n",
    "    U(\\Delta t) = R(z) U_0\n",
    "$$\n",
    "\n",
    "where $z=\\lambda\\Delta t\\in\\mathbb{C}$, which will be a discrete approximation to the true solution\n",
    "\n",
    "$$\n",
    "    U(\\Delta t) = e^{z} U_0\n",
    "$$\n",
    "* Absolute stability will require that $|R(z)| < 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute Stability of the backward Euler Method\n",
    "\n",
    "Now try this on backward Euler.\n",
    "\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t f(t_{n+1}, U_{n+1}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which for our model problem is\n",
    "\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t\\lambda U_{n+1}\n",
    "$$\n",
    "or \n",
    "\n",
    "$$\n",
    "    (1 - \\Delta t\\lambda) U_{n+1} = U_n\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or rearranging and solving for $U_{n+1}$ gives\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "    U_{n+1} &= \\frac{1}{1-\\Delta t\\lambda} U_n\\\\\n",
    "            &= \\frac{1}{1-z}U_n\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "so\n",
    "$$\n",
    "    R(z) = \\frac{1}{1-z}\n",
    "$$ \n",
    "and absolute stability requires \n",
    "\n",
    "$$\n",
    "    |R(z)| \\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$ for $\\lambda<0$ (and even some positive values of $\\lambda$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Backward Euler\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "axes.set_facecolor((255.0/255.0,145.0/255.0,0/255.0,1.0))\n",
    "axes.fill(numpy.cos(t) + 1.0, numpy.sin(t), 'w')\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": [
    "However, stability doesn't imply accuracy either.  In this case, backwards Euler will still fail after long times for the simple harmonic oscillator as well as it will damp towards zero.  For the general Linear autonomous dynamical system\n",
    "\n",
    "$$\n",
    "    \\frac{d\\mathbf{u}}{dt} = A\\mathbf{u}, \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Backwards Euler method can be written as\n",
    "\n",
    "$$\n",
    "        \\mathbf{U}_{n+1} = \\mathbf{U}_n + \\Delta t A\\mathbf{U}_{n+1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or rearranging, requires solving the linear algebraic problem at each step\n",
    "\n",
    "$$\n",
    "    (I - \\Delta t A)\\mathbf{U}_{n+1} = \\mathbf{U}_n\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Backwards Euler for a linear system $u' = Au$\n",
    "def beuler(A, t_span, u0, N):\n",
    "    \"\"\" simple implementation of constant step-size forward euler method for vector valued f\n",
    "        This doc string should have so much more in it\n",
    "    \"\"\"\n",
    "    t = numpy.linspace(t_span[0], t_span[1],N)\n",
    "    delta_t = t[1] - t[0]\n",
    "\n",
    "    u = numpy.empty((len(t),len(u0)))\n",
    "    B = numpy.eye(len(u0)) - delta_t*A\n",
    "    u[0] = u0\n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        u[n + 1] = numpy.linalg.solve(B,u[n])\n",
    "    return t, u.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.array([[0, 1], [-1, 0]])\n",
    "t_span = [0., 4*numpy.pi]\n",
    "u_0 = numpy.array([1., 0.])\n",
    "N = 10000\n",
    "t, u = beuler(A,t_span,u_0, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.quiver(X, V, V, -X)\n",
    "axes.plot(1.*numpy.cos(t), -numpy.sin(t))\n",
    "axes.plot(u[0],u[1],'r')\n",
    "axes.plot([-2, 2],[0.0, 0.0],'k--')\n",
    "axes.plot([0.0, 0.0],[-2, 2],'k--')\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('v',fontsize=16)\n",
    "\n",
    "axes.set_title(\"Backward Euler: simple Harmonic oscillator\", fontsize=16)\n",
    "axes.set_aspect('equal')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## General Stability Regions for Linear Methods\n",
    "\n",
    "We can also map out  stability regions for Linear Multi-step methods (although the analysis is slightly more complicated).\n",
    "\n",
    "All of the Linear Multi-step methods from the previous lecture can be written\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j U_{n+j} = \\Delta t \\sum^r_{j=0} \\beta_j f(t,U_{n+j})\n",
    "$$\n",
    "\n",
    "which for our model problem, $f(t,u) = \\lambda u$ becomes\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j U_{n+j} = \\Delta t \\sum^r_{j=0} \\beta_j \\lambda U_{n+j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or  rearranging and again defining $z=\\Delta t\\lambda$ becomes\n",
    "\n",
    "$$\n",
    "    \\sum^r_{j=0} (\\alpha_j - \\beta_j z) U_{n+j} = 0.\n",
    "$$\n",
    "which takes the form of a \"Linear Difference Equation\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Characteristic Polynomials and Linear Difference Equations\n",
    "\n",
    "To help us analyze stability regions we need to take a small aside and consider linear difference equations and their solutions.  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}$.  This expression 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",
    "by dividing by $\\xi^n$.  If $\\xi$ then is a root of the polynomial\n",
    "\n",
    "$$\\rho(\\xi) \\equiv \\sum^r_{j=0} \\alpha_j \\xi^j$$\n",
    "\n",
    "then $\\xi$ solves the original difference equation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  the Fibonacci Numbers\n",
    "\n",
    "A classic example of  linear difference equation is the Fibonnacci sequence, $0,1,1,2,3,5,\\ldots$. \n",
    "\n",
    "Which are solutions to the two-term difference equation\n",
    "\n",
    "$$\n",
    "    U_{n+2} = U_{n+1} + U_n, \\quad U_0=0,\\, U_1 = 1\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "    U_{n+2} - U_{n+1} - U_n = 0\n",
    "$$\n",
    "\n",
    "or $\\alpha = \\begin{bmatrix} -1 & -1 & 1 \\\\\\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Substituting $U_n = \\xi^n$ into \n",
    "$$\n",
    "U_{n+2} - U_{n+1} - U_n =0\n",
    "$$\n",
    "yields\n",
    "$$\n",
    "    \\xi^n\\left(\\xi^2 - \\xi -1\\right) = 0\n",
    "$$\n",
    "which has solutions $\\xi =  \\frac{1\\pm\\sqrt{5}}{2}$.  Where the positive root is our favorite number $\\phi=1.618033$ which is the golden ratio.  The negative root is $1-\\phi$. Note, these roots are related to the eigenvalues of a related problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Thus all Fibonacci numbers can be written as a linear combination of the two roots\n",
    "\n",
    "$$\n",
    "    U_n = A\\phi^n + B(1-\\phi)^n\n",
    "$$\n",
    "\n",
    "where $A$ and $B$ are determined by the first two terms in the sequence [0,1].\n",
    "\n",
    "Note, in this case $\\phi > 1$ (while $|(1-\\phi)|<1$) thus $U_n\\rightarrow\\infty$ as $n\\rightarrow\\infty$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stability for LMM's\n",
    "Anyway,  Back to the Linear difference equations arising from applying any Linear multistep method to the model problem\n",
    "$$\n",
    "    \\sum^r_{j=0} (\\alpha_j - \\beta_j z) U_{n+j} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Letting\n",
    "$$\n",
    "    \\rho(\\xi) = \\sum^r_{j=0} \\alpha_j \\xi^j, \\quad \\sigma(\\xi) = \\sum^r_{j=0} \\beta_j \\xi^j\n",
    "$$\n",
    "and substituting $U_n = \\xi^n$ into the difference equations reduces to \n",
    "$$\n",
    "    \\pi(\\xi, z) = \\rho(\\xi) - z \\sigma(\\xi) = 0\n",
    "$$\n",
    "where $\\pi(\\xi, z)$ is the **stability polynomial** of the method, whose roots are solutions of the difference method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Absolute stability of the linear multi-step methods\n",
    "\n",
    "For $U_n$ to stay bounded as $n\\rightarrow\\infty$,  it requires that all roots $\\xi_i$ of $\\pi(\\xi, z)$ satisfy\n",
    "$$\n",
    "    |\\xi_i| \\leq 1\n",
    "$$\n",
    "If this is the case, then the multi-step method is **zero-stable**.  Roughly you can think of these roots as multiplying the past truncation errors so that if they are less than or equal to 1 they do not grow in magnitude.  We then define the region of **absolute stability** as the values for $z$ for which this is true.  This approach reduces to the $R$ method for  one-step methods as a special case."
   ]
  },
  {
   "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",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "setting $U_n = \\xi^n$ implies\n",
    "$$\n",
    "\\begin{aligned}\n",
    "   0 &= \\xi^n\\left(\\xi -  (1 + z)\\right)\\\\\n",
    "\\end{aligned}$$\n",
    "\n",
    "or $\\pi(\\xi, z) = \\xi -  (1 + z)$ 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": [
    "### Example:  Backwards Euler's Method\n",
    "\n",
    "Similarly for backwards Euler\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    0 &= U_{n+1} - U_n - \\Delta t \\lambda U_{n+1} \\\\\n",
    "    &= U_{n+1}(1 - \\Delta t\\lambda) - U_n \\\\\n",
    "    &= \\xi(1 - z)  - 1 \\\\\n",
    "    &=\\pi(\\xi, z)\n",
    "\\end{aligned}$$\n",
    "\n",
    "whose root is $\\xi = 1/(1 - z)$  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example: Adams-Bashforth 2\n",
    "\n",
    "As an example consider the Adams-Bashforth 2-stage LMM\n",
    "$$\n",
    "    U_{n+2} = U_{n+1} + \\frac{\\Delta t}{2} (-f(U_n) + 3 f(U_{n+1}))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let $f(U) = \\lambda U$,  $z = \\Delta t \\lambda$\n",
    "\n",
    "$$\n",
    "    U_{n+2} - U_{n+1} - \\frac{1}{2} (-z U_n + 3 z U_{n+1}) = 0 \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "substituting $U_n=\\xi^n$, and factoring out $\\xi^n$ yields\n",
    "$$\n",
    "    \\pi(\\xi, z) = \\xi^{2} - \\xi - \\frac{1}{2} (-z  + 3 z \\xi)\n",
    "$$\n",
    "with roots when\n",
    "\n",
    "$$\n",
    "     2 \\xi^2 - (2 + 3z) \\xi + z =0\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plotting Stability Regions\n",
    "\n",
    "Given a stability polynomial $\\pi(\\xi,z)$, absolute stability requires that \n",
    "\n",
    "$$\n",
    "    \\max{|\\xi_i|} <= 1\n",
    "$$\n",
    "where $\\xi_i$ are the (generally complex) roots of $\\pi(\\xi, z)$ and therefore a function of $z$.  In the case of single step schemes, the problem reduces to finding the region of the complex plane where $|R(z)| < 1$.  \n",
    "\n",
    "Either way,  the simplest way to visualize the regions of the complex plane where the scheme is stable, is to calculate either $|R(z)|$ or $\\max|\\xi_i(z)|$ and plot the 1-contour and determine all regions that are $<1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Some useful codes\n",
    "\n",
    "study the following cells for some useful routines for plotting stability regions for both Single-step multi-stage schemes and Linear Multi-step schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def stability_plot(X, Y, C, axes, title=None, continuous=True):\n",
    "    \"\"\"\n",
    "    Utility function to make stability diagram given complex stability scalar C\n",
    "    \n",
    "    parameters:\n",
    "    -----------\n",
    "    \n",
    "    X, Y:  numpy.meshgrids for complex plane\n",
    "    C: numpy array\n",
    "        Field to plot,  either |R(z)| for a single step scheme, or max(|xi_i(z)|) for a LMM scheme\n",
    "    axes: matplotlib axes object\n",
    "        subplot or plot to draw in. \n",
    "    title: string\n",
    "        subplot title if not None \n",
    "    continuous: boolean\n",
    "        if True, plot a continous coloring of C\n",
    "        if False, plot Heaviside(C)\n",
    "    \"\"\"\n",
    "    if  continuous:\n",
    "        Ch = C\n",
    "    else:\n",
    "        Ch = numpy.heaviside(C-1,0.)\n",
    "    pcolor_plot = axes.pcolor(X, Y, Ch, vmin=0, vmax=1, cmap=plt.get_cmap('Greens_r'), shading='auto')\n",
    "    axes.contour(X, Y, C, 'k', levels=[1.0])\n",
    "    fig = plt.gcf()\n",
    "    fig.colorbar(pcolor_plot)\n",
    "    axes.plot(x, numpy.zeros(x.shape),'k--')\n",
    "    axes.plot(numpy.zeros(y.shape), y,'k--')\n",
    "    \n",
    "    axes.set_xlabel('Re', fontsize=16)\n",
    "    axes.set_ylabel('Im', fontsize=16)\n",
    "    if title is not None:\n",
    "        axes.set_title(title, fontsize=16)\n",
    "    \n",
    "    axes.set_aspect('equal')    \n",
    "    \n",
    "def plot_stability_ssms(R, x, y, axes=None, title=None, continuous=True):\n",
    "    \"\"\" \n",
    "    plot stability regions for single-step multi-stage ODE schemes given the function R(z)\n",
    "    such that U_{n+1} = R(z)U_n   and z is complex\n",
    "    \n",
    "    parameters:\n",
    "    -----------\n",
    "    \n",
    "    R: calleable\n",
    "        function of a complex variable z such that if |R|<=1, the scheme is absolutely stable\n",
    "    x: numpy array\n",
    "        array of values for the real axis\n",
    "    y: numpy array\n",
    "        values to plot for the imaginary axis\n",
    "    axes: matplotlib axes object\n",
    "        subplot or plot to draw in.  If axes=None create a new figure\n",
    "    title: string\n",
    "        subplot title if \n",
    "    continuous: boolean\n",
    "        if True, plot a continous coloring of C\n",
    "        if False, plot Heaviside(C)\n",
    "    \"\"\"\n",
    "    \n",
    "    X,Y = numpy.meshgrid(x,y)\n",
    "    Z = X + 1j * Y\n",
    "    if axes is None:\n",
    "        fig = plt.figure(figsize=(8,6))\n",
    "        axes = fig.add_subplot(1,1,1)\n",
    "    \n",
    "    abs_R = numpy.abs(R(Z))\n",
    "    stability_plot(X, Y, abs_R, axes, title, continuous)\n",
    "    \n",
    "def plot_stability_lmm(pi_coeff, x, y, axes=None, title=None, continuous=True):\n",
    "    \"\"\" \n",
    "    plot stability regions for linear multi-step  ODE schemes given the coefficients of the stability polynomial\n",
    "    pi(xi, z)\n",
    "    \n",
    "    parameters:\n",
    "    -----------\n",
    "    \n",
    "    pi_coeff: calleable (function of z)\n",
    "        function that returns array of stability polynomial pi(z)\n",
    "    x: numpy array\n",
    "        array of values for the real axis\n",
    "    y: numpy array\n",
    "        values to plot for the imaginary axis\n",
    "    axes: matplotlib axes object\n",
    "        subplot or plot to draw in.  If axes=None create a new figure\n",
    "    title: string\n",
    "        subplot title if not None   \n",
    "    continuous: boolean\n",
    "        if True, plot a continous coloring of C\n",
    "        if False, plot Heaviside(C)\n",
    "    \"\"\"\n",
    "       \n",
    "    X,Y = numpy.meshgrid(x,y)\n",
    "    Z = X + 1j * Y\n",
    "    if axes is None:\n",
    "        fig = plt.figure(figsize=(8,6))\n",
    "        axes = fig.add_subplot(1,1,1)\n",
    "    \n",
    "    norm_max = numpy.empty(Z.shape)\n",
    "    for i,row in enumerate(Z):\n",
    "        for j, z in enumerate(row):\n",
    "            norm_max[i,j] = max(numpy.abs(numpy.roots(pi_coeff(z))))\n",
    "    \n",
    "    stability_plot(X, Y, norm_max, axes, title, continuous)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-3,3,100)\n",
    "y = numpy.linspace(-3,3,100)\n",
    "\n",
    "R_euler = lambda z: 1 + z\n",
    "R_beuler = lambda z: 1./(1. - z)\n",
    "R_trap = lambda z: (1. + 0.5*z)/(1. - 0.5*z)\n",
    "R_RK2 = lambda z: 1 + z + z**2/2\n",
    "R_Taylor4 = lambda z: 1 + z + z**2/2 + z**3/6. + z**4/24.\n",
    "\n",
    "continuous=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,12))\n",
    "axes = fig.add_subplot(2,2,1)\n",
    "plot_stability_ssms(R_euler, x, y, axes=axes, title='Stability of Eulers method', continuous=continuous)\n",
    "\n",
    "axes = fig.add_subplot(2,2,2)\n",
    "plot_stability_ssms(R_beuler, x, y, axes=axes, title='Stability of Backwards Eulers', continuous=continuous)\n",
    "\n",
    "axes = fig.add_subplot(2,2,3)\n",
    "plot_stability_ssms(R_trap, x, y, axes=axes, title='Stability of Trapezoidal method', continuous=continuous)\n",
    "\n",
    "axes = fig.add_subplot(2,2,4)\n",
    "plot_stability_ssms(R_Taylor4, x, y, axes=axes, title='Stability of Taylor4 method', continuous=continuous)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-3,3,100)\n",
    "y = numpy.linspace(-3,3,100)\n",
    "\n",
    "pi_euler = lambda z: numpy.array([1., -(1 + z)])\n",
    "pi_AB2 = lambda z: numpy.array([2, -(2. + 3*z), z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "plot_stability_lmm(pi_euler, x, y, axes=axes, title='Stability of Eulers method')\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "plot_stability_lmm(pi_AB2, x, y, axes=axes, title='Stability of Adams-Bashforth2 method')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Application to Stiff ODEs\n",
    "\n",
    "Consider  the ODE IVP \n",
    "\n",
    "$$u'(t) = \\lambda (u - \\cos t) - \\sin t,\\quad u(t_0)=u_0$$\n",
    "\n",
    "The general solution of the ODE is\n",
    "\n",
    "$$u(t) = e^{\\lambda (t - t_0)} (u_0 - \\cos t_0) + \\cos t$$\n",
    "\n",
    "For $u_0=1$ at $t_0=0$,  the solution is simply $\\cos t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What happens to solutions that are slightly different from $u_0 = 1$ or $t_0 = 0$? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$u(t) = e^{\\lambda (t - t_0)} (u_0 - \\cos t_0) + \\cos t$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Plot \"hairy\" solutions to the ODE\n",
    "u = lambda t_0, u0, lam, t: numpy.exp(lam * (t - t_0)) * (u0 - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 10))\n",
    "\n",
    "for i,lam in enumerate([-1, -10, 1, 10]):\n",
    "    axes = fig.add_subplot(2, 2, i+1)\n",
    "    for u0 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, u0, lam, t),'b')\n",
    "    t = numpy.linspace(0.0,10.0,100)\n",
    "    axes.plot(t, numpy.cos(t), 'r', linewidth=5)\n",
    "    axes.set_ylim((-2, 2))\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((-2, 2))\n",
    "plt.show()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot \"inverse hairy\" solutions to the ODE\n",
    "u = lambda t_0, u0, lam, t: numpy.exp(lam * (t - t_0)) * (u0 - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\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()\n",
    "        "
   ]
  },
  {
   "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",
    "We can write this problem as a system of Linear ODE's\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{dA}{dt} &= -K_1 A\\\\\n",
    "    \\frac{dB}{dt} &= K_1 A -K_2B\\\\\n",
    "    \\frac{dC}{dt} &= K_2B\\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "with initial condition $A_0, B_0, C_0$\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If we let\n",
    "\n",
    "$$\\mathbf{u} = \\begin{bmatrix} A \\\\ B \\\\ C \\end{bmatrix}$$\n",
    "\n",
    "then we can rewrite this as a linear dynamical system\n",
    "$$\n",
    " \\frac{d\\mathbf{u}}{dt} = M\\mathbf{u},  \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$ \n",
    "where \n",
    "\n",
    "$$ M = \n",
    "\\begin{bmatrix}\n",
    "    -K_1 & 0 & 0 \\\\\n",
    "    K_1 & -K_2 & 0 \\\\\n",
    "    0 & K_2 & 0\n",
    "\\end{bmatrix} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $A$ is diagonalizable, the general solution of a linear dynamical system  can be written in terms of the eigenvalues ($\\lambda_i$) and eigenvectors $\\mathbf{x}_i$ of $M$ as\n",
    "\n",
    "$$\\mathbf{u}(t) = c_{1} e^{\\lambda_1 t}\\mathbf{x}_1 + c_{2}e^{\\lambda_2t}\\mathbf{x}_2 + c_{3}e^{\\lambda_3t}\\mathbf{x}_3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And given our matrix \n",
    "$$ M = \n",
    "\\begin{bmatrix}\n",
    "    -K_1 & 0 & 0 \\\\\n",
    "    K_1 & -K_2 & 0 \\\\\n",
    "    0 & K_2 & 0\n",
    "\\end{bmatrix} \n",
    "$$\n",
    "\n",
    "The eigenvalues are ??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And the general solution is \n",
    "\n",
    "$$\\mathbf{u}(t) = c_{1} e^{-K_1t}\\mathbf{x}_1 + c_{2}e^{-K_2t}\\mathbf{x}_2 + c_{3}\\mathbf{x}_3$$\n",
    "\n",
    "and the constants $c_i$ can be found by solving $X\\mathbf{c} = \\mathbf{u}_0$ which describe the initial condition as a linear combination of the eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Solve the chemical systems example using a Forward Euler scheme\n",
    "# Problem parameters\n",
    "K_1 = 1\n",
    "K_2 = 90.\n",
    "N =  1000\n",
    "\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",
    "u_0 = [3, 5.0, 0.0]\n",
    "t_span = [0, 3.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "f = lambda t, u : A.dot(u)\n",
    "t, U = euler(f, t_span, u_0, N)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t, U.transpose())\n",
    "axes.legend([\"[A]\", \"[B]\", \"[C]\"])\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",
    "#axes.set_xlim((0.0, 2.0))\n",
    "axes.grid()\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 linear dynamical systems like the chemical decay problem,  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_2 / K_1$ if $K_2 > K_1$.  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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# set the eigenvalues and time step\n",
    "K_1 = 1.0\n",
    "K_2 = 3.0\n",
    "delta_t = 1.0\n",
    "eigenvalues = [-K_1, -K_2]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the region of absolute stability for Forward Euler\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "t = numpy.linspace(0.0, 2.0 * numpy.pi, 100)\n",
    "\n",
    "\n",
    "axes.fill(numpy.cos(t) - 1.0, numpy.sin(t), color=(255.0/255.0,145.0/255.0,0/255.0,1.0))\n",
    "for lam in eigenvalues:\n",
    "    print('z = {}'.format(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": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# solve using backwards Euler\n",
    "K_1 = 1.0\n",
    "K_2 = 100.0\n",
    "A = numpy.array([[-K_1, 0, 0], [K_1, -K_2, 0], [0, K_2, 0]])\n",
    "\n",
    "N = 10000 \n",
    "t_span = [0., 3.]\n",
    "U0 = [3., 5.0, 0.]\n",
    "\n",
    "t, U = beuler(A, t_span, U0, N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "    \n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t, U.T)\n",
    "axes.legend([\"[A]\", \"[B]\", \"[C]\"])\n",
    "axes.set_title(\"Chemical System\")\n",
    "axes.set_xlabel(\"$t$\")\n",
    "axes.set_title(\"$[A], [B], [C]$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "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 as a stability region of\n",
    "\n",
    "$$|1 - z| \\geq 1$$\n",
    "\n",
    "which covers the full left-half plane.  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.\n",
    "\n",
    "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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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((255.0/255.0,145.0/255.0,0/255.0,1.0))\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": "slide"
    }
   },
   "source": [
    "# L-Stability\n",
    "\n",
    "It turns out not all A-stable methods are alike.  Consider the backward Euler method and the trapezoidal method defined by\n",
    "\n",
    "$$\\frac{U_{n+1} - U_n}{\\Delta t} = \\frac{1}{2}(f(U_n) + f(U_{n+1}))$$\n",
    "\n",
    "whose stability polynomial (or $R(z))$ is\n",
    "\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",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or \n",
    "$$\n",
    "      \\pi(\\xi,z) = \\left(1 - \\frac{z}{2}\\right ) \\xi - \\left(1 + \\frac{z}{2}\\right)\n",
    "$$ \n",
    "with single root\n",
    "$$\n",
    "\\xi(z) = R(z) = \\frac{1 + \\frac{z}{2}}{1 - \\frac{z}{2}}\n",
    "$$\n",
    "\n",
    "which shows that it is A-stable (i.e. $|R(z)| < 1$ for all $\\Re(z)<0$)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Let's apply both these methods the stiff problem\n",
    "\n",
    "$$u'(t) = \\lambda (u - \\cos t) - \\sin t,\\quad u(t_0)=u_0$$\n",
    "\n",
    "With solution \n",
    "\n",
    "$$u(t) = e^{\\lambda (t - t_0)} (u_0 - \\cos t_0) + \\cos t$$\n",
    "\n",
    "and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot \"hairy\" solutions to the ODE\n",
    "u = lambda t_0, u0, lam, t: numpy.exp(lam * (t - t_0)) * (u0 - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "\n",
    "for i,lam in enumerate([-2, -20]):\n",
    "    axes = fig.add_subplot(1, 2, i+1)\n",
    "    for u0 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, u0, lam, t),'b')\n",
    "    t = numpy.linspace(0.0,10.0,100)\n",
    "    axes.plot(t, numpy.cos(t), 'r', linewidth=5)\n",
    "    axes.set_ylim((-1.5, 1.5))\n",
    "    axes.grid()\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((-1.5, 1.5))\n",
    "plt.show()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler and Trapezoidal for the stiff ODE\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)) * (u0 - numpy.cos(t_0)) + numpy.cos(t)\n",
    "\n",
    "t_0 = 0.0\n",
    "t_f = 1.0\n",
    "eta = 2.0\n",
    "lam = -1e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "num_steps = [10, 20, 40, 50]\n",
    "# num_steps = numpy.arange(100, 1000, 100)\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "    \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='Trapezoidal')\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-stable**, where we require that the method is A-stable and that\n",
    "\n",
    "$$\\lim_{z \\rightarrow \\infty} |R(z)| = 0.$$\n",
    "\n",
    "Turns out that backwards Euler is L-stable while trapezoidal rule is not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "R_beuler = lambda z: 1./(1. - z)\n",
    "R_trap = lambda z: (1. + 0.5*z)/(1. - 0.5*z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,6))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "plot_stability_ssms(R_beuler, x, y, axes=axes, title='Stability of Backwards Eulers')\n",
    "\n",
    "axes = fig.add_subplot(1,2,2)\n",
    "plot_stability_ssms(R_trap, x, y, axes=axes, title='Stability of Trapezoidal')"
   ]
  },
  {
   "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 t\\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",
    "    \\text{BDF-1}~& (r = 1):& & U_{n+1} - U_n = \\Delta t f(U_{n+1}) \\\\\n",
    "    \\text{BDF-2}~& (r = 2):& &3 U_{n+2} - 4 U_{n+1} + U_n = 2 \\Delta t f(U_{n+2}) \\\\\n",
    "    \\text{BDF-3}~& (r = 3):& &11U_{n+3} - 18U_{n+2} + 9U_{n+1} - 2 U_n = 6 \\Delta t f(U_{n+3}) \\\\\n",
    "    \\text{BDF-4}~& (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": "fragment"
    }
   },
   "source": [
    "* Inspection of BDF-1 shows that it's equivalent to ??\n",
    "* You should also recognize BDF-2 from our Numerical differentiation notes$\\ldots$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from fdcoeffV import fdcoeffV\n",
    "# order of the BDF scheme (i.e. N=2 implies BDF-2)\n",
    "N=3\n",
    "# order of the derivative\n",
    "k = 1\n",
    "x = numpy.array(range(N+1))\n",
    "w = fdcoeffV(k, x[N], x)\n",
    "print('weights = {}'.format(w))\n",
    "print('weights x 6 = {}'.format(6*w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "## Finish the code off here for calculating BDF coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pi_BDF1 = lambda z: numpy.array([ 1. - z, -1.])\n",
    "pi_BDF2 = lambda z: numpy.array([ 3. - 2*z, -4. , 1.])\n",
    "pi_BDF3 = lambda z: numpy.array([ 11 - 6.*z, -18., 9., -2.])\n",
    "pi_BDF4 = lambda z: numpy.array([ 25. - 12*z, -48., 36, -16., 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-10,10,50)\n",
    "y = numpy.linspace(-10,10,50)\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "axes = fig.add_subplot(2,2,1)\n",
    "plot_stability_lmm(pi_BDF1, x, y, axes=axes, title='Stability of BDF1 method')\n",
    "\n",
    "axes = fig.add_subplot(2,2,2)\n",
    "plot_stability_lmm(pi_BDF2, x, y, axes=axes, title='Stability of BDF2 method')\n",
    "\n",
    "axes = fig.add_subplot(2,2,3)\n",
    "plot_stability_lmm(pi_BDF3, x, y, axes=axes, title='Stability of BDF3 method')\n",
    "\n",
    "axes = fig.add_subplot(2,2,4)\n",
    "plot_stability_lmm(pi_BDF4, x, y, axes=axes, title='Stability of BDF4 method')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### BDF Schemes\n",
    "\n",
    "#### Pro's\n",
    "\n",
    "* LMM's so only require 1 (non-linear) solve per time step\n",
    "* relatively easy to derive and code methods of variable order\n",
    "* L-Stable and A-Stable\n",
    "\n",
    "#### Con's\n",
    "* Still need to deal solving for system's of non-linear equations\n",
    "* LMM's are harder to make adaptive (but can do so using interpolation techniques)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example,  BDF-1 for a general non-linear ODE\n",
    "\n",
    "Consider the general $n$-dimensional nonlinear system of ODE's (e.g van-der pol oscillator or SIR problem) which can be written as \n",
    "\n",
    "$$\n",
    "    \\frac{d\\mathbf{u}}{dt} = \\mathbf{f}(t,\\mathbf{u}), \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$\n",
    "\n",
    "where\n",
    "$$\n",
    "    \\mathbf{u}(t) = \\begin{bmatrix} u_1(t) \\\\ u_2(t) \\\\ \\vdots \\\\ u_n(t) \\end{bmatrix},\\quad\n",
    "    \\mathbf{f}(t,\\mathbf{u}) = \\begin{bmatrix} f_1(t,\\mathbf{u}) \\\\ f_2(t,\\mathbf{u}) \\\\ \\vdots \\\\ f_n(t,\\mathbf{u}) \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and $\\mathbf{f}$ is a vector of scalar, possibly non-linear  functions of the solution vector $\\mathbf{u}$\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "for example, in the Van-der-pol oscillator $n=2$ and\n",
    "\n",
    "$$\n",
    "    \\mathbf{f}(t,\\mathbf{u}) = \\begin{bmatrix} u_2 \\\\ \\mu(1 - u_1^2)u_2 - u_1 \\end{bmatrix}\n",
    "$$\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## BDF-1 discretization (Backwards-Euler)\n",
    "\n",
    "So the BDF-1 (Backwards Euler) discretization of this system can be written\n",
    "\n",
    "$$\n",
    "    \\mathbf{U}_{n+1} - \\mathbf{U}_{n} - \\Delta t\\mathbf{f}(t,\\mathbf{U}_{n+1}) = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "Which itself is a more general system of non-linear Equations for $\\mathbf{U}_{n+1}$ of form\n",
    "\n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{U}_{n+1}) = \\begin{bmatrix} F_1(t,\\mathbf{U}_{n+1},\\mathbf{U}_n) \\\\\n",
    "    F_2(t,\\mathbf{U}_{n+1},\\mathbf{U}_n) \\\\\n",
    "    \\vdots\\\\\n",
    "    F_n(t,\\mathbf{U}_{n+1},\\mathbf{U}_n)\\\\ \\end{bmatrix} =\n",
    "    \\mathbf{0}\n",
    "$$\n",
    "\n",
    "where, in this case\n",
    "$$\n",
    "    F_i = {U_{i}}_{(n+1)} - {U_{i}}_{(n)} - \\Delta tf_i(t,\\mathbf{U}_{n+1})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Solution of systems of Non-linear equations:  Newton's method\n",
    "\n",
    "If our system had only one scalar variable and function $f(t,u)$, then each step would reduce to solving the scalar non-linear problem\n",
    "$$\n",
    "    F(u_{n+1}) = 0\n",
    "$$\n",
    "which is a rootfinding problem, and we could use any of our rootfinders to solve this (in particular we could use robust bracketing schemes like brent).  Unfortunately once we move to higher-dimensions the problems become significantly harder (and solutions may not exist).  However we can still attempt to apply Newton's method to the general non-linear problem\n",
    "\n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{x}) = \\mathbf{0}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##  Newton's method: preliminaries\n",
    "\n",
    "As a simpler non-linear system, we will consider the system of equations\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    x_1^2 + 2x_2^2 &=\\gamma_1 \\\\\n",
    "    x_1 + 5x_2 &=\\gamma_2 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "which geometrically is the intersection of an ellipse and a straight line, which, depending on the values of $\\gamma_1$ and $\\gamma_2$ should have 0,1, or 2 solutions\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-4,4,100)\n",
    "y = numpy.linspace(-4,4,100)\n",
    "X,Y = numpy.meshgrid(x,y)\n",
    "\n",
    "Z1 = lambda gamma: X**2 + 2*Y**2 - gamma\n",
    "Z2 = lambda gamma: X + 5*Y - gamma\n",
    "\n",
    "gamma_1 = 10*numpy.ones(3)\n",
    "gamma_2 = [-2, 10, -15]\n",
    "fig = plt.figure(figsize=(24,6))\n",
    "\n",
    "for i in range(3):\n",
    "    axes = fig.add_subplot(1,3,i+1)\n",
    "    axes.contour(X,Y,Z1(gamma_1[i]),[0],colors='r')\n",
    "    axes.contour(X,Y,Z2(gamma_2[i]),[0],colors='b')\n",
    "    axes.axis('equal')\n",
    "\n",
    "    axes.set_xlabel('$x_1$')\n",
    "    axes.set_ylabel('$x_2$')\n",
    "    axes.set_title('$\\gamma_1={}, \\gamma_2={}$'.format(gamma_1[i],gamma_2[i]))\n",
    "    axes.grid()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now this system is sufficiently simple that we could solve it using substitution and solve\n",
    "a quadratic equation for the roots (or even use a 1-D rootfinder).  But in general we would like a more general method where we can start from some initial guess and rapidly converge to one of the roots.  Newton's method provides a general way to solve these sorts of problems (with all the issues that go with Newton for scalar problems)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##  Newton's method: Review $f(x)=0$\n",
    "\n",
    "Newton's method for a scalar function of a single variable says, given any $x$ such that $f(x)\\neq 0$, there should be a correction $\\delta$ such that\n",
    "\n",
    "$$\n",
    "    f(x+\\delta) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Or expanding in a Taylor's series around $x$\n",
    "\n",
    "$$\n",
    "    f(x+\\delta) = f(x) + f'(x)\\delta + O(f''\\delta^2) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Neglecting the higher order terms this becomes  \n",
    "\n",
    "$$\n",
    "    f(x) + f'(x)\\delta = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or solving for the correction\n",
    "$$\n",
    "    \\delta =\\frac{-f(x)}{f'(x)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If the correction were exact it would follow that $f(x+\\delta)=0$ and we would have our solution. \n",
    "\n",
    "In general,  the higher-order terms are important so the residual is not zero, but this gives us an iterative scheme that we hope will converge\n",
    "\n",
    "Set an initial guess $x_0$. Then iterate for $k=0,1,2\\ldots$ \n",
    "\n",
    "$$  \n",
    "\\begin{align}\n",
    "    \\text{Solve:} & & \\delta_k &= \\frac{-f(x_k)}{f'(x_k)}\\\\\n",
    "    & & x_{k+1} &= x_k + \\delta_k\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "until $|f(x_k)| < \\text{tol}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "If $f$ is well behaved with simple roots and the initial guess is within the basin of attraction of a given root,  Newton can find the root with quadratic convergence.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "It can also behave very badly and fly off to infinity..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "##  Newton's method: $\\mathbf{F}(\\mathbf{x})=\\mathbf{0}$\n",
    "\n",
    "Newton's method for a vector valued function is actually very similar (without the comfort of things such as brackets)\n",
    "\n",
    "Again for general $\\mathbf{x}$, $\\mathbf{F}(\\mathbf{x})\\neq\\mathbf{0}$ but is a vector often called the \"residual\".\n",
    "\n",
    "By the rules of vector norms \n",
    "\n",
    "$$\n",
    "    ||\\mathbf{F}(\\mathbf{x})|| \\geq 0\n",
    "$$\n",
    "\n",
    "and is only zero at a root of $\\mathbf{F}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If the residual is non-zero for some arbitrary $\\mathbf{x}$, Newton's method assumes that there is some correction *vector* $\\boldsymbol{\\delta}$ such that\n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{x}+\\boldsymbol{\\delta})=\\mathbf{0}\n",
    "$$\n",
    "    \n",
    "And we simply need to expand a vector-valued function in its Taylor series and truncate at linear terms to solve for the correction.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Preliminaries:  Taylor Series expansion of   $f(\\mathbf{x}):\\mathbb{R}^n\\rightarrow\\mathbb{R}$.  \n",
    "\n",
    "if $f(\\mathbf{x})$ is a scalar function of a vector valued variable (e.g. a surface like $f_1(\\mathbf{x}$), then the Taylor series of\n",
    "\n",
    "$$\n",
    "    f(\\mathbf{x} +\\boldsymbol{\\delta}) = f(\\mathbf{x}) + \\sum_{i=1}^n \\frac{\\partial f}{\\partial x_i}(\\mathbf{x})\\delta_i + \\frac{1}{2}\\sum_{j=1}^n\\sum_{i=1}^n \\frac{\\partial^2 f}{\\partial x_i\\partial x_j}(\\mathbf{x})\\delta_i\\delta_j + HOT\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or more beautifully in vector notation\n",
    "$$\n",
    "    f(\\mathbf{x} +\\boldsymbol{\\delta}) = f(\\mathbf{x}) + (\\nabla f)^T\\boldsymbol{\\delta} + \\frac{1}{2}\\boldsymbol{\\delta}^TH(\\mathbf{x})\\boldsymbol{\\delta} +  HOT\n",
    "$$\n",
    "\n",
    "where $\\nabla f$ is the gradient of $f$,  and $H$ is the Hessian, a symmetric matrix whose components are\n",
    "$$\n",
    "    H_{ij} = \\frac{\\partial^2 f}{\\partial x_i\\partial x_j}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "####  Taylor Series expansion of   $\\mathbf{F}(\\mathbf{x}):\\mathbb{R}^n\\rightarrow\\mathbb{R}^n$.  \n",
    "\n",
    "This is readily extended to vector valued functions as\n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{x}+\\boldsymbol{\\delta}) = \\begin{bmatrix}\n",
    "        F_1(\\mathbf{x}+\\boldsymbol{\\delta}) \\\\\n",
    "        F_2(\\mathbf{x}+\\boldsymbol{\\delta}) \\\\\n",
    "        \\vdots \\\\\n",
    "        F_n(\\mathbf{x}+\\boldsymbol{\\delta}) \\\\\n",
    "        \\end{bmatrix} \\approx \\begin{bmatrix}\n",
    "        F_1(\\mathbf{x}) + (\\nabla F_1)^T\\boldsymbol{\\delta} \\\\\n",
    "        F_2(\\mathbf{x}) + (\\nabla F_2)^T\\boldsymbol{\\delta} \\\\\n",
    "        \\vdots \\\\\n",
    "        F_n(\\mathbf{x}) + (\\nabla F_n)^T\\boldsymbol{\\delta} \\\\\n",
    "\\end{bmatrix} + HOT\n",
    "$$\n",
    "                        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or more succinctly as \n",
    "\n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{x}+\\boldsymbol{\\delta}) \\approx \\mathbf{F}(\\mathbf{x}) + J(\\mathbf{x})\\boldsymbol{\\delta}\n",
    "$$\n",
    "\n",
    "where $J$ is the ``Jacobian'' of $\\mathbf{F}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "####  The Jacobian\n",
    "\n",
    "The Jacobian of $\\mathbf{F}(\\mathbf{x})$ is the matrix whose rows are the gradients of the individual functions $F_i(\\mathbf{x})$\n",
    "\n",
    "$$\n",
    " J = \n",
    "    \\begin{bmatrix}\n",
    "        \\nabla F_1^T \\\\\n",
    "        \\nabla F_2^T \\\\\n",
    "        \\vdots \\\\\n",
    "        \\nabla F_n^T\\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or whose components are \n",
    "\n",
    "$$\n",
    "    J_{ij} = \\frac{\\partial F_i}{\\partial x_j}\n",
    "$$\n",
    "\n",
    "i.e. all the partial derivatives of all the functions with respect to all the variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Our toy non-linear problem \n",
    "$$\n",
    "\\begin{align}\n",
    "    x_1^2 + 2x_2^2 &=\\gamma_1 \\\\\n",
    "    x_1 + 5x_2 &=\\gamma_2 \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Can be rewritten in residual form as \n",
    "$$\n",
    "\\mathbf{F}(\\mathbf{x}) = \\begin{bmatrix}\n",
    "    x_1^2 + 2x_2^2 - \\gamma_1 \\\\\n",
    "    x_1 + 5x_2 - \\gamma_2 \\\\\n",
    " \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "with Jacobian?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "    J(\\mathbf{x}) = \\begin{bmatrix}\n",
    "    2x_1 &  4x_2  \\\\\n",
    "    1 & 5  \\\\\n",
    " \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Newton's method for $\\mathbf{F}(\\mathbf{x}) = \\mathbf{0}$ (finally)\n",
    "\n",
    "Given all the pieces,  raw Newton for $\\mathbf{F}(\\mathbf{x}):\\mathbb{R}^n\\rightarrow\\mathbb{R}^n$ is simply\n",
    "\n",
    "given an initial guess $\\mathbf{x}$ such that $\\mathbf{F}(\\mathbf{x})\\neq\\mathbf{0}$: seek a correction such that\n",
    "\n",
    "\n",
    "$$\n",
    "   \\mathbf{F}(\\mathbf{x}+\\boldsymbol{\\delta})\\approx \\mathbf{F}(\\mathbf{x}) + J\\boldsymbol{\\delta} = \\mathbf{0} \n",
    "$$ \n",
    "\n",
    "or solve\n",
    "\n",
    "$$\n",
    "    J\\boldsymbol{\\delta} = -\\mathbf{F}(\\mathbf{x})\n",
    "$$ \n",
    "\n",
    "and correct $\\mathbf{x}\\leftarrow\\mathbf{x}+\\boldsymbol{\\delta}$\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "or as an iterative scheme:\n",
    "\n",
    "Set an initial guess $x_0$. Then iterate for $k=0,1,2\\ldots$ \n",
    "\n",
    "$$  \n",
    "\\begin{align}\n",
    "    \\text{Solve:} & &  J(\\mathbf{x}_k)\\boldsymbol{\\delta_k} &= -\\mathbf{F}(\\mathbf{x}_k)\\\\\n",
    "    & &\\mathbf{x}_{k+1} &= \\mathbf{x}_k + \\boldsymbol{\\delta}_k\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "until $||\\mathbf{F}(\\mathbf{x}_k)|| < \\text{tol}$\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Minor digression: Linear systems\n",
    "\n",
    "We can actually use Newton to solve linear systems as well.  Consider the general linear problem\n",
    "\n",
    "$$\n",
    "    A\\mathbf{x} = \\mathbf{b}\n",
    "$$\n",
    "\n",
    "which we can rewrite in residual form as \n",
    "\n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{x}) =  A\\mathbf{x} - \\mathbf{b} = \\mathbf{0}\n",
    "$$\n",
    "\n",
    "i.e. the solution of our original problem is also a root of $\\mathbf{F}$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For any other initial guess $\\mathbf{x}_0$,  $\\mathbf{F}(\\mathbf{x}_0)\\neq\\mathbf{0}$.  However we can easily find the correction $\\boldsymbol{\\delta}$ as\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    " \\mathbf{F}(\\mathbf{x}_0+\\boldsymbol{\\delta})&= A\\mathbf{x_0} + A\\boldsymbol{\\delta} -\\mathbf{b} = \\mathbf{0}\\\\\n",
    "      &=A\\boldsymbol{\\delta} + \\mathbf{F}(\\mathbf{x}_0) = \\mathbf{0}\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "which we can solve immediately for the correction by solving the linear problem\n",
    "$$\n",
    "    A\\boldsymbol{\\delta} = -\\mathbf{F}(\\mathbf{x}_0)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which is the same as the general Newton problem if $A=J$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But it is!\n",
    "\n",
    "if \n",
    "$$\n",
    "    \\mathbf{F}(\\mathbf{x}) = A\\mathbf{x} - \\mathbf{b}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Then in index notation\n",
    "$$\n",
    "    F_i(\\mathbf{x}) = \\sum_{j} a_{ij}x_j - b_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "And the Jacobian\n",
    "$$\n",
    "    J_{ij} = \\frac{\\partial F_i}{\\partial x_j} = a_{ij}\n",
    "$$\n",
    "or $J = A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Or the solution of our original problem is \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\mathbf{x} &= \\mathbf{x}_0 + \\boldsymbol{\\delta} \\\\\n",
    "    &= \\mathbf{x}_0 - A^{-1}\\mathbf{F}(\\mathbf{x_0})\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which for a simple linear function $\\mathbf{F}(\\mathbf{x}) = A\\mathbf{x} - \\mathbf{b}$ becomes\n",
    "\n",
    "$$ \n",
    "\\begin{align*}\n",
    "    \\mathbf{x} &= \\mathbf{x}_0 - A^{-1}\\left(A\\mathbf{x}_0 - \\mathbf{b}\\right) \\\\\n",
    "                &= A^{-1}\\mathbf{b}\\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Thus we can solve for either the solution itself or the correction with the residual.  These are equivalent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Some quick and dirty Code\n",
    "def newton(F,J,x0,tol=1.e-6,MAX_ITS=100,verbose=True):\n",
    "    \"\"\" Solve F(x) = 0 using Newton's method until ||F(x)|| < tol or reaches Max iterations\n",
    "    \n",
    "    Params:\n",
    "    -------\n",
    "        F: calleable: \n",
    "            Function returning residual F(x)\n",
    "        J: calleable\n",
    "            Function returning Jacobian J(x)\n",
    "        tol: float\n",
    "            stopping criteria for ||F|| < tol\n",
    "        MAX_ITS: integer\n",
    "            maximum number of iterations\n",
    "        verbose: bool\n",
    "            if true spits out norm of the residual at each iteration\n",
    "            \n",
    "    Returns:\n",
    "    --------\n",
    "        x: list\n",
    "            list of points for each Newton iteration (just used for plotting intermediate results)\n",
    "            the solution is x[-1]\n",
    "    Raises:\n",
    "    -----------\n",
    "        ValueError if number of its exceeds MAX_ITS\n",
    "        \n",
    "    \"\"\"\n",
    "    x = [ x0 ]\n",
    "    for k in range(MAX_ITS+1):\n",
    "        xk = x[k]\n",
    "        res = numpy.linalg.norm(F(xk))\n",
    "        if verbose:\n",
    "            print('k = {}, ||F|| = {}'.format(k,res))\n",
    "        delta = numpy.linalg.solve(J(xk),-F(xk))\n",
    "        x.append( xk + delta)\n",
    "        \n",
    "        if res < tol: \n",
    "            return numpy.array(x), k\n",
    "                \n",
    "    raise ValueError('Maximum number of iterations exceeded {}'.format(MAX_ITS))\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Let's test with our toy problem\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gamma = numpy.array([10, -2])\n",
    "\n",
    "F = lambda x: numpy.array([ x[0]**2 + 2*x[1]**2 - gamma[0],\n",
    "                            x[0]    + 5*x[1]    - gamma[1] ])\n",
    "\n",
    "J = lambda x: numpy.array([[ 2*x[0],  4*x[1] ],\n",
    "                           [ 1    ,    5   ]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "x0 = numpy.array([-.5, 2])\n",
    "\n",
    "x, nits = newton(F,J,x0)\n",
    "print('\\n',x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "xa = numpy.linspace(-4,4,100)\n",
    "ya = numpy.linspace(-4,4,100)\n",
    "X,Y = numpy.meshgrid(xa,ya)\n",
    "\n",
    "Z1 = lambda gamma: X**2 + 2*Y**2 - gamma\n",
    "Z2 = lambda gamma: X + 5*Y - gamma\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.contour(X,Y,Z1(gamma[0]),[0],colors='r')\n",
    "axes.contour(X,Y,Z2(gamma[1]),[0],colors='b')\n",
    "for xi in x:    \n",
    "    axes.plot(xi[0],xi[1],'gx',markersize=8)\n",
    "axes.plot(xi[0],xi[1],'go',markersize=10)\n",
    "\n",
    "axes.axis('equal')\n",
    "\n",
    "axes.set_xlabel('$x_1$')\n",
    "axes.set_ylabel('$x_2$')\n",
    "axes.set_title('$\\gamma_1={}, \\gamma_2={}$'.format(gamma[0],gamma[1]))\n",
    "axes.grid()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A project idea:  \n",
    "\n",
    "Write an adaptive BDF-2/BDF-3  scheme that can solve a general non-linear problem and test it using the vander-pol oscillator.  This turns out to be a lot more challenging than it looks for the vander-pol problem."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}