07_ivp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Methods for Initial Value Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We now turn towards time dependent PDEs.  Before moving to the full PDEs we will explore numerical methods for systems of ODEs that are initial value problems of the general form\n",
    "$$\n",
    "    \\frac{\\text{d} \\vec{u}}{\\text{d}t} = \\vec{f}(t, \\vec{u}) \\quad \\vec{u}(0) = \\vec{u}_0\n",
    "$$\n",
    "where\n",
    " - $\\vec{u}(t)$ is the state vector\n",
    " - $\\vec{f}(t, \\vec{u})$ is a vector-valued function that controls the growth of $\\vec{u}$ with time\n",
    " - $\\vec{u}(0)$ is the initial condition at time $t = 0$\n",
    "Note that the right hand side function $f$ could in actuality be the discretization in space of a PDE, i.e. a system of equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples:  Simple radioactive decay\n",
    "$\\vec{u} = [c]$\n",
    "   \n",
    "$$\\frac{\\text{d} c}{\\text{d}t} = -\\lambda c \\quad c(0) = c_0$$\n",
    "   \n",
    "\n",
    "which has solutions of the form $c(t) = c_0 e^{-\\lambda t}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, 1.0 * numpy.exp(-decay_constant * t))\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples:  Complex radioactive decay (or chemical system).\n",
    "\n",
    "Chain of decays from one species to another.\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{\\text{d} c_1}{\\text{d}t} &= -\\lambda_1 c_1 \\\\\n",
    "    \\frac{\\text{d} c_2}{\\text{d}t} &= \\lambda_1 c_1 - \\lambda_2 c_2 \\\\\n",
    "    \\frac{\\text{d} c_2}{\\text{d}t} &= \\lambda_2 c_3 - \\lambda_3 c_3 \n",
    "\\end{aligned}$$\n",
    "\n",
    "$$\\frac{\\text{d} \\vec{u}}{\\text{d}t} = \\frac{\\text{d}}{\\text{d}t}\\begin{bmatrix} c_1 \\\\ c_2 \\\\ c_3 \\end{bmatrix} = \n",
    "\\begin{bmatrix} \n",
    "    -\\lambda_1 & 0 & 0 \\\\\n",
    "    \\lambda_1 & -\\lambda_2 & 0 \\\\\n",
    "    0 & \\lambda_2 & -\\lambda_3\n",
    "\\end{bmatrix} \\begin{bmatrix} c_1 \\\\ c_2 \\\\ c_3 \\end{bmatrix}$$\n",
    "\n",
    "$$\\frac{\\text{d} \\vec{u}}{\\text{d}t} = A \\vec{u}$$\n",
    "\n",
    "For systems of equations like this the general solution to the ODE is the matrix exponential:\n",
    "\n",
    "$$\\vec{u}(t) = \\vec{u}_0 e^{A t}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples: Van der Pol Oscillator\n",
    "\n",
    "$$y'' - \\mu (1 - y^2) y' + y = 0~~~~~\\text{with}~~~~ y(0) = y_0, ~~~y'(0) = v_0$$\n",
    " \n",
    "$$\\vec{u} = \\begin{bmatrix} y \\\\ y' \\end{bmatrix} = \\begin{bmatrix} u_1 \\\\ u_2 \\end{bmatrix}$$\n",
    "   \n",
    "$$\\frac{\\text{d}}{\\text{d}t} \\begin{bmatrix} u_1 \\\\ u_2 \\end{bmatrix} = \\begin{bmatrix} u_2 \\\\ \\mu (1 - u_1^2) u_2 - u_1 \\end{bmatrix} = \\vec{f}(t, \\vec{u})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "\n",
    "def f(t, u, mu=5):\n",
    "    return numpy.array([u[1], mu * (1.0 - u[0]**2) * u[1] - u[0]])\n",
    "\n",
    "t = numpy.linspace(0.0, 100, 1000)\n",
    "u = numpy.empty((2, t.shape[0]))\n",
    "u[:, 0] = [0.1, 0.0]\n",
    "\n",
    "integrator = integrate.ode(f)\n",
    "integrator.set_integrator(\"dopri5\")\n",
    "integrator.set_initial_value(u[:, 0])\n",
    "\n",
    "for (n, t_n) in enumerate(t[1:]):\n",
    "    integrator.integrate(t_n)\n",
    "    if not integrator.successful():\n",
    "        break\n",
    "    u[:, n + 1] = integrator.y\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t, u[0,:])\n",
    "axes.set_title(\"Solution to Van der Pol Oscillator\")\n",
    "axes.set_xlabel(\"t\")\n",
    "axes.set_ylabel(\"y(t)\")\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(u[0,:], u[1, :])\n",
    "axes.set_title(\"Phase Diagram for Van der Pol Oscillator\")\n",
    "axes.set_xlabel(\"y(t)\")\n",
    "axes.set_ylabel(\"y'(t)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples:  Heat Equation\n",
    "\n",
    "Let's try to construct a system of ODEs that represents the heat equation\n",
    "$$\n",
    "    u_t = u_{xx}.\n",
    "$$\n",
    "If we discretize the right hand side with second order, centered differences with $m$ points we would have\n",
    "$$\n",
    "    \\frac{\\text{d}}{\\text{d} t} U_i(t) = \\frac{U_{i+1}(t) - 2 U_i(t) + U_{i-1}(t)}{\\Delta x^2}\n",
    "$$\n",
    "where we now have $m$ unknown, time dependent functions to solve for.  This approach to discretizing a PDE is sometimes called a method-of-lines approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence and Uniqueness of Solutions \n",
    "\n",
    "One important step before diving into the numerical methods for IVP ODE problems is to understand what the behavior of the solutions are, whether they exist, and if they might be unique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Linear Systems\n",
    "\n",
    "For linear ODEs we have the generic system \n",
    "$$\n",
    "    u'(t) = A(t) u + g(t)\n",
    "$$\n",
    "where $A$ is time-dependent matrix and $g$ a vector.  Note that linear systems always have a unique solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $g(t) = 0$ for all $t$ we say the ODE is *homogeneous* and the matrix $A$ is time independent (then implying that it is also *autonomous*) then the solution to this ODE is \n",
    "$$\n",
    "    u(t) = u(t_0) e^{A(t - t_0)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the case where $g(t) \\neq 0$ for all $t$ the ODE is *inhomogeneous* we can use Duhamel's principle which tells us\n",
    "$$\n",
    "    u(t) = u(t_0) e^{A(t-t_0)} + \\int^t_{t_0} e^{A(t - \\tau)} g(\\tau) d\\tau.\n",
    "$$\n",
    "We can think of the operator $e^{A(t-\\tau)}$ as the solution operator for the homogeneous ODE which can map the solution at time $\\tau$ to the solution at time $t$ giving this form of the solution a Green's function type property."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Non-linear Existance and Uniqueness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Lipschitz Continuity\n",
    "\n",
    "Generalizing uniqueness to non-linear ODEs requires a special type of continuity called *Lipschitz continuity*.  Consider the  ODE\n",
    "$$\n",
    "    u'(t) = f(u,t),  \\quad \\quad \\quad u(t_0) = u_0,\n",
    "$$\n",
    "we will require a certain amount of smoothness in the right hand side function $f(u,t)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We say that $f$ is Lipshitz continuous in $u$ over some domain\n",
    "$$\n",
    "    \\Omega = \\{(u,t) : |u - u_0| \\leq a, t_0 \\leq t \\leq t_1 \\}\n",
    "$$\n",
    "if there exists a constant $L > 0$ such that\n",
    "$$\n",
    "    |f(u,t) - f(u^\\ast, t)| \\leq L |u - u^\\ast| \\quad \\quad \\forall (u,t) ~\\text{and}~ (u^\\ast,t) \\in \\Omega.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If $f(u,t)$ is differentiable with respect to $u$ in $\\Omega$, i.e. the Jacobian $f_u = \\partial f / \\partial u$ exists, and is bounded then we can say\n",
    "$$\n",
    "    L = \\max_{(u,t) \\in \\Omega} |f_u(u,t)|.\n",
    "$$\n",
    "We can use this bound since\n",
    "$$\n",
    "    f(u,t) = f(u^\\ast, t) + f_u(v,t)(u-u^\\ast)\n",
    "$$\n",
    "for some $v$ chosen to be in-between $u$ and $u^\\ast$ which is effectively the Taylor series error bound and implies smoothness of $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "With Lipshitz continuity of $f$ we can guarantee a unique solution the IVP at least to time $T = \\min(t_1, t_0 + a/S)$ where \n",
    "$$\n",
    "    S = \\max_{(u,t)\\in\\Omega} |f(u,t)|.\n",
    "$$\n",
    "This value $S$ is the modulus of the maximum slope that the solution $u(t)$ can obtain in $\\Omega$ and guarantees that we remain in $\\Omega$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Consider $u'(t) = (u(t))^2, u(0) = u_0 > 0$.  If we define our domain of interest as above we can compute the Lipshitz constant as\n",
    "$$\n",
    "    L = \\max_{(u,t) \\in \\Omega} | 2 u | = 2 (u_0 + a)\n",
    "$$\n",
    "where we have used the restriction from $\\Omega$ that $|u - u_0| \\leq a$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly we can compute $S$ to find\n",
    "$$\n",
    "    S = \\max_{(u,t)\\in\\Omega} |f(u,t)| = (u_0 + a)^2\n",
    "$$\n",
    "so that we can guarantee a unique solution up until $T  = a / (u_0 + a)^2$.  Given that we can choose $a$ we can simply choose a value that maximized $T$, in this case $a = u_0$ does this and we conclude that we have a unique solution up until $T = 1 / 4 u_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Since we also know the exact solution to the ODE above,\n",
    "$$\n",
    "    u(t) = \\frac{1}{1/u_0 - t},\n",
    "$$\n",
    "we can see that $|u(t)| < \\infty$ as long as $t \\neq 1/u_0$.  Note that once we reach the pole in the denominator there is no longer a solution possible for the IVP past this point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "Consider the IVP\n",
    "$$\n",
    "    u' = \\sqrt{u} \\quad \\quad u(0) = 0.\n",
    "$$\n",
    "\n",
    "Where is this $f$ Lipshitz continuous?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Computing the derivative we find\n",
    "$$\n",
    "    f_u = \\frac{1}{2\\sqrt{u}}\n",
    "$$\n",
    "which goes to infinity as $u \\rightarrow 0$.  We can therefore not guarantee a unique solution near the given initial condition.  In fact we know this as the ODE has two solutions\n",
    "$$\n",
    "    u(t) = 0 \\quad \\text{and} \\quad u(t) = \\frac{1}{4} t^2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Systems of Equations\n",
    "\n",
    "A similar notion for Lipschitz continuity exists in a particular norm $||\\cdot||$ if there is a constant $L$ such that\n",
    "$$\n",
    "    ||f(u,t) - f(u^\\ast,t)|| \\leq L ||u - u^\\ast||\n",
    "$$\n",
    "for all $(u,t)$ and $(u^\\ast,t)$ in the domain $\\Omega = \\{(u,t) : ||u-u_0|| \\leq a, t_0 \\leq t \\leq t_1 \\}$.  Note that if the function $f$ is Lipschitz continuous in one norm it is continuous in any norm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Basic Stepping Schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Looking back at our work on numerical differentiation why not approximate the derivative as a finite difference:\n",
    "$$\n",
    "    \\frac{u(t + \\Delta t) - u(t)}{\\Delta t} = f(t, u)\n",
    "$$\n",
    "We still need to decide how to evaluate the $f(t, u)$ term however.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Lets look at this from a perspective of quadrature, take the integral of both sides:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\int^{t + \\Delta t}_t \\frac{\\text{d} u}{\\text{d}\\tilde{t}} d\\tilde{t} &= \\int^{t + \\Delta t}_t f(t, u) d\\tilde{t} \\\\ ~ \\\\\n",
    "u(t + \\Delta t) - u(t) &= \\Delta t ~f(t, u(t)) \\\\ ~ \\\\\n",
    "\\frac{u(t + \\Delta t) - u(t)}{\\Delta t} &= f(t, u(t))\n",
    "\\end{aligned}$$\n",
    "\n",
    "where we have used a left-sided quadrature rule for the integral on the right.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Introducing some notation to simplify things\n",
    "$$\n",
    "    t_0 = 0  \\quad \\quad t_1 = t_0 + \\Delta t  \\quad \\quad t_n = t_{n-1} + \\Delta t = n \\Delta t + t_0\n",
    "$$\n",
    "\n",
    "$$\n",
    "    U^0 = u(t_0) \\quad \\quad U^1 = u(t_1) \\quad \\quad U^n = u(t_n)\n",
    "$$\n",
    "we can rewrite our scheme as\n",
    "$$\n",
    "    \\frac{U^{n+1} - U^n}{\\Delta t} = f(t_n, U^n)\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    U^{n+1} = U^n + \\Delta t f(t_n, U^n)\n",
    "$$\n",
    "which is known as the *forward Euler method*.  In essence we are approximating the derivative with the value of the function at the point we are at $t_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(-decay_constant * t), label=\"True Solution\")\n",
    "\n",
    "# Plot Euler step\n",
    "dt = 1e3\n",
    "u_np = c_0 + dt * (-decay_constant * c_0)\n",
    "axes.plot((0.0, dt), (c_0, u_np), 'k')\n",
    "axes.plot((dt, dt), (u_np, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, 0.0), (c_0, u_np), 'k--')\n",
    "axes.plot((0.0, dt), (u_np, u_np), 'k--')\n",
    "axes.text(400, u_np - 0.05, '$\\Delta t$', fontsize=16)\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_xlim(-1e2, 1.6e3)\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Note where we expect error due to the approximation and how it manifests in the example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: -decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = c_0 * numpy.exp(-decay_constant * t_exact)\n",
    "\n",
    "# Implement Euler\n",
    "t = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = c_0\n",
    "for (n, t_n) in enumerate(t[:-1]):\n",
    "    U[n + 1] = U[n] + delta_t * f(t_n, U[n])\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U, 'or', label=\"Euler\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Forward Euler\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_xlabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A similar method can be derived if we consider instead using the second order accurate central difference:\n",
    "\n",
    "$$\\frac{U^{n+1} - U^{n-1}}{2\\Delta t} = f(t_{n}, U^{n})$$\n",
    "\n",
    "this method is known as the leap-frog method.  Note that the way we have written this method requires a previous function evaluation and technically is a \"multi-step\" method although we do not actually use the current evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(-decay_constant * t), label=\"True Solution\")\n",
    "\n",
    "# Plot Leap-Frog step\n",
    "dt = 1e3\n",
    "u_np = c_0 + dt * (-decay_constant * c_0 * numpy.exp(-decay_constant * dt / 2.0))\n",
    "axes.plot((0.0, dt), (c_0, u_np), 'k')\n",
    "axes.plot((dt, dt), (u_np, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, 0.0), (c_0, u_np), 'k--')\n",
    "axes.plot((0.0, dt), (u_np, u_np), 'k--')\n",
    "axes.text(400, u_np - 0.05, '$\\Delta t$', fontsize=16)\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_xlim(-1e2, 1.6e3)\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: -decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = c_0 * numpy.exp(-decay_constant * t_exact)\n",
    "\n",
    "# Implement leap-frog\n",
    "t = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = c_0\n",
    "# First evaluation use Euler to get us going\n",
    "U[1] = U[0] + delta_t * f(t[0], U[0])\n",
    "for n in range(1, t.shape[0] - 1):\n",
    "    U[n + 1] = U[n - 1] + 2.0 * delta_t * f(t[n], U[n])\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U, 'or', label=\"Leap-Frog\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Leap-Frog\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_xlabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similar to forward Euler is the *backward Euler* method which evaluates the function $f$ at the updated time (right hand quadrature rule) so that\n",
    "$$\n",
    "    U^{n+1} = U^n + \\Delta t f(t_{n+1}, U^{n+1}).\n",
    "$$\n",
    "Schemes where the function $f$ is evaluated at the unknown time are called *implicit methods*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(-decay_constant * t), label=\"True Solution\")\n",
    "\n",
    "# Plot Euler step\n",
    "dt = 1e3\n",
    "u_np = c_0 + dt * (-decay_constant * c_0 * numpy.exp(-decay_constant * dt))\n",
    "axes.plot((0.0, dt), (c_0, u_np), 'k')\n",
    "axes.plot((dt, dt), (u_np, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, 0.0), (c_0, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, dt), (c_0 * numpy.exp(-decay_constant * dt), c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.text(400, u_np - 0.05, '$\\Delta t$', fontsize=16)\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_xlim(-1e2, 1.6e3)\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Write code that implements the backward Euler method\n",
    "$$\n",
    "    U^{n+1} = U^n + \\Delta t f(t_{n+1}, U^{n+1}).\n",
    "$$\n",
    "\n",
    "Note in the following what bias the error tends to have and relate this back to the approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: -decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = c_0 * numpy.exp(-decay_constant * t_exact)\n",
    "\n",
    "# Implement backwards Euler\n",
    "t = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = c_0\n",
    "for n in range(0, t.shape[0] - 1):\n",
    "    U[n + 1] = U[n] / (1.0 + decay_constant * delta_t)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U, 'or', label=\"Backwards Euler\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Backwards Euler\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_xlabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A modification on the Euler methods involves using the approximated midpoint to evaluate $f(t, u)$ called the midpoint method.  The scheme is\n",
    "$$\n",
    "    \\frac{U^{n+1} - U^{n}}{\\Delta t} = f\\left(\\frac{U^n + U^{n+1}}{2} \\right).\n",
    "$$\n",
    "This is the simplest example of a *symplectic integrator* which has special properties for integrating Hamiltonian systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(-decay_constant * t), label=\"True Solution\")\n",
    "\n",
    "# Plot Midpoint step\n",
    "dt = 1e3\n",
    "u_np = c_0 * (1.0 - decay_constant * dt / 2.0) / (1.0 + decay_constant * dt / 2.0)\n",
    "axes.plot((0.0, dt), (c_0, u_np), 'k')\n",
    "axes.plot((dt, dt), (u_np, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, 0.0), (c_0, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, dt), (c_0 * numpy.exp(-decay_constant * dt), c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.text(400, u_np - 0.05, '$\\Delta t$', fontsize=16)\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_xlim(-1e2, 1.6e3)\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Implement the midpoint method\n",
    "$$\n",
    "    \\frac{U^{n+1} - U^{n}}{\\Delta t} = f\\left(\\frac{U^n + U^{n+1}}{2} \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: -decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = c_0 * numpy.exp(-decay_constant * t_exact)\n",
    "\n",
    "# Implement midpoint\n",
    "t = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = c_0\n",
    "integration_constant = (1.0 - decay_constant * delta_t / 2.0) / (1.0 + decay_constant * delta_t / 2.0)\n",
    "for n in range(0, t.shape[0] - 1):\n",
    "    U[n + 1] = U[n] * integration_constant\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U, 'or', label=\"Midpoint Method\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Midpoint Method\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_xlabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Another simple implicit method is based on integration using the trapezoidal method.  The scheme is\n",
    "$$\n",
    "    \\frac{U^{n+1} - U^{n}}{\\Delta t} = \\frac{1}{2} (f(U^n) + f(U^{n+1}))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0.0, 1.6e3, 100)\n",
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(-decay_constant * t), label=\"True Solution\")\n",
    "\n",
    "# Plot Trapezoidal step\n",
    "dt = 1e3\n",
    "u_np = c_0 * (1.0 - decay_constant * dt / 2.0) / (1.0 + decay_constant * dt / 2.0)\n",
    "axes.plot((0.0, dt), (c_0, u_np), 'k')\n",
    "axes.plot((dt, dt), (u_np, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, 0.0), (c_0, c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.plot((0.0, dt), (c_0 * numpy.exp(-decay_constant * dt), c_0 * numpy.exp(-decay_constant * dt)), 'k--')\n",
    "axes.text(400, u_np - 0.05, '$\\Delta t$', fontsize=16)\n",
    "\n",
    "axes.set_title(\"Radioactive Decay with $t_{1/2} = 1600$ years\")\n",
    "axes.set_xlabel('t (years)')\n",
    "axes.set_ylabel('$c$')\n",
    "axes.set_xlim(-1e2, 1.6e3)\n",
    "axes.set_ylim((0.5,1.0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Again implement the trapezoidal method\n",
    "$$\n",
    "    \\frac{U^{n+1} - U^{n}}{\\Delta t} = \\frac{1}{2} (f(U^n) + f(U^{n+1}))\n",
    "$$\n",
    "What is this method equivalent to and why?  Is this generally true?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "c_0 = 1.0\n",
    "decay_constant = numpy.log(2.0) / 1600.0\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = c_0 * numpy.exp(-decay_constant * t_exact)\n",
    "\n",
    "# Implement trapezoidal method\n",
    "t = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = c_0\n",
    "integration_constant = (1.0 - decay_constant * delta_t / 2.0) / (1.0 + decay_constant * delta_t / 2.0)\n",
    "for n in range(t.shape[0] - 1):\n",
    "    U[n + 1] = U[n] * integration_constant\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U, 'or', label=\"Trapezoidal\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Trapezoidal\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_xlabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Error Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Truncation Errors\n",
    "We can define truncation errors the same as we did before where we insert the true solution of the ODE into the difference equation and use Taylor series expansions.  It is important that at this juncture we use the form of the difference that models the derivative directly as otherwise we will find something different.\n",
    "\n",
    "Define the finite difference approximation to the derivative as $D(U^{n+1}, U^n, U^{n-1}, \\ldots)$ and write the schemes above as\n",
    "\n",
    "$$\n",
    "    D(U^{n+1}, U^n, U^{n-1}, \\ldots) = F(t^{n+1}, t^n, t^{n-1}, \\ldots, U^{n+1}, U^n, U^{n-1}, \\ldots)\n",
    "$$\n",
    "\n",
    "where $F(\\cdot)$ now has some relation to evaluations of the function $f(t, u)$.  Then the local truncation error can be defined as\n",
    "\n",
    "$$\n",
    "    \\tau^n = D(u(t^{n+1}), u(t^n), u(t^{n-1}), \\ldots) - F(t^{n+1}, t^n, t^{n-1}, \\ldots, u(t^{n+1}), u(t^n), u(t^{n-1}), \\ldots)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similarly if we know\n",
    "$$\n",
    "    \\lim_{\\Delta t \\rightarrow 0} \\tau^n = 0\n",
    "$$\n",
    "then the discretized equation is considered consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Order of accuracy is also defined the same way as before.  If\n",
    "$$\n",
    "    || \\tau || \\leq C \\Delta t^p\n",
    "$$\n",
    "uniformly on $t \\in [0, T]$ then the discretization is $p$th order accurate.  Note that a method is consistent if $p > 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Error Analysis of Forward Euler\n",
    "\n",
    "We can analyze the error and convergence order of forward Euler by considering the Taylor series centered at $t_n$:\n",
    "$$\n",
    "    u(t) = u(t_n) + (t - t_n) u'(t_n) + \\frac{u''(t_n)}{2} (t - t_n)^2 + \\mathcal{O}((t-t_n)^3)\n",
    "$$\n",
    "\n",
    "Try to compute the LTE for forward Euler's method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Evaluating this series at $t_{n+1}$ gives\n",
    "$$\\begin{aligned}\n",
    "u(t_{n+1}) &= u(t_n) + (t_{n+1} - t_n) u'(t_n) + \\frac{u''(t_n)}{2} (t_{n+1} - t_n)^2 + \\mathcal{O}((t_{n+1}-t_n)^3)\\\\\n",
    "&=u(t_n) + \\Delta t f(t_n, u(t_n)) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "From the definition of truncation error we can use our Taylor series expression and find the truncation error to be\n",
    "$$\\begin{aligned}\n",
    "    \\tau^n &= \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - f(t_n, u(t_n)) \\\\\n",
    "    &= \\frac{1}{\\Delta t} \\left [u(t_n) + \\Delta t ~ f(t_n, u(t_n)) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3) - u(t_n) - \\Delta t ~ f(t_n, u(t_n)) \\right ]\\\\\n",
    "    &= \\frac{1}{\\Delta t} \\left [ \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3) \\right ] \\\\\n",
    "    &= \\frac{u''(t_n)}{2} \\Delta t + \\mathcal{O}(\\Delta t^2)\n",
    "\\end{aligned}$$\n",
    "This implies that forwar Euler is first order accurate and therefore consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Error Analysis of Leap-Frog Method\n",
    "\n",
    "To easily analyze this method we will expand the Taylor series from before to another order and evaluate at both the needed positions:\n",
    "$$\n",
    "    u(t) = u(t_n) + (t - t_n) u'(t_n) + (t - t_n)^2 \\frac{u''(t_n)}{2}  + (t - t_n)^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}((t-t_n)^4)\n",
    "$$\n",
    "leading to \n",
    "$$\\begin{aligned}\n",
    "u(t_{n+1}) &= u(t_n) + \\Delta t f_n + \\Delta t^2 \\frac{u''(t_n)}{2}  + \\Delta t^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^4)\\\\\n",
    "u(t_{n-1}) &= u(t_n) - \\Delta t f_n + \\Delta t^2 \\frac{u''(t_n)}{2}  - \\Delta t^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^4)\n",
    "\\end{aligned}$$\n",
    "See if you can compute the LTE in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Plugging this into our definition of the truncation error along with the leap-frog method definition leads to\n",
    "$$\\begin{aligned}\n",
    "    \\tau^n &= \\frac{u(t_{n+1}) - u(t_{n-1})}{2 \\Delta t} - f(t_n, u(t_n)) \\\\\n",
    "    &=\\frac{1}{\\Delta t} \\left[\\frac{1}{2}\\left( u(t_n) + \\Delta t f_n + \\Delta t^2 \\frac{u''(t_n)}{2}  + \\Delta t^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^4)\\right) \\right . \\\\\n",
    "    &\\quad \\quad\\left . - \\frac{1}{2} \\left ( u(t_n) - \\Delta t f_n + \\Delta t^2 \\frac{u''(t_n)}{2}  - \\Delta t^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^4)\\right ) - \\Delta t~ f(t_n, u(t_n)) \\right ] \\\\\n",
    "    &= \\frac{1}{\\Delta t} \\left [\\Delta t^3 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^5)\\right ] \\\\\n",
    "    &= \\Delta t^2 \\frac{u'''(t_n)}{6} + \\mathcal{O}(\\Delta t^4)\n",
    "\\end{aligned}$$\n",
    "Therefore the method is second order accurate and is consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler and Leap-Frog\n",
    "f = lambda t, u: -u\n",
    "u_exact = lambda t: numpy.exp(-t)\n",
    "u_0 = 1.0\n",
    "\n",
    "t_f = 10.0\n",
    "num_steps = [2**n for n in range(4,10)]\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error = numpy.empty((5, len(num_steps)))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    \n",
    "    # Note that in the cases below we can instantiate this array now\n",
    "    # rather than every time as none of the implicit methods require\n",
    "    # the space to store the future solution\n",
    "    U = numpy.empty(t.shape)\n",
    "    \n",
    "    # Compute ForwardEuler solution\n",
    "    U[0] = u_0\n",
    "    for n in range(t.shape[0] - 1):\n",
    "        U[n+1] = U[n] + delta_t[i] * f(t[n], U[n])\n",
    "    error[0, i] = numpy.linalg.norm(delta_t[i] * (U - u_exact(t)), ord=1)\n",
    "        \n",
    "    # Compute Leap-Frog\n",
    "    U[0] = u_0\n",
    "    U[1] = U[0] + delta_t[i] * f(t[0], U[0])\n",
    "    for n in range(1, t.shape[0] - 1):\n",
    "        U[n+1] = U[n-1] + 2.0 * delta_t[i] * f(t[n], U[n])\n",
    "    error[1, i] = numpy.linalg.norm(delta_t[i] * (U - u_exact(t)), ord=1)\n",
    "    \n",
    "    # Compute Backward Euler\n",
    "    U[0] = u_0\n",
    "    for n in range(0, t.shape[0] - 1):\n",
    "        U[n + 1] = U[n] / (1.0 + delta_t[i])\n",
    "    error[2, i] = numpy.linalg.norm(delta_t[i] * (U - u_exact(t)), ord=1)\n",
    "    \n",
    "    # Compute mid-pointU[0] = c_0\n",
    "    U[0] = u_0\n",
    "    integration_constant = (1.0 - delta_t[i] / 2.0) / (1.0 + delta_t[i] / 2.0)\n",
    "    for n in range(0, t.shape[0] - 1):\n",
    "        U[n + 1] = U[n] * integration_constant\n",
    "    error[3, i] = numpy.linalg.norm(delta_t[i] * (U - u_exact(t)), ord=1)\n",
    "    \n",
    "    # Compute trapezoidal\n",
    "    U[0] = u_0\n",
    "    integration_constant = (1.0 - delta_t[i] / 2.0) / (1.0 + delta_t[i] / 2.0)\n",
    "    for n in range(t.shape[0] - 1):\n",
    "        U[n + 1] = U[n] * integration_constant\n",
    "    error[4, i] = numpy.linalg.norm(delta_t[i] * (U - u_exact(t)), ord=1)\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "style = ['bo', 'go', 'ro', 'cs', 'yo']\n",
    "label = ['Forward Euler', \"Leap-Frog\", \"Backward Euler\", \"Mid-Point\", \"Trapezoidal\"]\n",
    "order = [1, 2, 1, 2, 2]\n",
    "\n",
    "for k in range(5):\n",
    "    axes.loglog(delta_t, error[k, :], style[k], label=label[k])\n",
    "    axes.loglog(delta_t, order_C(delta_t[2], error[k, 2], order[k]) * delta_t**order[k], 'k--')\n",
    "\n",
    "# axes.legend(loc=2)\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### One-Step Errors\n",
    "\n",
    "There is another definition of local truncation error sometimes used in ODE numerical methods called the *one-step error* which is slightly different than our local truncation error definition.  Our definition uses the direct discretization of the derivatives to find the LTE where as this alternative bases the error on a form that looks like it is updating the previous value.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As an example consider the leap-frog method, the LTE we found before was based on\n",
    "$$\n",
    "    \\frac{U_{n+1} - U_{n-1}}{2 \\Delta t} = f(U_n)\n",
    "$$\n",
    "leading us to a second order LTE.  \n",
    "\n",
    "For the one-step error we consider instead \n",
    "$$\n",
    "    U_{n+1} = U_{n-1} + 2 \\Delta t f(U_n)\n",
    "$$\n",
    "which leads to the one-step error $\\mathcal{O}(\\Delta t^3)$ instead!\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\mathcal{L}^n &= u(t_{n+1}) - u(t_{n-1}) - 2 \\Delta t f(u(t_n)) \\\\\n",
    "    &= \\frac{1}{3} \\Delta t^3 u'''(t_n) + \\mathcal{O}(\\Delta t^5) \\\\\n",
    "    &= 2 ~\\Delta t ~\\tau^n.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This one-step error is suggestively named to indicate that perhaps this is the error for one time step where as the global error may be higher.  To remain consistent with our previous discussion of convergences we will continue to use our previous definition of the LTE.  We will show that with the appropriate definition of stability and a $p$ order LTE we can expect a $p$th order global error.  In general for a $p+1$th order one-step error the global error will be $p$th order."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Taylor Series Methods\n",
    "\n",
    "A **Taylor series method** can be derived by direct substitution of the right-hand-side function $f(t, u)$ and it's appropriate derivatives into the Taylor series expansion for $u(t_{n+1})$.  For a $p$th order method we would look at the Taylor series up to that order and replace all the derivatives of $u$ with derivatives of $f$ instead.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the general case we have\n",
    "$$\\begin{aligned}\n",
    "    u(t_{n+1}) = u(t_n) + \\Delta t u'(t_n) + \\frac{\\Delta t^2}{2} u''(t_n) + \\frac{\\Delta t^3}{6} u'''(t_n) + \\cdots + \\frac{\\Delta t^p}{p!} u^{(p)}(t_n)\n",
    "\\end{aligned}$$\n",
    "which contains derivatives of $u$ up to $p$th order.  We then replace these derivatives with the appropriate derivative of $f$ which will always be one less than the derivative of $u$ (due to the original ODE)\n",
    "$$\n",
    "    u^{(p)}(t_n) = f^{(p-1)}(t_n, u(t_n))\n",
    "$$\n",
    "leading to the method\n",
    "$$\n",
    "    u(t_{n+1}) = u(t_n) + \\Delta t f(t_n, u(t_n)) + \\frac{\\Delta t^2}{2} f'(t_n, u(t_n)) + \\frac{\\Delta t^3}{6} f''(t_n, u(t_n)) + \\cdots + \\frac{\\Delta t^p}{p!} f^{(p-1)}(t_n, u(t_n)).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The drawback to these methods is that we have to derive a new one each time we have a new $f$ and we also need $p-1$ derivatives of $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### 2nd Order Taylor Series Method\n",
    "\n",
    "We want terms up to second order so we need to take the derivative of $u' = f(t, u)$ once to find $u'' = f'(t, u)$.  See if you can derive the method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\\begin{align*}\n",
    "    u(t_{n+1}) &= u(t_n) + \\Delta t u'(t_n) + \\frac{\\Delta t^2}{2} u''(t_n) \\\\\n",
    "    &=u(t_n) + \\Delta t f(t_n, u(t_n)) + \\frac{\\Delta t^2}{2} f'(t_n, u(t_n)) \\quad \\text{or} \\\\\n",
    "    U^{n+1} &= U^n + \\Delta t f(t_n, U^n) + \\frac{\\Delta t^2}{2} f'(t_n, U^n).\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Runge-Kutta Methods\n",
    "\n",
    "One way to derive higher-order ODE solvers is by computing intermediate stages.  These are not *multi-step* methods as they still only require information from the current time step but they raise the order of accuracy by adding *stages*.  These types of methods are called **Runge-Kutta** methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Two-stage Runge-Kutta Methods\n",
    "\n",
    "The basic idea behind the simplest of the Runge-Kutta methods is to approximate the solution at $t_n + \\Delta t / 2$ via Euler's method and use this in the function evaluation for the final update.\n",
    "$$\\begin{aligned}\n",
    "    U^* &= U^n + \\frac{1}{2} \\Delta t f(U^n) \\\\\n",
    "    U^{n+1} &= U^n + \\Delta t f(U^*) = U^n + \\Delta t f(U^n + \\frac{1}{2} \\Delta t f(U^n))\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The truncation error can be computed similarly to how we did so before but we do need to figure out how to compute the derivative inside of the function.  Note that due to \n",
    "$$\n",
    "    f(u(t_n)) = u'(t_n)\n",
    "$$ \n",
    "that differentiating this leads to \n",
    "$$\n",
    "    f'(u(t_n)) u'(t_n) = u''(t_n)\n",
    "$$\n",
    "leading to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f\\left(u(t_n) + \\frac{1}{2} \\Delta t f(u(t_n)) \\right ) &= f\\left(u(t_n) +\\frac{1}{2} \\Delta t u'(t_n) \\right ) \\\\\n",
    "    &= f(u(t_n)) + \\frac{1}{2} \\Delta t u'(t_n) f'(u(t_n)) + \\frac{1}{8} \\Delta t^2 (u'(t_n))^2 f''(u(t_n)) + \\mathcal{O}(\\Delta t^3) \\\\\n",
    "    &=u'(t_n) + \\frac{1}{2} \\Delta t u''(t_n) + \\mathcal{O}(\\Delta t^2)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Going back to the truncation error we have\n",
    "$$\\begin{aligned}\n",
    "    \\tau^n &= \\frac{1}{\\Delta t} \\left[u(t_n) + \\Delta t f\\left(u(t_n) + \\frac{1}{2} \\Delta t f(u(t_n))\\right) - \\left(u(t_n) + \\Delta t f(t_n, u(t_n)) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3) \\right ) \\right] \\\\\n",
    "    &=\\frac{1}{\\Delta t} \\left[\\Delta t u'(t_n) + \\frac{1}{2} \\Delta t^2 u''(t_n) + \\mathcal{O}(\\Delta t^3) - \\Delta t u'(t_n) - \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3) \\right] \\\\\n",
    "    &= \\mathcal{O}(\\Delta t^2)\n",
    "\\end{aligned}$$\n",
    "so this method is second order accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  4-stage Runge-Kutta Method\n",
    "\n",
    "$\\begin{aligned}\n",
    "    Y_1 &= U^n \\\\\n",
    "    Y_2 &= U^n + \\frac{1}{2} \\Delta t f(Y_1, t_n) \\\\\n",
    "    Y_3 &= U^n + \\frac{1}{2} \\Delta t f(Y_2, t_n + \\Delta t / 2) \\\\\n",
    "    Y_4 &= U^n + \\Delta t f(Y_3, t_n + \\Delta t / 2) \\\\\n",
    "    U^{n+1} &= U^n + \\frac{\\Delta t}{6} \\left [f(Y_1, t_n) + 2 f(Y_2, t_n + \\Delta t / 2) + 2 f(Y_3, t_n + \\Delta t/2) + f(Y_4, t_n + \\Delta t) \\right ]\n",
    "\\end{aligned}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Implement and compare the two-stage and 4-stage Runge-Kutta methods\n",
    "f = lambda t, u: -u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 10.0, 100)\n",
    "u_exact = numpy.exp(-t_exact)\n",
    "\n",
    "N = 10\n",
    "t = numpy.linspace(0, 10.0, N)\n",
    "delta_t = t[1] - t[0]\n",
    "\n",
    "# RK 2\n",
    "U_2 = numpy.empty(t.shape)\n",
    "U_2[0] = 1.0\n",
    "for (n, t_n) in enumerate(t[1:]):\n",
    "    U_2[n+1] = U_2[n] + 0.5 * delta_t * f(t_n, U_2[n])\n",
    "    U_2[n+1] = U_2[n] + delta_t * f(t_n + 0.5 * delta_t, U_2[n+1])\n",
    "\n",
    "# RK4\n",
    "U_4 = numpy.empty(t.shape)\n",
    "U_4[0] = 1.0\n",
    "for (n, t_n) in enumerate(t[1:]):\n",
    "    y_1 = U_4[n]\n",
    "    y_2 = U_4[n] + 0.5 * delta_t * f(t_n, y_1)\n",
    "    y_3 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5 * delta_t, y_2)\n",
    "    y_4 = U_4[n] + delta_t * f(t_n + 0.5 * delta_t, y_3)\n",
    "    U_4[n+1] = U_4[n] + delta_t / 6.0 * (f(t_n, y_1) + 2.0 * f(t_n + 0.5 * delta_t, y_2) + 2.0 * f(t_n + 0.5 * delta_t, y_3) + f(t_n + delta_t, y_4))\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t_exact, u_exact, 'k', label=\"True\")\n",
    "axes.plot(t, U_2, 'ro', label=\"2-Stage\")\n",
    "axes.plot(t, U_4, 'bo', label=\"4-Stage\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler and RK\n",
    "f = lambda t, u: -u\n",
    "u_exact = lambda t: numpy.exp(-t)\n",
    "\n",
    "t_f = 10.0\n",
    "num_steps = [2**n for n in range(5,12)]\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_euler = numpy.empty(len(num_steps))\n",
    "error_2 = numpy.empty(len(num_steps))\n",
    "error_4 = numpy.empty(len(num_steps))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    \n",
    "    # Compute Euler solution\n",
    "    U_euler = numpy.empty(t.shape)\n",
    "    U_euler[0] = 1.0\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U_euler[n+1] = U_euler[n] + delta_t[i] * f(t_n, U_euler[n])\n",
    "        \n",
    "    # Compute 2 and 4-stage\n",
    "    U_2 = numpy.empty(t.shape)\n",
    "    U_4 = numpy.empty(t.shape)\n",
    "    U_2[0] = 1.0\n",
    "    U_4[0] = 1.0\n",
    "    for (n, t_n) in enumerate(t[1:]):\n",
    "        U_2[n+1] = U_2[n] + 0.5 * delta_t[i] * f(t_n, U_2[n])\n",
    "        U_2[n+1] = U_2[n] + delta_t[i] * f(t_n, U_2[n+1])\n",
    "        y_1 = U_4[n]\n",
    "        y_2 = U_4[n] + 0.5 * delta_t[i] * f(t_n, y_1)\n",
    "        y_3 = U_4[n] + 0.5 * delta_t[i] * f(t_n + 0.5 * delta_t[i], y_2)\n",
    "        y_4 = U_4[n] + delta_t[i] * f(t_n + 0.5 * delta_t[i], y_3)\n",
    "        U_4[n+1] = U_4[n] + delta_t[i] / 6.0 * (f(t_n, y_1) + 2.0 * f(t_n + 0.5 * delta_t[i], y_2) + 2.0 * f(t_n + 0.5 * delta_t[i], y_3) + f(t_n + delta_t[i], y_4))\n",
    "        \n",
    "    # Compute error for each\n",
    "    error_euler[i] = numpy.abs(U_euler[-1] - u_exact(t_f)) / numpy.abs(u_exact(t_f))\n",
    "    error_2[i] = numpy.abs(U_2[-1] - u_exact(t_f)) / numpy.abs(u_exact(t_f))\n",
    "    error_4[i] = numpy.abs(U_4[-1] - u_exact(t_f)) / numpy.abs(u_exact(t_f))\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_t, error_euler, 'bo', label='Forward Euler')\n",
    "axes.loglog(delta_t, error_2, 'ro', label='2-stage')\n",
    "axes.loglog(delta_t, error_4, 'go', label=\"4-stage\")\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_euler[1], 1.0) * delta_t**1.0, '--b')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_2[1], 2.0) * delta_t**2.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_4[1], 4.0) * delta_t**4.0, '--g')\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear Multi-Step Methods\n",
    "\n",
    "Multi-step methods (as introduced via the leap-frog method) are ODE methods that require multiple time step evaluations to work.  Some of the advanatages of using a multi-step method rather than one-step method included\n",
    "\n",
    " - Taylor series methods require differentiating the given equation which can be cumbersome and difficult to impelent\n",
    " - One-step methods at higher order often require the evaluation of the function $f$ many times\n",
    " \n",
    "Disadvantages\n",
    "\n",
    " - Methods are not self-starting, i.e. they require other methods to find the initial values\n",
    " - The time step $\\Delta t$ in one-step methods can be changed at any time while multi-step methods this is much more complex\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### General Linear Multi-Step Methods\n",
    "\n",
    "All linear multi-step methods can be written as the linear combination of past, present and future solutions:\n",
    "\n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j U^{n+j} = \\Delta t \\sum^r_{j=0} \\beta_j f(U^{n+j}, t_{n+j})\n",
    "$$\n",
    "\n",
    "If $\\beta_r = 0$ then the method is explicit (only requires previous time steps).  Note that the coefficients are not unique as we can multiply both sides by a constant.  In practice a normalization of $\\alpha_r = 1$ is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: Adams Methods\n",
    "\n",
    "$$\n",
    "    U^{n+r} = U^{n+r-1} + \\Delta t \\sum^r_{j=0} \\beta_j f(U^{n+j})\n",
    "$$\n",
    "\n",
    "All these methods have $\\alpha_r = 1$, $\\alpha_{r-1} = -1$ and $\\alpha_j=0$ for $j < r - 1$ leaving the method to be specified by how the evaluations of $f$ is done determining the $\\beta_j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Adams-Bashforth Methods\n",
    "The **Adams-Bashforth** methods are explicit solvers that maximize the order of accuracy given a number of steps $r$.  This is accomplished by looking at the Taylor series and picking the coefficients $\\beta_j$ to elliminate as many terms in the Taylor series as possible.\n",
    "$$\\begin{aligned}\n",
    "    \\text{1-step:} & & U_{n+1} &= U_n +\\Delta t f(U_n) \\\\\n",
    "    \\text{2-step:} & & U_{n+2} &= U_{n+1} + \\frac{\\Delta t}{2} (-f(U_n) + 3 f(U_{n+1})) \\\\\n",
    "    \\text{3-step:} & & U_{n+3} &= U_{n+2} + \\frac{\\Delta t}{12} (5 f(U_n) - 16 f(U_{n+1}) + 23 f(U_{n+2})) \\\\\n",
    "    \\text{4-step:} & & U_{n+4} &= U_{n+3} + \\frac{\\Delta t}{24} (-9 f(U_n) + 37 f(U_{n+1}) -59 f(U_{n+2}) + 55 f(U_{n+3}))\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "run_control": {
     "marked": false
    },
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Use 2-step Adams-Bashforth to compute solution\n",
    "f = lambda t, u: -u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 10.0, 100)\n",
    "u_exact = numpy.exp(-t_exact)\n",
    "\n",
    "N = 50\n",
    "t = numpy.linspace(0, 10.0, N)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "\n",
    "# Use RK-2 to start the method\n",
    "U[0] = 1.0\n",
    "U[1] = U[0] + 0.5 * delta_t * f(t[0], U[0])\n",
    "U[1] = U[0] + delta_t * f(t[0], U[1])\n",
    "for n in range(0,len(t)-2):\n",
    "    U[n+2] = U[n + 1] + delta_t / 2.0 * (-f(t[n], U[n]) + 3.0 * f(t[n+1], U[n+1]))\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t_exact, u_exact, 'k', label=\"True\")\n",
    "axes.plot(t, U, 'ro', label=\"2-step A-B\")\n",
    "\n",
    "axes.set_title(\"Adams-Bashforth Method\")\n",
    "axes.set_xlabel(\"t\")\n",
    "axes.set_xlabel(\"u(t)\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Adams-Moulton Methods\n",
    "The **Adams-Moulton** methods are the implicit versions of the Adams-Bashforth methods.  Since this gives one additional parameter to use $\\beta_r$ these methods are generally one order of accuracy greater than their counterparts.\n",
    "$$\\begin{aligned}\n",
    "    \\text{1-step:} &  & U_{n+1} &= U_n + \\frac{\\Delta t}{2} (f(U_n) + f(U_{n+1})) \\\\\n",
    "    \\text{2-step:} &  & U_{n+2} &= U_{n+1} + \\frac{\\Delta t}{12} (-f(U_n) + 8f(U_{n+1}) + 5f(U_{n+2})) \\\\\n",
    "    \\text{3-step:} &  & U_{n+3} &= U_{n+2} + \\frac{\\Delta t}{24} (f(U_n) - 5f(U_{n+1}) + 19f(U_{n+2}) + 9f(U_{n+3})) \\\\\n",
    "    \\text{4-step:} &  & U_{n+4} &= U_{n+3} + \\frac{\\Delta t}{720}(-19 f(U_n) + 106 f(U_{n+1}) -264 f(U_{n+2}) + 646 f(U_{n+3}) + 251 f(U_{n+4}))\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Use 2-step Adams-Moulton to compute solution\n",
    "# u' = - decay u\n",
    "decay_constant = 1.0\n",
    "f = lambda t, u: -decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 10.0, 100)\n",
    "u_exact = numpy.exp(-t_exact)\n",
    "\n",
    "N = 20\n",
    "t = numpy.linspace(0, 10.0, N)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "U[0] = 1.0\n",
    "U[1] = U[0] + 0.5 * delta_t * f(t[0], U[0])\n",
    "U[1] = U[0] + delta_t * f(t[0], U[1])    \n",
    "integration_constant = 1.0 / (1.0 + 5.0 * decay_constant * delta_t / 12.0)\n",
    "for n in range(t.shape[0] - 2):\n",
    "    U[n+2] = (U[n+1] + decay_constant * delta_t / 12.0 * (U[n] - 8.0 * U[n+1])) * integration_constant\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t_exact, u_exact, 'k', label=\"True\")\n",
    "axes.plot(t, U, 'ro', label=\"2-step A-M\")\n",
    "\n",
    "axes.set_title(\"Adams-Moulton Method\")\n",
    "axes.set_xlabel(\"t\")\n",
    "axes.set_xlabel(\"u(t)\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Truncation Error for Multi-Step Methods\n",
    "\n",
    "We can again find the truncation error in general for linear multi-step methods:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\tau^n &= \\frac{1}{\\Delta t} \\left [\\sum^r_{j=0} \\alpha_j u(t_{n+j}) - \\Delta t \\sum^r_{j=0} \\beta_j f(t_{n+j}, u(t_{n+j})) \\right ]\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the general expansion and evalution of the Taylor series about $t_n$ we have\n",
    "$$\\begin{aligned}\n",
    "    u(t_{n+j}) &= u(t_n) + j \\Delta t u'(t_n) + \\frac{1}{2} (j \\Delta t)^2 u''(t_n) + \\mathcal{O}(\\Delta t^3) \\\\\n",
    "    u'(t_{n+j}) &= u'(t_n) + j \\Delta t u''(t_n) + \\frac{1}{2} (j \\Delta t)^2 u'''(t_n) + \\mathcal{O}(\\Delta t^3)\n",
    "\\end{aligned}$$\n",
    "leading to\n",
    "$$\\begin{aligned}\n",
    "\\tau^n &= \\frac{1}{\\Delta t}\\left( \\sum^r_{j=0} \\alpha_j\\right) u(t_{n}) + \\left(\\sum^r_{j=0} (j\\alpha_j - \\beta_j)\\right) u'(t_n) + \\Delta t \\left(\\sum^r_{j=0} \\left (\\frac{1}{2}j^2 \\alpha_j - j \\beta_j \\right) \\right) u''(t_n) \\\\\n",
    "&\\quad \\quad + \\cdots + \\Delta t^{q - 1} \\left (\\frac{1}{q!} \\left(j^q \\alpha_j - \\frac{1}{(q-1)!} j^{q-1} \\beta_j \\right) \\right) u^{(q)}(t_n) + \\cdots\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The method is *consistent* if the first two terms of the expansion vanish, i.e. \n",
    "$$\n",
    "    \\sum^r_{j=0} \\alpha_j = 0\n",
    "$$ \n",
    "and \n",
    "$$\n",
    "    \\sum^r_{j=0} j \\alpha_j = \\sum^r_{j=0} \\beta_j.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between RK-2, AB-2 and AM-2\n",
    "f = lambda t, u: -u\n",
    "u_exact = lambda t: numpy.exp(-t)\n",
    "\n",
    "t_f = 10.0\n",
    "num_steps = [2**n for n in range(4,10)]\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_rk = numpy.empty(len(num_steps))\n",
    "error_ab = numpy.empty(len(num_steps))\n",
    "error_am = numpy.empty(len(num_steps))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "        \n",
    "    # Compute RK2\n",
    "    U_rk = numpy.empty(t.shape)\n",
    "    U_rk[0] = 1.0\n",
    "    for n in range(t.shape[0]-1):\n",
    "        U_rk[n+1] = U_rk[n] + 0.5 * delta_t[i] * f(t[n], U_rk[n])\n",
    "        U_rk[n+1] = U_rk[n] + delta_t[i] * f(t[n], U_rk[n+1])\n",
    "        \n",
    "    # Compute Adams-Bashforth 2-stage\n",
    "    U_ab = numpy.empty(t.shape)\n",
    "    U_ab[:2] = U_rk[:2]\n",
    "    for n in range(t.shape[0] - 2):\n",
    "        U_ab[n+2] = U_ab[n + 1] + delta_t[i] / 2.0 * (-f(t[n], U_ab[n]) + 3.0 * f(t[n+1], U_ab[n+1]))\n",
    "    \n",
    "    # Compute Adama-Moulton 2-stage\n",
    "    U_am = numpy.empty(t.shape)\n",
    "    U_am[:2] = U_rk[:2]\n",
    "    decay_constant = 1.0\n",
    "    integration_constant = 1.0 / (1.0 + 5.0 * decay_constant * delta_t[i] / 12.0)\n",
    "    for n in range(t.shape[0] - 2):\n",
    "        U_am[n+2] = (U_am[n+1] + decay_constant * delta_t[i] / 12.0 * (U_am[n] - 8.0 * U_am[n+1])) * integration_constant\n",
    "        \n",
    "    # Compute error for each\n",
    "    error_rk[i] = numpy.linalg.norm(delta_t[i] * (U_rk - u_exact(t)), ord=1)\n",
    "    error_ab[i] = numpy.linalg.norm(delta_t[i] * (U_ab - u_exact(t)), ord=1)\n",
    "    error_am[i] = numpy.linalg.norm(delta_t[i] * (U_am - u_exact(t)), ord=1)\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_t, error_rk, 'bo', label='RK-2')\n",
    "axes.loglog(delta_t, error_ab, 'ro', label='AB-2')\n",
    "axes.loglog(delta_t, error_am, 'go', label=\"AM-2\")\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_rk[1], 2.0) * delta_t**2.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_ab[1], 2.0) * delta_t**2.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_am[1], 3.0) * delta_t**3.0, '--g')\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t) - u(t)|$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Predictor-Corrector Methods\n",
    "\n",
    "One way to simplify the Adams-Moulton methods so that implicit evaluations are not needed is by estimating the required implicit function evaluations with an explicit method.  These are often called **predictor-corrector** methods as the explicit method provides a *prediction* of what the solution might be and the not explicit *corrector* step works to make that estimate more accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example: One-Step Adams-Bashforth-Moulton\n",
    "\n",
    "Use the One-step Adams-Bashforth method to predict the value of $U^{n+1}$ and then use the Adams-Moulton method to correct that value:\n",
    "\n",
    "$\\hat{U}^{n+1} = U^n + \\Delta t f(U^n)$\n",
    "\n",
    "$U^{n+1} = U^n + \\frac{1}{2} \\Delta t (f(U^n) + f(\\hat{U}^{n+1})$\n",
    "\n",
    "This method is second order accurate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# One-step Adams-Bashforth-Moulton\n",
    "f = lambda t, u: -u\n",
    "\n",
    "t_exact = numpy.linspace(0.0, 10.0, 100)\n",
    "u_exact = numpy.exp(-t_exact)\n",
    "\n",
    "N = 100\n",
    "t = numpy.linspace(0, 10.0, N)\n",
    "delta_t = t[1] - t[0]\n",
    "U = numpy.empty(t.shape)\n",
    "\n",
    "U[0] = 1.0\n",
    "for n in range(t.shape[0] - 1):\n",
    "    U[n+1] = U[n] + delta_t * f(t[n], U[n])\n",
    "    U[n+1] = U[n] + 0.5 * delta_t * (f(t[n], U[n]) + f(t[n+1], U[n+1]))\n",
    "    \n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(t_exact, u_exact, 'k', label=\"True\")\n",
    "axes.plot(t, U, 'ro', label=\"2-step A-B\")\n",
    "\n",
    "axes.set_title(\"Adams-Bashforth-Moulton P/C Method\")\n",
    "axes.set_xlabel(\"t\")\n",
    "axes.set_xlabel(\"u(t)\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between RK-2, AB-2 and AM-2\n",
    "f = lambda t, u: -u\n",
    "u_exact = lambda t: numpy.exp(-t)\n",
    "\n",
    "t_f = 10.0\n",
    "num_steps = [2**n for n in range(4,10)]\n",
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_ab = numpy.empty(len(num_steps))\n",
    "error_am = numpy.empty(len(num_steps))\n",
    "error_pc = numpy.empty(len(num_steps))\n",
    "\n",
    "for (i, N) in enumerate(num_steps):\n",
    "    t = numpy.linspace(0, t_f, N)\n",
    "    delta_t[i] = t[1] - t[0]\n",
    "    \n",
    "    # RK-2 bootstrap for AB and AM\n",
    "    U_rk = numpy.empty(2)\n",
    "    U_rk[0] = 1.0\n",
    "    U_rk[1] = U_rk[0] + 0.5 * delta_t[i] * f(t[0], U_rk[0])\n",
    "    U_rk[1] = U_rk[0] + delta_t[i] * f(t[0], U_rk[1])\n",
    "        \n",
    "    # Compute Adams-Bashforth 2-stage\n",
    "    U_ab = numpy.empty(t.shape)\n",
    "    U_ab[:2] = U_rk[:2]\n",
    "    for n in range(t.shape[0] - 2):\n",
    "        U_ab[n+2] = U_ab[n + 1] + delta_t[i] / 2.0 * (-f(t[n], U_ab[n]) + 3.0 * f(t[n+1], U_ab[n+1]))\n",
    "    \n",
    "    # Compute Adams-Moulton 2-stage\n",
    "    U_am = numpy.empty(t.shape)\n",
    "    U_am[:2] = U_ab[:2]\n",
    "    decay_constant = 1.0\n",
    "    integration_constant = 1.0 / (1.0 + 5.0 * decay_constant * delta_t[i] / 12.0)\n",
    "    for n in range(t.shape[0] - 2):\n",
    "        U_am[n+2] = (U_am[n+1] + decay_constant * delta_t[i] / 12.0 * (U_am[n] - 8.0 * U_am[n+1])) * integration_constant\n",
    "        \n",
    "    # Compute Adams-Bashforth-Moulton\n",
    "    U_pc = numpy.empty(t.shape)\n",
    "    U_pc[0] = 1.0\n",
    "    for n in range(t.shape[0] - 1):\n",
    "        U_pc[n+1] = U_pc[n] + delta_t[i] * f(t[n], U_pc[n])\n",
    "        U_pc[n+1] = U_pc[n] + 0.5 * delta_t[i] * (f(t[n], U_pc[n]) + f(t[n+1], U_pc[n+1]))\n",
    "        \n",
    "    # Compute error for each\n",
    "    error_ab[i] = numpy.linalg.norm(delta_t[i] * (U_ab - u_exact(t)), ord=1)\n",
    "    error_am[i] = numpy.linalg.norm(delta_t[i] * (U_am - u_exact(t)), ord=1)\n",
    "    error_pc[i] = numpy.linalg.norm(delta_t[i] * (U_pc - u_exact(t)), ord=1)\n",
    "    \n",
    "# Plot error vs. delta_t\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.loglog(delta_t, error_pc, 'bo', label='PC')\n",
    "axes.loglog(delta_t, error_ab, 'ro', label='AB-2')\n",
    "axes.loglog(delta_t, error_am, 'go', label=\"AM-2\")\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_pc[1], 2.0) * delta_t**2.0, '--b')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_ab[1], 2.0) * delta_t**2.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_am[1], 3.0) * delta_t**3.0, '--g')\n",
    "\n",
    "axes.legend(loc=4)\n",
    "axes.set_title(\"Comparison of Errors\")\n",
    "axes.set_xlabel(\"$\\Delta t$\")\n",
    "axes.set_ylabel(\"$|U(t) - u(t)|$\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}