09_ODE_ivp_part1.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": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Numerical Solution to ODE Initial Value Problems - Part 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Many physical, biological, and societal systems can be written as a system of ordinary differential equations (ODEs).  In the case where the initial state (value) is know the problems can be written as\n",
    "\n",
    "$$\n",
    "    \\frac{\\text{d} \\vec{\\!u}}{\\text{d}t} = \\vec{\\!f}(t, \\vec{\\!u}) \\quad \\vec{\\!u}(0) = \\vec{\\!u}_0\n",
    "$$\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$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples:  Simple radioactive decay\n",
    "$$\n",
    "    \\vec{\\!u} = [c]\n",
    "$$\n",
    "   \n",
    "$$\n",
    "    \\frac{\\text{d} c}{\\text{d}t} = \\lambda c \\quad c(0) = c_0\n",
    "$$\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": "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, 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_xlim((0.0, 1650))\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_3}{\\text{d}t} &= \\lambda_2 c_3 - \\lambda_3 c_3 \n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\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:  Particle tracking in a fluid\n",
    "\n",
    "$$\\frac{\\text{d} \\vec{\\!X}}{\\text{d}t} = \\vec{\\!V}(t, \\vec{\\!X})$$\n",
    "\n",
    "In fact all ODE IVP systems can be thought of as tracking particles through a flow field (dynamical system).  In 1-dimension the flow \"manifold\" we are on is fixed by the initial condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples: Van der Pol Oscillator\n",
    "\n",
    "$$y'' - \\mu (1 - y^2) y' + y = 0 \\quad \\quad \\text{with} \\quad \\quad  y(0) = y_0, \\quad 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",
    "# N = 100\n",
    "N = 500\n",
    "# N = 1000\n",
    "t = numpy.linspace(0.0, 100, N)\n",
    "u = numpy.empty((2, N))\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.plot(u[0,:], u[1, :], 'ro')\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": "slide"
    }
   },
   "source": [
    "## Basic Stepping Schemes\n",
    "\n",
    "Introducing some notation to simpify things\n",
    "$$\\begin{aligned}\n",
    "    t_0 &= 0 \\\\\n",
    "    t_1 &= t_0 + \\Delta t \\\\\n",
    "    t_n &= t_{n-1} + \\Delta t = n \\Delta t + t_0 \\\\\n",
    "    u_0 &= u(t_0) \\approx U_0 \\\\\n",
    "    u_1 &= u(t_1) \\approx U_1 \\\\\n",
    "    u_n &= u(t_n) \\approx U_2 \\\\\n",
    "\\end{aligned}$$\n",
    "where lower-case letters are \"exact\".  "
   ]
  },
  {
   "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",
    "\n",
    "We still need to decide how to evaluate the $f(t, u)$ term however.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let us look at this from a perspective of quadrature, take the integral of both sides:\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) &\\approx \\Delta t f(t, u(t))\n",
    "\\end{aligned}$$\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": [
    "We can rewrite our scheme\n",
    "$$\n",
    "    u(t + \\Delta t) - u(t) = \\Delta t f(t, u(t))\n",
    "$$\n",
    "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": {
    "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, 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": "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 Forward Euler\n",
    "t_euler = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t_euler[1] - t_euler[0]\n",
    "u_euler = numpy.empty(t_euler.shape)\n",
    "u_euler[0] = c_0\n",
    "for (n, t_n) in enumerate(t_euler[:-1]):\n",
    "    u_euler[n + 1] = u_euler[n] + delta_t * f(t_n, u_euler[n])\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t_euler, u_euler, '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",
    "N = 25\n",
    "\n",
    "# Stable example\n",
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "t_leapfrog = numpy.linspace(0.0, 1.6e3, 25)\n",
    "\n",
    "# Unstable example\n",
    "# decay_constant = -1.0\n",
    "# t_exact = numpy.linspace(0.0, 5.0, 100)\n",
    "# t_leapfrog = numpy.linspace(0.0, 5.0, N)\n",
    "\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "\n",
    "u_exact = c_0 * numpy.exp(decay_constant * t_exact)\n",
    "\n",
    "# Implement leap-frog\n",
    "delta_t = t_leapfrog[1] - t_leapfrog[0]\n",
    "u_leapfrog = numpy.empty(t_leapfrog.shape)\n",
    "u_leapfrog[0] = c_0\n",
    "u_leapfrog[1] = u_leapfrog[0] + delta_t * f(t_leapfrog[0], u_leapfrog[0])\n",
    "for n in range(1, t_leapfrog.shape[0] - 1):\n",
    "    u_leapfrog[n + 1] = u_leapfrog[n - 1] + 2.0 * delta_t * f(t_leapfrog[n], u_leapfrog[n])\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t_leapfrog, u_leapfrog, '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.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Similar to forward Euler is the *backward Euler* method which, as you may have guessed, evaluates the function $f$ at the updated time 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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For some cases we can solve the equation before hand.  For instance in the case of our example problem we have:\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t f(t_{n+1}, U_{n+1}) = U_n + \\Delta t (\\lambda U_{n+1})\n",
    "$$\n",
    "which can be solved for $U_{n+1}$ to find\n",
    "$$\\begin{aligned}\n",
    "    U_{n+1} &= U_n + \\Delta t (\\lambda U_{n+1}) \\\\\n",
    "    U_{n+1} \\left[ 1 - \\Delta t \\lambda \\right ] &= U_n \\\\\n",
    "    U_{n+1} &= \\frac{U_n}{1 - \\Delta t \\lambda}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "It's also useful to be able to do this in the case of systems of ODEs.  Let $f(U) = A U$, then\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    U_{n+1} &= U_n + \\Delta t (A U_{n+1}) \\\\\n",
    "    U_{n+1} \\left [ I - \\Delta t A \\right ] &= U_n \\\\\n",
    "    U_{n+1} &= \\left [ I - \\Delta t A \\right]^{-1} U_n\n",
    "\\end{aligned}$$\n",
    "\n",
    "In general however we are often not able to do this with arbitrary $f$."
   ]
  },
  {
   "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 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": "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_backwards = numpy.linspace(0.0, 1.6e3, 10)\n",
    "delta_t = t_backwards[1] - t_backwards[0]\n",
    "u_backwards = numpy.empty(t_backwards.shape)\n",
    "u_backwards[0] = c_0\n",
    "for n in range(0, t_backwards.shape[0] - 1):\n",
    "    u_backwards[n + 1] = u_backwards[n] / (1.0 - decay_constant * delta_t)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t_backwards, u_backwards, '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": [
    "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",
    "$$\n",
    "\n",
    "In this case what is the update scheme?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    U_{n+1} &= U_{n} + \\frac{\\Delta t}{2} (f(U_n) + f(U_{n+1})) \\\\\n",
    "    U_{n+1} &= U_{n} + \\frac{\\Delta t}{2} (\\lambda U_n + \\lambda U_{n+1}) \\\\\n",
    "    U_{n+1} \\left[1 - \\frac{\\Delta t \\lambda}{2}  \\right] &= U_{n} \\left[1 + \\frac{\\Delta t \\lambda}{2} \\right] \\\\\n",
    "    U_{n+1} &= U_{n} \\frac{1 + \\frac{\\Delta t \\lambda}{2}}{1 - \\frac{\\Delta t \\lambda}{2}} \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "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",
    "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 of ODE Methods\n",
    "\n",
    "At this point it is also helpful to introduce more notation to distinguish between the true solution to the ODE $u(t_n)$ and the approximated value which we will denote $U_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Definition:** We define the *truncation error* of a scheme by replacing the $U_n$ with the true solution $u(t_n)$ in the finite difference formula and looking at the difference from the exact solution.\n",
    "\n",
    "For example we will use the difference form of forward Euler\n",
    "$$\n",
    "    \\frac{U_{n+1} - U_n}{\\Delta t} = f(t_n)\n",
    "$$\n",
    "and define the truncation error as\n",
    "$$\n",
    "    T(t, u; \\Delta t) = \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - f(t_n, u(t_n)).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Definition:** A method is called *consistent* if \n",
    "$$\n",
    "    \\lim_{\\Delta t \\rightarrow 0} T(t, u; \\Delta t) = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Definition:** We say that a method is *order* $p$ accurate if\n",
    "\n",
    "$$\n",
    "    \\lVert T(t, u; \\Delta t) \\rVert \\leq C \\Delta t^p\n",
    "$$\n",
    "\n",
    "uniformally on $t \\in [0, T]$.  This can also be written as $T(t, u; \\Delta t) = \\mathcal{O}(\\Delta t^p)$.  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",
    "$$"
   ]
  },
  {
   "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_n + \\Delta t f(t_n, u_n) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the definition of truncation error we can use our Taylor series expression and find the truncation error.  Take the finite difference form of forward Euler\n",
    "$$\n",
    "    \\frac{U_{n+1} - U_n}{\\Delta t} = f(t_n)\n",
    "$$\n",
    "and replacing the derivative formulation with $u(t_n)$ to find\n",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - f(t_n) \\\\\n",
    "    &= \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - u'(t_n).\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From here we use the Taylor series centered at $t_n$ and evaluated at $t_{n+1}$ to find\n",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - u'(t_n) \\\\\n",
    "    &= \\frac{1}{\\Delta t} \\left[ u(t_n) + u'(t_n) (t - t_n) + \\frac{u''(t_n)}{2} (t - t_n)^2 + \\mathcal{O}((t-t_n)^3) - u(t_n) \\right] - u'(t_n) \\\\\n",
    "    &=  u'(t_n) + \\frac{u''(t_n)}{2} \\Delta t + \\mathcal{O}(\\Delta t^2) - u'(t_n) \\\\\n",
    "    &= \\frac{u''(t_n)}{2} \\Delta t + \\mathcal{O}(\\Delta t^2).\n",
    "\\end{aligned}$$\n",
    "\n",
    "This implies that forward Euler is first order accurate and therefore consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Another equivalent definition of the truncation error uses the form\n",
    "$$\n",
    "    U_{n+1} = u(t_n) + \\Delta t f(t_n)\n",
    "$$\n",
    "and the definition\n",
    "$$\n",
    "    T(t, u; \\Delta t) = \\frac{1}{\\Delta t} \\left [ U_{n+1} - u(t_{n+1}) \\right]\n",
    "$$\n",
    "to find\n",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{1}{\\Delta t} [U_{n+1} - u(t + \\Delta t)] \\\\\n",
    "    &= \\frac{1}{\\Delta t} \\left[ \\underbrace{u_n + \\Delta t f(t_n, u_n)}_{U_{n+1}} - \\underbrace{\\left( u_n + \\Delta t f(t_n, u_n) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3) \\right )}_{u(t_{n+1})}\\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}$$"
   ]
  },
  {
   "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:\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_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}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We need one more expansion however due to leap-frog.  Recall that leap-frog has the form\n",
    "$$\n",
    "    \\frac{U_{n+1} - U_{n-1}}{2 \\Delta t} = f(t_n, U_n)\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    U_{n+1} = U_{n-1} + 2 \\Delta t f(t_n, U_n).\n",
    "$$\n",
    "To handle the $U_{n-1}$ term we need to write this with relation to $u(t_n)$.  Again we use the Taylor series\n",
    "$$\n",
    "    u(t_{n-1}) = u_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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    u(t_{n+1}) &= u_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_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",
    "\n",
    "Plugging these into our definition of the truncation error along with the leap-frog method definition leads to\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{1}{\\Delta t} \\left [\\underbrace{U_{n-1} + 2 \\Delta t f_n}_{U_{n+1}} - \\underbrace{\\left(u_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 )}_{u(t + \\Delta t)} \\right ] \\\\\n",
    "    &=\\frac{1}{\\Delta t} \\left [\\underbrace{\\left(u_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)}_{U_{n-1}} + 2\\Delta t f_n - \\left(u_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",
    "    &=\\frac{1}{\\Delta t} \\left [- \\Delta t^3 \\frac{u'''(t_n)}{3} + \\mathcal{O}(\\Delta t^4) \\right ] \\\\\n",
    "    &=- \\Delta t^2 \\frac{u'''(t_n)}{3} + \\mathcal{O}(\\Delta t^3)\n",
    "\\end{aligned}$$\n",
    "Therefore the method is second order accurate and is consistent theoretically.  In practice it's a bit more complicated than that."
   ]
  },
  {
   "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_euler = numpy.empty(len(num_steps))\n",
    "error_trap = numpy.empty(len(num_steps))\n",
    "error_leapfrog = 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] = u_0\n",
    "    for n in range(t.shape[0] - 1):\n",
    "        u_euler[n+1] = u_euler[n] + delta_t[i] * f(t[n], u_euler[n])\n",
    "        \n",
    "    # Compute trapezoidal\n",
    "    u_trap = numpy.empty(t.shape)\n",
    "    u_trap[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_trap[n + 1] = u_trap[n] * integration_constant\n",
    "        \n",
    "    # Compute Leap-Frog\n",
    "    u_leapfrog = numpy.empty(t.shape)\n",
    "    u_leapfrog[0] = 1.0\n",
    "    u_leapfrog[1] = u_euler[1]\n",
    "    for n in range(1, t.shape[0] - 1):\n",
    "        u_leapfrog[n+1] = u_leapfrog[n-1] + 2.0 * delta_t[i] * f(t[n], u_leapfrog[n])\n",
    "        \n",
    "    # Compute error for each\n",
    "    error_euler[i] = numpy.linalg.norm(delta_t[i] * (u_euler - u_exact(t)), ord=1)\n",
    "    error_trap[i] = numpy.linalg.norm(delta_t[i] * (u_trap - u_exact(t)), ord=1)\n",
    "    error_leapfrog[i] = numpy.linalg.norm(delta_t[i] * (u_leapfrog - 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",
    "axes.loglog(delta_t, error_euler, 'bo', label='Forward Euler')\n",
    "axes.loglog(delta_t, error_trap, 'ro', label='Trapezoidal')\n",
    "axes.loglog(delta_t, error_leapfrog, 'go', label=\"Leap-Frog\")\n",
    "\n",
    "axes.loglog(delta_t, order_C(delta_t[2], error_euler[2], 1.0) * delta_t**1.0, '--b')\n",
    "axes.loglog(delta_t, order_C(delta_t[2], error_trap[2], 2.0) * delta_t**2.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[2], error_leapfrog[2], 2.0) * delta_t**2.0, '--r')\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": "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{align*}\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{align*}$$\n",
    "which contains derivatives of $u$ up to $p$th order.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "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": [
    "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",
    "$$\n",
    "\n",
    "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)$ and therefore\n",
    "$$\\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)) ~~~ \\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^*) \\\\\n",
    "    &= 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",
    "$$\\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",
    "    T(t, u; \\Delta t) &= \\frac{1}{\\Delta t} \\left[u_n + \\Delta t f\\left(u_n + \\frac{1}{2} \\Delta t f(u_n)\\right) - \\left(u_n + \\Delta t f(t_n, u_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",
    "\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",
    "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",
    "\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, u_2[n+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",
    "    \\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",
    "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",
    "All these methods have $\\alpha_r = 1$, $\\alpha_{r-1} = -1$ and $\\alpha_j=0$ for $j < r - 1$."
   ]
  },
  {
   "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": {
    "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 = 20\n",
    "# N = 10\n",
    "t = numpy.linspace(0, 10.0, N)\n",
    "delta_t = t[1] - t[0]\n",
    "u_ab2 = numpy.empty(t.shape)\n",
    "\n",
    "# Use RK-2 to start the method\n",
    "u_ab2[0] = 1.0\n",
    "u_ab2[1] = u_ab2[0] + 0.5 * delta_t * f(t[0], u_ab2[0])\n",
    "u_ab2[1] = u_ab2[0] + delta_t * f(t[0], u_ab2[1])\n",
    "for n in range(0,len(t)-2):\n",
    "    u_ab2[n+2] = u_ab2[n + 1] + delta_t / 2.0 * (-f(t[n], u_ab2[n]) + 3.0 * f(t[n+1], u_ab2[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_ab2, 'ro', label=\"2-step A-B\")\n",
    "\n",
    "axes.set_title(\"Adams-Bashforth Method\")\n",
    "axes.set_xlabel(\"t\")\n",
    "axes.set_ylabel(\"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",
    "# N = 10\n",
    "# N = 5\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_ylabel(\"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",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{1}{\\Delta t} \\left [\\sum^r_{j=0} \\alpha_j u_{n+j} - \\Delta t \\sum^r_{j=0} \\beta_j f(u_{n+j}, 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",
    "    T(t, u; \\Delta t) &= \\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 (\\sum^r_{j=0} \\left(\\frac{1}{q!} j^q \\alpha_j - \\frac{1}{(q-1)!} j^{q-1} \\beta_j \\right) \\right) u^{(q)}(t_n) + \\cdots\n",
    "\\end{aligned}$$\n",
    "\n",
    "The method is *consistent* if the first two terms of the expansion vanish, i.e. $\\sum^r_{j=0} \\alpha_j = 0$ and $\\sum^r_{j=0} j \\alpha_j = \\sum^r_{j=0} \\beta_j$."
   ]
  },
  {
   "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, 'ko', label='RK-2')\n",
    "axes.loglog(delta_t, error_ab, 'bo', 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[0], error_ab[0], 1.0) * delta_t**1.0, '--r')\n",
    "axes.loglog(delta_t, order_C(delta_t[1], error_ab[1], 2.0) * delta_t**2.0, '--b')\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",
    "$$\\begin{aligned}\n",
    "    \\hat{U}_{n+1} &= U_n + \\Delta t f(U_n) \\\\\n",
    "    U_{n+1} &= U_n + \\frac{1}{2} \\Delta t (f(U_n) + f(\\hat{U}_{n+1}) \n",
    "\\end{aligned}$$\n",
    "leading to a second order accurate method."
   ]
  },
  {
   "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_ylabel(\"u(t)\")\n",
    "axes.legend(loc=1)\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.1"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}