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": {
    "hide_input": true,
    "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}\\mathbf{u}}{\\text{d}t} = \\mathbf{f}(t, \\mathbf{u}) \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$\n",
    "\n",
    "where\n",
    " - $\\mathbf{\\!u}(t)$ is the state vector\n",
    " - $\\mathbf{\\!f}(t, \\mathbf{\\!u})$ is a vector-valued function that describes the evolution of $\\mathbf{u}$ with time\n",
    " - $\\mathbf{\\!u}(0)$ is the initial condition at time $t = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Example for our time:  A non-linear  model of Epidemics\n",
    "\n",
    "Classical [Kermack and McKendrick (1927)](https://royalsocietypublishing.org/doi/10.1098/rspa.1927.0118) SIR model of epidemics (with reinfection)\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\frac{ds}{dt} &= -si + kr \\\\\n",
    "    \\frac{di}{dt} &= si -\\sigma i \\\\\n",
    "    \\frac{dr}{dt} &= \\sigma i - kr\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Where the variable $s$ represents the fraction of a population that is **Susceptible** to infection, $i$ is the proportion **Infected** and $r$ the fraction **Recovered**.  For this model $s+i+r =1$.  The parameters $\\sigma, k \\geq 0$ control the relative rates of infection and recovery.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For this problem\n",
    "$$\n",
    "\\mathbf{u}(t) = \\begin{bmatrix} s(t)\\\\ i(t)\\\\ r(t)\\\\\\end{bmatrix},\\quad\\quad\n",
    "\\mathbf{f}(t,\\mathbf{u}) = \\begin{bmatrix} -si + kr \\\\ si -\\sigma i \\\\ \\sigma i - kr\\\\\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Numerical Solutions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Solve using SciPy's ODE integrator solve_ivp\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "# define the RHS of our system of ODE's\n",
    "def f_sir(t, u, sigma, k):\n",
    "    s,i,r = u\n",
    "    return numpy.array([-s*i + k*r,\n",
    "                        (s - sigma)*i,\n",
    "                        sigma*i - k*r ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sigma = .5\n",
    "k = 0.025\n",
    "t_max = 100\n",
    "u_0 = [0.999, 0.001, 0.]\n",
    "sol = solve_ivp(f_sir, [0, t_max] , u_0, args=(sigma, k), rtol=1.e-6, atol=1.e-9,dense_output = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0, t_max, 300)\n",
    "z = sol.sol(t)\n",
    "\n",
    "fig = plt.figure(figsize=(20,7))\n",
    "axes = fig.add_subplot(1,2,1)\n",
    "axes.plot(t,z[0],'r',label='s', linewidth=2)\n",
    "axes.plot(t,z[1],'b',label='i', linewidth=2)\n",
    "axes.plot(t,z[2],'g',label='r',  linewidth=2)\n",
    "axes.plot(t,sigma*numpy.ones(t.shape),'k--',label='$\\sigma$')\n",
    "\n",
    "axes.legend(loc='best',shadow=True, fontsize=14)\n",
    "axes.set_xlabel('Time',fontsize=16)\n",
    "axes.set_ylabel('Population',fontsize=16)\n",
    "axes.grid()\n",
    "axes.set_title('SIR system: $\\sigma={}$, $k={}$'.format(sigma,k),fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Questions an epidemiologist might ask:\n",
    " - What are the dynamics of this system?  does it predict steady or oscillatory solutions?\n",
    " - Can we estimate critical parameters ($\\sigma$, $k$) from data?\n",
    " - Can we reliably use this model to predict the future?\n",
    " - How do we evaluate whether this is a *useful* model?\n",
    " - How might I modify/improve this model.\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " \n",
    "### Questions a Computational Mathematician might ask:\n",
    " - Does a solution to the model even exist and is it unique?\n",
    " - Is our approximate numerical solution accurate?\n",
    " - What are the dynamics of this system?  does it predict steady or oscillatory solutions?\n",
    " - how do we understand the sensitivity to parameters?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Existence and Uniqueness of solutions (n-D autonomous systems)\n",
    "\n",
    "For proof see [Hirsch, Smale, Devaney, Dynamical Systems](https://www.amazon.com/Differential-Equations-Dynamical-Systems-Introduction/dp/0123820103)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Theorem: (Picard-Lindelhof)\n",
    "\n",
    "Given an Autonomous, dynamical system\n",
    "\n",
    "$$\n",
    "    \\frac{\\text{d}\\mathbf{u}}{\\text{d}t} = \\mathbf{f}(\\mathbf{u}) \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$\n",
    "\n",
    "with $\\mathbf{u}\\in\\mathbb{R}^n$ and $\\mathbf{f}:\\mathbb{R}^n\\rightarrow\\mathbb{R}^n$\n",
    "\n",
    "Consider a \"spherical\" domain of radius $\\rho$ around the initial condition $\\mathbf{u}_0$.  \n",
    "\n",
    "If, within this domain,  $\\mathbf{f}$ is \n",
    "\n",
    "* Bounded: $|\\mathbf{f}| < M$\n",
    "* Lipshitz Continuous: $$|\\mathbf{f}(\\mathbf{x}) - \\mathbf{f}(\\mathbf{y})| < K|\\mathbf{x}- \\mathbf{y} |$$\n",
    "\n",
    "Then a unique solution exists to the ODE IVP for some interval of time $t\\in[-a,a]$ where $0 < a <  \\min(\\rho/M, 1/K)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Geometric Picture\n",
    "\n",
    "<table>\n",
    "    <tr align=center><td><img align=left src=\"./images/Picard-Lindelhof-figure.png\" width=800></td>\n",
    "        <td>$$\n",
    "    \\frac{\\text{d}\\mathbf{u}}{\\text{d}t} = \\mathbf{f}(\\mathbf{u}), \\quad \\mathbf{u}(0) = \\mathbf{u}_0\n",
    "$$</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* **Short Version**: If $\\mathbf{f}$ is sufficiently smooth, then a local solution to the ODE exists and is unique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* **Caveat**: The theorem itself gives *NO* constructive way to find that solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Other Examples:  Simple radioactive decay\n",
    "$$\n",
    "    \\mathbf{\\!u} = [c], \\quad \\mathbf{f} = [-\\lambda 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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.)/1600.\n",
    "t = numpy.linspace(0., 4800, 11)\n",
    "y = numpy.linspace(0., 1.2, 11)\n",
    "T, Y = numpy.meshgrid(t,y)\n",
    "dt = numpy.ones(T.shape)\n",
    "dy = -dt*Y\n",
    "\n",
    "tp = numpy.linspace(0., 4800, 100)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.quiver(T,Y, dt,dy,linewidth=0.1,color='gray')\n",
    "axes.plot(tp,numpy.exp(decay_constant*tp),linewidth=3)\n",
    "axes.plot(0.,1.,'ro', markersize=10)\n",
    "axes.plot(1600., 0.5,'rx',markersize=10)\n",
    "axes.grid()\n",
    "axes.set_title(\"Radioactive Decay, $u' = - \\lambda u$, $u(0)=1$, $t_{1/2}=1600$ yr\", fontsize=18)\n",
    "axes.set_xlabel('t (years)', fontsize=16)\n",
    "axes.set_ylabel('u', fontsize=16)\n",
    "\n",
    "axes.set_ylim((-.1,1.2))\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "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_2 - \\lambda_3 c_3 \n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\frac{\\text{d} \\mathbf{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} \\mathbf{u}}{\\text{d}t} = A \\mathbf{u}$$\n",
    "\n",
    "For systems of equations like this the general solution to the ODE is the matrix exponential:\n",
    "\n",
    "$$\\mathbf{u}(t) = e^{A t}\\mathbf{u}_0 = \\sum_{i=1}^3 c_i\\mathbf{s}_i e^{-\\lambda_i t}$$\n",
    "\n",
    "which can be solved given the eigenvalues and eigenvectors of $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Examples:  Particle tracking in a fluid\n",
    "\n",
    "$$\\frac{\\text{d} \\mathbf{X}}{\\text{d}t} = \\mathbf{V}(t, \\mathbf{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": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0., 1., 11)\n",
    "y = numpy.linspace(0., 1., 11)\n",
    "x_fine = numpy.linspace(0., 1.)\n",
    "y_fine = numpy.linspace(0., 1.)\n",
    "\n",
    "X, Y = numpy.meshgrid(x,y)\n",
    "X_fine, Y_fine = numpy.meshgrid(x_fine, y_fine)\n",
    "\n",
    "pi = numpy.pi\n",
    "psi = numpy.sin(pi*X_fine)*numpy.sin(pi*Y_fine)\n",
    "U = pi*numpy.sin(pi*X)*numpy.cos(pi*Y)\n",
    "V = -pi*numpy.cos(pi*X)*numpy.sin(pi*Y)\n",
    "\n",
    "x0 = 0.75\n",
    "y0 = 0.75\n",
    "psi0 = numpy.sin(pi*x0)*numpy.sin(pi*y0)\n",
    "\n",
    "fig = plt.figure(figsize=(8,8))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.quiver(X,Y, U, V)\n",
    "axes.plot(.75, 0.75,'ro')\n",
    "axes.contour(X_fine, Y_fine, psi, [ psi0 ])\n",
    "axes.grid()\n",
    "axes.set_title(\"Particle tracking\", fontsize=18)\n",
    "axes.set_xlabel('y', fontsize=16)\n",
    "axes.set_ylabel('x', fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "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",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\\mathbf{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} = \\mathbf{f}(t, \\mathbf{u})$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "def f_vanderpol(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",
    "t_span = (0., 50.)\n",
    "u0 = [ 1., 1. ]\n",
    "f = lambda t, u: f_vanderpol(t, u, mu=5)\n",
    "sol = solve_ivp(f, t_span, u0,method='BDF',rtol=1.e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "\n",
    "axes.plot(sol.t, sol.y[0])\n",
    "axes.set_title(\"Solution to Van der Pol Oscillator\", fontsize=18)\n",
    "axes.set_xlabel(\"t\", fontsize=16)\n",
    "axes.set_ylabel(\"y(t)\", fontsize=16)\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "\n",
    "axes.plot(sol.y[0],sol.y[1],'r')\n",
    "axes.set_title(\"Phase Diagram for Van der Pol Oscillator\", fontsize=18)\n",
    "axes.set_xlabel(\"y(t)\", fontsize=16)\n",
    "axes.set_ylabel(\"y'(t)\", fontsize=16)\n",
    "axes.axis('equal')\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Big Questions\n",
    "\n",
    "Given a RHS $\\mathbf{f}(t,\\mathbf{u})$ and initial condition $\\mathbf{u}(t_0) \\ldots$\n",
    "\n",
    "*  How do you find a discrete numerical solution that approximates the trajectory $\\mathbf{u}(t)$?\n",
    "*  How do you control the accuracy of the approximation?\n",
    "*  How do you improve the efficiency of the approximation?\n",
    "*  How do you understand  stability and convergence?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Some Notation:  Basic Stepping schemes\n",
    "\n",
    "Introducing some notation to simplify 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",
    "$$\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": "fragment"
    }
   },
   "source": [
    "One obvious way to do this, is to just use $f(t, u(t))$ and write the update scheme as\n",
    "\n",
    "$$\n",
    "    u(t + \\Delta t) = u(t) + \\Delta t f(t,u(t))\n",
    "$$\n",
    "\n",
    "Which is our first integration scheme (which goes by the name Euler's method).  As usual, though the first scheme is often the worst scheme, but with a bit of understanding we can do much better with not a lot more work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Integral form of ODE IVP's: the  Relationship to quadrature\n",
    "\n",
    "Euler's method is an example of a \"Single Step, multi-stage\" scheme of which there are many.  However, to derive them it is actually more instructive to work with the integral form of an ODE, which will put the equations in a form where we can use our ideas from quadrature to make progress"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Given a a system of ODE's\n",
    "$$\n",
    "    \\frac{d\\mathbf{u}}{dt} = \\mathbf{f}(t,\\mathbf{u})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can integrate both sides\n",
    "$$\n",
    "    \\int^{t + \\Delta t}_t \\frac{d\\mathbf{u}}{d \\tau} d\\tau = \\int^{t + \\Delta t}_t \\mathbf{f}(\\tau, \\mathbf{u}) d\\tau\n",
    "$$\n",
    "\n",
    "which is equivalent to the differential form. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    " However using the fundamental theorem of calculus tells us that the LHS is $u(t + \\Delta t) - u(t)$  so we can write the ODE as\n",
    "\n",
    "$$ \n",
    "   \\mathbf{u}(t + \\Delta t) = \\mathbf{u}(t) + \\int^{t + \\Delta t}_t \\mathbf{f}(\\tau, \\mathbf{u}(\\tau)) d\\tau\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Single-Step Multi-Stage Schemes\n",
    "\n",
    "The integral form of an ODE initial value problem can be written\n",
    "\n",
    "$$ \n",
    "    u(t + \\Delta t) = u(t) + \\int^{t + \\Delta t}_t f(\\tau, u(\\tau)) d\\tau\n",
    "$$\n",
    "\n",
    "Which says that our solution $u$, if it exists  at some time $\\Delta t$ in the future,  is $u(t)$ plus a *number* \n",
    "\n",
    "$$\n",
    "    K = \\int^{t + \\Delta t}_t f(\\tau, u(\\tau) )d\\tau\n",
    "$$\n",
    "\n",
    "which is a definite *line integral* (along an unknown solution). \n",
    "\n",
    "An important class of ODE solvers are called *Single Step, Multi-stage schemes* which can be most easily understood as extensions of the  Newton-Cotes quadrature schemes for approximating $K$ (plus an error term that will scale as $\\Delta t^p$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The Geometric picture\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = numpy.linspace(0., 4800, 11)\n",
    "y = numpy.linspace(0., 1.2, 11)\n",
    "T, Y = numpy.meshgrid(t,y)\n",
    "dt = numpy.ones(T.shape)\n",
    "dy = -dt*Y\n",
    "\n",
    "tK = 2000.\n",
    "uK = numpy.exp(decay_constant*tK)\n",
    "K = uK -1.\n",
    "tp = numpy.linspace(0., 4800, 100)\n",
    "tk = numpy.linspace(0., tK, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(10,8))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.quiver(T,Y, dt,dy, color='gray')\n",
    "axes.plot(tp,numpy.exp(decay_constant*tp))\n",
    "axes.plot(0.,1.,'ro')\n",
    "axes.plot(tk,numpy.exp(decay_constant*tk),'r--')\n",
    "axes.plot(tK, uK, 'ro')\n",
    "axes.plot([0.,0.], [1., uK], 'r--')\n",
    "axes.text(10., 0.72, '$K$', fontsize=24, color='red')\n",
    "axes.plot([0.,tK],[uK, uK], 'r--')\n",
    "axes.text(900., uK - .1, '$\\Delta t$', fontsize=24, color='red')\n",
    "axes.text(-10, 1., '$U_0$', fontsize=24, color='blue')\n",
    "axes.text(tK+10, uK, '$U_1$', fontsize=24, color='blue')\n",
    "\n",
    "\n",
    "axes.grid()\n",
    "axes.set_title(\"Direction Set, $u' = - \\lambda u$, $u(0)=1$\", fontsize=18)\n",
    "axes.set_xlabel('t (years)', fontsize=16)\n",
    "axes.set_ylabel('u', fontsize=16)\n",
    "\n",
    "axes.set_ylim((-.1,1.2))\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "#### Exercise (left for the reader)\n",
    "\n",
    "Show that \n",
    "$$\n",
    "    u(\\Delta t) = u_0 + \\int_0^{\\Delta t} f(\\tau, u(\\tau)) d\\tau\n",
    "$$\n",
    "\n",
    "when $f(\\tau, u(\\tau)) = -\\lambda u(\\tau)$ and $u(\\tau) = u_0e^{-\\lambda \\tau}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Forward Euler scheme\n",
    "For example, if we approximate $K$ with a left-sided quadrature rule\n",
    "$$\n",
    "    K = \\int^{t + \\Delta t}_t f(\\tau, u(\\tau)) d\\tau \\approx \\Delta t f(t, u(t))\n",
    "$$\n",
    "\n",
    "then our first ODE algorithm can be written\n",
    "$$\n",
    "u(t + \\Delta t) =  u(t) + \\Delta t f(t, u(t))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which is exactly Euler's method that we derived previously"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "in terms of our discrete approximation $U$\n",
    "$$\n",
    "\\begin{align}\n",
    "    K_1 &= \\Delta t f(t_n, U_n)\\\\\n",
    "    U_{n+1} &= U_n + K_1\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    " 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",
    "# Euler step\n",
    "dt = 1e3\n",
    "u_np = c_0 + dt * (decay_constant * c_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, c_0 * numpy.exp(decay_constant * t), label=\"True Solution\")\n",
    "axes.plot(0., 1., 'ro')\n",
    "axes.text(0., 1.01, '$U_0$', fontsize=16)\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.text(10., 0.75, '$K_1$', fontsize=16)\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.05))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Implement Forward Euler\n",
    "def euler(f, t_span, u0, N):\n",
    "    \"\"\" simple implementation of constant step-size forward euler method\n",
    "        This doc string should have so much more in it\n",
    "    \"\"\"\n",
    "    t = numpy.linspace(t_span[0], t_span[1],N)\n",
    "    u = numpy.empty(t.shape)\n",
    "    u[0] = u0\n",
    "    delta_t = t[1] - t[0]\n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        K1 = delta_t * f(t_n, u[n])\n",
    "        u[n + 1] = u[n] + K1        \n",
    "    return t, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "t_span = [0.0, 1.6e3]\n",
    "u0 = 1.\n",
    "N = 21\n",
    "t_euler, u_euler = euler(f, t_span, u0, N)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "u_exact = lambda t : c_0 * numpy.exp(decay_constant * t)\n",
    "\n",
    "fig = plt.figure(figsize=(16, 6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t_euler, u_euler, 'or', label=\"Euler\")\n",
    "axes.plot(t_exact, u_exact(t_exact), 'k--', label=\"True Solution\")\n",
    "\n",
    "axes.set_title(\"Forward Euler\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.grid()\n",
    "axes.legend()\n",
    "\n",
    "abs_err = numpy.abs(u_euler - u_exact(t_euler))\n",
    "rel_err = abs_err/u_exact(t_euler)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(t_euler,abs_err,'ro',label='absolute error')\n",
    "axes.plot(t_euler,rel_err,'bo',label='relative error')\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"error\")\n",
    "axes.set_title('Error')\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Backward's Euler\n",
    "\n",
    "Similar to forward Euler is the *backward Euler* method which,   uses a right-rectangle rule to estimate $K$ given $f$ at a future time. i.e. \n",
    "\n",
    "$$\n",
    "    K\\approx \\Delta t f(t_{n+1}, U_{n+1})\n",
    "$$\n",
    "\n",
    "However, the update scheme now becomes\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t f(t_{n+1}, U_{n+1}).\n",
    "$$\n",
    "\n",
    "which requires a (usually non-linear) solve for $U_{n+1}$. 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 by hand.  For instance in the case of our example problem, $f=\\lambda U$, we have:\n",
    "\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",
    "\n",
    "which can be solved for $U_{n+1}$ to find\n",
    "\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": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "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(figsize=(8,6))\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 Backwards 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, u_np, 'ro')\n",
    "axes.text(dt+ 10., u_np, '$U_1$', fontsize=16)\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",
    "axes.text(10., 0.85, '$K_1$', fontsize=16)\n",
    "\n",
    "axes.grid()\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.05))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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_steps = 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t_exact = numpy.linspace(0.0, 1.6e3, 100)\n",
    "\n",
    "# Implement backwards Euler\n",
    "t_backwards = numpy.linspace(0.0, 1.6e3, n_steps)\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(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t_backwards, u_backwards, 'or', label=\"Backwards Euler\")\n",
    "axes.plot(t_exact, u_exact(t_exact), 'k--', label=\"True Solution\")\n",
    "axes.grid()\n",
    "axes.set_title(\"Backwards Euler\")\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"$c(t)$\")\n",
    "axes.set_ylim((0.4,1.1))\n",
    "axes.legend()\n",
    "\n",
    "abs_err = numpy.abs(u_backwards - u_exact(t_backwards))\n",
    "rel_err = abs_err/u_exact(t_backwards)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(t_backwards,abs_err,'ro',label='absolute error')\n",
    "axes.plot(t_backwards,rel_err,'bo',label='relative error')\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"error\")\n",
    "axes.set_title('Error')\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "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",
    "     \\left [ I - \\Delta t A \\right ]U_{n+1} &= 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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Another simple implicit method is based on quadrature 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 for $f(u) = \\lambda u$?"
   ]
  },
  {
   "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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "n_steps = 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "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",
    "\n",
    "# Implement trapezoidal method\n",
    "t = numpy.linspace(0.0, 1.6e3, n_steps)\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(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t, u, 'or', label=\"Trapezoidal\")\n",
    "axes.plot(t_exact, u_exact(t_exact), 'k--', label=\"True Solution\")\n",
    "axes.grid()\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",
    "\n",
    "abs_err = numpy.abs(u  - u_exact(t))\n",
    "rel_err = abs_err/u_exact(t)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(t,abs_err,'ro',label='absolute error')\n",
    "axes.plot(t,rel_err,'bo',label='relative error')\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"error\")\n",
    "axes.set_title('Error')\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "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,U_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, \\tau]$.  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",
    "$$\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": "fragment"
    }
   },
   "source": [
    "Evaluating this series at $t_{n+1}$ gives\n",
    "\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": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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",
    "$$\n",
    "    \\frac{U_{n+1} - U_n}{\\Delta t} = f(t_n, U_n)\n",
    "$$\n",
    "\n",
    "and replacing the derivative formulation with $u(t_n)$ to find\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    T(t, u; \\Delta t) &= \\frac{u(t_{n+1}) - u(t_n)}{\\Delta t} - f(t_n, u_n) \\\\\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Given the Taylor's series expansion for $u(t_{n+1})$\n",
    "\n",
    "$$\n",
    "u(t_{n+1}) =u(t_n) + \\Delta t f(t_n, u_n) + \\frac{u''(t_n)}{2} \\Delta t^2 + \\mathcal{O}(\\Delta t^3)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We substitute to find \n",
    "$$\n",
    "    T(t, u; \\Delta t)  = \\frac{u''(t_n)}{2} \\Delta t + \\mathcal{O}(\\Delta t^2).\n",
    "$$\n",
    "\n",
    "This implies that forward Euler is first order accurate and therefore consistent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "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": [
    "#### Truncation Error vs Step Error\n",
    "\n",
    "Sometimes we will also consider the \"Step Error\"  which is the error that is introduced over one step\n",
    "\n",
    "$$\n",
    "    E_h = | U_{n+1} - u_{n+1} |\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This leads to an alternate definition of the truncation error as\n",
    "\n",
    "$$\n",
    "    T(t,u;\\Delta t) = \\frac{E_h}{\\Delta t} = \\frac{1}{\\Delta t} [U_{n+1} - u_{n+1}]\n",
    "$$\n",
    "\n",
    "so if the Truncation error is $O(\\Delta t^p)$ then the step error will be order $O(\\Delta t^{p+1})$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So for Forward (or Backward's) Euler the step error \n",
    "\n",
    "$$\n",
    "    E_h = O(\\Delta t^2)\n",
    "$$\n",
    "\n",
    "The step error can be very useful in *adaptive stepping* schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Runge-Kutta Methods\n",
    "\n",
    "One way to derive higher-order ODE solvers is to use higher order quadrature schemes that sample the function at a number of  intermediate stages to provide a more accurate estimate of $K$.  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 $K$ using a mid-point scheme (which should be 2nd order accurate).  Unforrunately, we don't know the value of the mid-point.  However we can use an Euler step of size $\\Delta t/2$ to estimate the mid-point.  \n",
    "\n",
    "We can write the algorithm as \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    K_1 &= \\Delta t f(U_n, t_n) \\\\\n",
    "    K_2 &= \\Delta t f(U_n + K_1/2, t_n + \\Delta t/2 )\\\\\n",
    "    U_{n+1} &= U_n + K_2 \\\\    \n",
    "\\end{aligned}$$\n",
    "\n",
    "Where we now evaluate the function in two stages $K_1$ and $K_2$.\n",
    "\n",
    "or for an autonomous ODE\n",
    "$$\n",
    "    U_{n+1} = U_n + \\Delta t f(U_n + \\frac{1}{2} \\Delta t f(U_n))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "# RK2 step\n",
    "dt = 1e3\n",
    "U0 = 1.0\n",
    "K1 = dt * f(0., U0)\n",
    "Y1 = U0 + K1/2\n",
    "K2 = dt * f(dt/2., Y1)\n",
    "U1 = U0 + K2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U0 * numpy.exp(decay_constant * t), label=\"True Solution\")\n",
    "axes.plot(0., U0, 'ro')\n",
    "axes.text(0., U0+.01, '$U_0$', fontsize=16)\n",
    "axes.plot((0.0, dt), (U0, U0 + K1), 'k--')\n",
    "axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k')\n",
    "\n",
    "axes.plot(dt/2., U0 + K1/2, 'ro')\n",
    "axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "axes.text(10., 0.85, '$\\\\frac{K_1}{2}$', fontsize=18)\n",
    "axes.plot((0.0, dt/2), (Y1, Y1), 'k--')\n",
    "axes.text(250, Y1 - 0.05, '$\\\\frac{\\Delta t}{2}$', fontsize=18)\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "axes.text(10., 0.85, '$\\\\frac{K_1}{2}$', fontsize=18)\n",
    "axes.plot((0.0, dt/2), (Y1, Y1), 'k--')\n",
    "axes.text(250, Y1 - 0.05, '$\\\\frac{\\Delta t}{2}$', fontsize=18)\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "axes.plot((0., dt), (U0, U1), 'k')\n",
    "axes.text(dt+20, U1, '$U_1$', fontsize=18)\n",
    "#axes.plot((0.0, 0.0), (U0, U1), 'g--')\n",
    "#axes.plot((0.0, dt), (U1, U1), 'g--')\n",
    "\n",
    "\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.05))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Error analysis RK2\n",
    "\n",
    "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 with respect to $t$ 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": [
    "Using our alternative definition of the truncation error we have\n",
    "\n",
    "$$\n",
    "        T(t, u; \\Delta t) = \\frac{1}{\\Delta t} \\left[U_{n+1} - u_{n+1} \\right] \n",
    "$$\n",
    "or"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\\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:  Improved Euler's method\n",
    "\n",
    "The Improved Euler's method is another RK2 scheme but instead of approximating a mid-point quadrature rule, it approximates a trapezoidal rule.   \n",
    "\n",
    "We can write the algorithm as \n",
    "\n",
    "$$\\begin{aligned}\n",
    "    K_1 &= \\Delta t f(U_n, t_n) \\\\\n",
    "    K_2 &= \\Delta t f(U_n + K1, t_n + \\Delta t )\\\\\n",
    "    U_{n+1} &= U_n + \\frac{1}{2}\\left[K_1 +K_2\\right] \\\\    \n",
    "\\end{aligned}$$\n",
    "\n",
    "Where we now use function evaluations at both the initial value, and at the euler point but take the average of those slopes.\n",
    "\n",
    "Again, error analysis shows that this scheme also has a truncation error $T(t,u:\\Delta t) = \\mathcal{O}(\\Delta t^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "# Improved Euler step\n",
    "dt = 1e3\n",
    "U0 = 1.0\n",
    "K1 = dt * f(0., U0)\n",
    "Y1 = U0 + K1\n",
    "K2 = dt * f(dt, Y1)\n",
    "U1 = U0 + 0.5*(K1 + K2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U0 * numpy.exp(decay_constant * t), label=\"True Solution\")\n",
    "axes.plot(0., U0, 'ro')\n",
    "axes.text(0., U0+.01, '$U_0$', fontsize=16)\n",
    "axes.plot((0.0, dt), (U0, U0 + K1), 'k--')\n",
    "axes.plot((0.0, dt), (U0, U0 + K1), 'k')\n",
    "\n",
    "axes.plot(dt, Y1, 'ro')\n",
    "axes.text(dt+10, Y1, '$Y_1$', fontsize=18)\n",
    "axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "axes.plot((0.0, dt), (Y1, Y1), 'k--')\n",
    "axes.text(350, Y1 - 0.05, '$\\\\frac{\\Delta t}{2}$', fontsize=18)\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "axes.plot(0., Y1, 'ko')\n",
    "axes.text(10., Y1, '$K_1$', fontsize=18)\n",
    "\n",
    "axes.plot((0., 0.), (U0, U0+K2),'b--')\n",
    "axes.plot(0., U0+K2,'bo--')\n",
    "axes.text(10., U0+K2, '$K_2$', fontsize=18)\n",
    "\n",
    "\n",
    "axes.plot(0., U1,'gx', markersize=15)\n",
    "axes.text(10., U1, '$0.5*(K_1 +K_2)$', fontsize=18)\n",
    "\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "axes.plot((0., dt), (U0, U1), 'k')\n",
    "axes.text(dt+20, U1, '$U_1$', fontsize=18)\n",
    "#axes.plot((0.0, 0.0), (U0, U1), 'g--')\n",
    "#axes.plot((0.0, dt), (U1, U1), 'g--')\n",
    "\n",
    "\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.05))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  4-stage Runge-Kutta Method\n",
    "\n",
    "If RK2 is related to a Mid-point quadrature scheme,  then the classic 4-stage, 4th order Runge-Kutta scheme should be reminiscent of Simpson's Quadrature rule.  It requires 4 samples of $f(t,u)$ at the beginning of the step, two-samples in the middle and one at the end, then a linear combination of those samples\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    K_1 &= \\Delta t f(t_n, U_n) \\\\\n",
    "    K_2 &= \\Delta t f(t_n + \\Delta t/2, U_n + K_1/2) \\\\\n",
    "    K_3 &= \\Delta t f(t_n + \\Delta t/2, U_n + K_2/2) \\\\\n",
    "    K_4 &= \\Delta t f(t_n + \\Delta t, U_n + K_3) \\\\\n",
    "        & \\\\\n",
    "    U_{n+1} &= U_n + \\frac{1}{6} \\left [K_1 + 2(K_2 + K_3)  + K_4) \\right ] \n",
    "\\end{aligned}$$\n",
    "\n",
    "With truncation error $T = O(\\Delta t^4)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "# RK4 step\n",
    "dt = 1e3\n",
    "U0 = 1.0\n",
    "K1 = dt * f(0., U0)\n",
    "K2 = dt * f(dt/2., U0 + K1/2)\n",
    "K3 = dt * f(dt/2., U0 + K2/2)\n",
    "K4 = dt * f(dt, U0 + K3)\n",
    "\n",
    "U1 = U0 + 1./6. *( K1 + 2 * (K2 + K3) + K4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, U0 * numpy.exp(decay_constant * t), label=\"True Solution\")\n",
    "axes.plot(0., U0, 'ro')\n",
    "axes.text(0.-20, U0-.04, '$K_1$', color='red',fontsize=16)\n",
    "\n",
    "axes.text(0., U0+.01, '$U_0$', fontsize=16)\n",
    "axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k--')\n",
    "axes.plot(dt/2., U0 + K1/2, 'ro')\n",
    "axes.text(dt/2-20, U0 + K1/2-.04, '$K_2$', color='red',fontsize=16)\n",
    "\n",
    "\n",
    "axes.plot((0.0, dt/2.), (U0, U0 + K2/2.), 'k--')\n",
    "axes.plot(dt/2., U0 + K2/2, 'ro')\n",
    "axes.text(dt/2-20, U0 + K2/2+.02, '$K_3$', color='red',fontsize=16)\n",
    "\n",
    "axes.plot((0.0, dt), (U0, U0 + K3), 'k--')\n",
    "axes.plot(dt, U0 + K3, 'ro')\n",
    "axes.text(dt-20, U0 + K3-.04, '$K_4$', color='red',fontsize=16)\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "#axes.plot((0., dt), (U0, U1), 'k')\n",
    "axes.text(dt+20, U1, '$U_1$', fontsize=18)\n",
    "\n",
    "\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.05))\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def RK2(f, t_span, u0, N):\n",
    "    \"\"\" implement constant step size 2 stage Runge-Kutta Method RK2\"\"\"\n",
    "    \n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    delta_t = t[1] - t[0]\n",
    "    u = numpy.empty(t.shape)\n",
    "    u[0] = u0 \n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        K_1 = delta_t * f(t_n, u[n])\n",
    "        K_2 = delta_t * f(t_n + delta_t/2., u[n] + K_1/2.)\n",
    "        u[n+1] = u[n] + K_2\n",
    "    return t, u\n",
    "\n",
    "def improved_euler(f, t_span, u0, N):\n",
    "    \"\"\" implement constant step size 2 stage Improved Euler Method trapezoidal rule\"\"\"\n",
    "    \n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    delta_t = t[1] - t[0]\n",
    "    u = numpy.empty(t.shape)\n",
    "    u[0] = u0\n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        K_1 = delta_t * f(t_n, u[n])\n",
    "        K_2 = delta_t * f(t_n + delta_t, u[n] + K_1)\n",
    "        u[n+1] = u[n] + 0.5 * (K_1 + K_2)\n",
    "    return t, u "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def RK4(f, t_span, u0, N):\n",
    "    \"\"\" implement constant step size 4 stage Runge-Kutta Method RK4\"\"\"\n",
    "    \n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    delta_t = t[1] - t[0]\n",
    "    u = numpy.empty(t.shape)\n",
    "    u[0] = u0\n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        K_1 = delta_t * f(t_n, u[n])\n",
    "        K_2 = delta_t * f(t_n + delta_t/2., u[n] + K_1/2.)\n",
    "        K_3 = delta_t * f(t_n + delta_t/2., u[n] + K_2/2.)\n",
    "        K_4 = delta_t * f(t_n + delta_t, u[n] + K_3)\n",
    "        u[n+1] = u[n] + 1./6. * (K_1 + 2.*( K_2 + K_3) + K_4)\n",
    "    return t, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Implement and compare the two-stage and 4-stage Runge-Kutta methods \n",
    "f = lambda t, u: -u\n",
    "N = 21\n",
    "t_span = [ 0., 5.0 ]\n",
    "u0 = 1.\n",
    "\n",
    "u_exact = lambda t: u0*numpy.exp(-t)\n",
    "\n",
    "t_exact = numpy.linspace(t_span[0], t_span[1], 100)\n",
    "t_euler, u_euler = euler(f, t_span, u0, N)\n",
    "t_ieuler, u_ieuler = improved_euler(f, t_span, u0, N)\n",
    "t_RK2, u_RK2 = RK2(f, t_span, u0, N)\n",
    "t_RK4, u_RK4 = RK4(f, t_span, u0, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(t_exact,u_exact(t_exact),'k',label='exact')\n",
    "axes.plot(t_euler, u_euler, 'ro', label='euler')\n",
    "axes.plot(t_ieuler, u_ieuler, 'co', label='improved euler')\n",
    "axes.plot(t_RK2, u_RK2, 'go', label='RK2')\n",
    "axes.plot(t_RK4, u_RK4, 'bo', label='RK4')\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('t', fontsize=16)\n",
    "axes.set_ylabel('u', fontsize=16)\n",
    "axes.legend(loc='best')\n",
    "\n",
    "err = lambda u, t: numpy.abs(u - u_exact(t))/u_exact(t)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "\n",
    "axes.semilogy(t_euler,err(u_euler,t_euler),'ro',label='euler')\n",
    "axes.semilogy(t_ieuler,err(u_ieuler,t_ieuler),'co',label='improved euler')\n",
    "axes.semilogy(t_RK2,err(u_RK2,t_RK2),'go',label='RK2')\n",
    "axes.semilogy(t_RK4,err(u_RK4,t_RK4),'bo',label='RK4')\n",
    "\n",
    "axes.set_xlabel(\"t (years)\")\n",
    "axes.set_ylabel(\"Rel. error\")\n",
    "axes.set_title('Error')\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence of Single Step Multi-Stage schemes\n",
    "\n",
    "All of the above schemes are consistent and have truncation errors $T\\propto\\Delta t^p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = numpy.array([ 2**n for n in range(4,10)])\n",
    "err_euler = numpy.zeros(len(N))\n",
    "err_ieuler = numpy.zeros(len(N))\n",
    "err_RK2 = numpy.zeros(len(N))\n",
    "err_RK4 = numpy.zeros(len(N))\n",
    "\n",
    "t_span = [ 0., 4.]\n",
    "dt = t_span[1]/N\n",
    "\n",
    "u0 = 1. \n",
    "u_exact = u0*numpy.exp(-t_span[1])\n",
    "\n",
    "for i, n in enumerate(N):\n",
    "    t, u_euler = euler(f, t_span, u0, n)\n",
    "    err_euler[i] = numpy.abs(u_euler[-1] - u_exact)\n",
    "    t, u_ieuler = improved_euler(f, t_span, u0, n)\n",
    "    err_ieuler[i] = numpy.abs(u_ieuler[-1] - u_exact)\n",
    "    t, u_RK2 = RK2(f, t_span, u0, n)\n",
    "    err_RK2[i] = numpy.abs(u_RK2[-1] - u_exact)\n",
    "    t, u_RK4 = RK4(f, t_span, u0, n)\n",
    "    err_RK4[i] = numpy.abs(u_RK4[-1] - u_exact)\n",
    "    \n",
    "err_fit = lambda dt, p: numpy.exp(p[1])*dt**p[0]\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# Euler\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_euler[2:]),1)\n",
    "line = axes.loglog(dt, err_euler, 'o', label='euler, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "# Improved Euler\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_ieuler[2:]),1)\n",
    "line = axes.loglog(dt, err_ieuler, 'o', label='improved euler, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "# RK2\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK2[2:]),1)\n",
    "line = axes.loglog(dt, err_RK2, 'o', label='rk2, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "#RK4\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK4[2:]),1)\n",
    "line = axes.loglog(dt, err_RK4, 'o', label='rk4, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('$\\Delta t$', fontsize=16)\n",
    "axes.set_ylabel('$Error$', fontsize=16)\n",
    "axes.set_title('Convergence: Single Step Schemes', fontsize=18)\n",
    "axes.legend(loc='best', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary: single-Step Multi-Stage Schemes\n",
    "\n",
    "The integral form of an ODE initial value problem can be written\n",
    "\n",
    "$$ \n",
    "    u(t + \\Delta t) = u(t) + \\int^{t + \\Delta t}_t f(\\tau, u(\\tau)) d\\tau\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which says that our solution $u$, if it exists  at some time $\\Delta t$ in the future,  is $u(t)$ plus a *number* \n",
    "\n",
    "$$\n",
    "    K = \\int^{t + \\Delta t}_t f(\\tau, u(\\tau) )d\\tau\n",
    "$$\n",
    "\n",
    "which is a definite *line integral* (along an unknown solution). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Single Step, Multi-stage schemes\n",
    "\n",
    "are most easily understood as extensions of the  Newton-Cotes quadrature schemes for approximating $K$ (plus an error term that will scale as $\\Delta t^p$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Explicit Schemes**\n",
    "<table width=\"80%\">\n",
    "    <tr align=\"center\"><th>Name</th> <th align=\"center\">Stages</th> <th align=\"center\">\"Quadrature\"</th><th align=\"center\">$$T$$</th></tr>\n",
    "     <tr align=\"center\"><td>Euler</td> <td align=\"center\">1</td> <td align=\"center\">Left-Rectangle</td><td align=\"center\">$$O(\\Delta t)$$</td></tr>\n",
    "    <tr align=\"center\"><td>Improved Euler</td> <td align=\"center\">2</td> <td align=\"center\">Trapezoidal</td><td align=\"center\">$$O(\\Delta t^2)$$</td></tr>\n",
    "    <tr align=\"center\"><td>RK2</td> <td align=\"center\">2</td> <td align=\"center\">Mid-Point</td><td align=\"center\">$$O(\\Delta t^2)$$</td></tr>\n",
    "    <tr align=\"center\"><td>RK4</td> <td align=\"center\">4</td> <td align=\"center\">Simpson</td><td align=\"center\">$$O(\\Delta t^4)$$</td></tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Implicit Schemes**\n",
    "<table width=\"80%\">\n",
    "    <tr align=\"center\"><th>Name</th> <th align=\"center\">Stages</th> <th align=\"center\">\"Quadrature\"</th><th align=\"center\">$$T$$</th></tr>\n",
    "     <tr align=\"center\"><td>Backwards-Euler</td> <td align=\"center\">1</td> <td align=\"center\">Right-Rectangle</td><td align=\"center\">$$O(\\Delta t)$$</td></tr>\n",
    "    <tr align=\"center\"><td>Trapezoidal</td> <td align=\"center\">2</td> <td align=\"center\">Trapezoidal</td><td align=\"center\">$$O(\\Delta t^2)$$</td></tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Convergence of Single Step Multi-Stage schemes\n",
    "\n",
    "All of the above schemes are consistent and have truncation errors $T\\propto\\Delta t^p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "N = numpy.array([ 2**n for n in range(4,10)])\n",
    "err_euler = numpy.zeros(len(N))\n",
    "err_ieuler = numpy.zeros(len(N))\n",
    "err_RK2 = numpy.zeros(len(N))\n",
    "err_RK4 = numpy.zeros(len(N))\n",
    "\n",
    "t_span = [ 0., 4.]\n",
    "dt = t_span[1]/N\n",
    "\n",
    "u0 = 1. \n",
    "u_exact = u0*numpy.exp(-t_span[1])\n",
    "\n",
    "for i, n in enumerate(N):\n",
    "    t, u_euler = euler(f, t_span, u0, n)\n",
    "    err_euler[i] = numpy.abs(u_euler[-1] - u_exact)\n",
    "    t, u_ieuler = improved_euler(f, t_span, u0, n)\n",
    "    err_ieuler[i] = numpy.abs(u_ieuler[-1] - u_exact)\n",
    "    t, u_RK2 = RK2(f, t_span, u0, n)\n",
    "    err_RK2[i] = numpy.abs(u_RK2[-1] - u_exact)\n",
    "    t, u_RK4 = RK4(f, t_span, u0, n)\n",
    "    err_RK4[i] = numpy.abs(u_RK4[-1] - u_exact)\n",
    "    \n",
    "err_fit = lambda dt, p: numpy.exp(p[1])*dt**p[0]\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# Euler\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_euler[2:]),1)\n",
    "line = axes.loglog(dt, err_euler, 'o', label='euler, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "# Improved Euler\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_ieuler[2:]),1)\n",
    "line = axes.loglog(dt, err_ieuler, 'o', label='improved euler, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "# RK2\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK2[2:]),1)\n",
    "line = axes.loglog(dt, err_RK2, 'o', label='rk2, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "#RK4\n",
    "p = numpy.polyfit(numpy.log(dt[2:]), numpy.log(err_RK4[2:]),1)\n",
    "line = axes.loglog(dt, err_RK4, 'o', label='rk4, p={:3.2f}'.format(p[0]))\n",
    "axes.loglog(dt, err_fit(dt,p),'--', color=line[0].get_color())\n",
    "\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('$\\Delta t$', fontsize=16)\n",
    "axes.set_ylabel('$Error$', fontsize=16)\n",
    "axes.set_title('Convergence: Single Step Schemes', fontsize=18)\n",
    "axes.legend(loc='best', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## General Explicit Runge-Kutta Schemes and ``Butcher Tableaus''\n",
    "\n",
    "\n",
    "The most general form of an explicit RK scheme with $s$ stages is\n",
    "\n",
    "$$\n",
    "    u_{n+1} = u_n +  \\sum_{i=1}^s b_i K_{i}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "K_{i} = \\Delta t f\\left(t_n + c_i \\Delta t, u_n + \\sum_{j=1}^{i-1} a_{ij} K_{j}\\right)\n",
    "$$\n",
    "\n",
    "$K_{i}$ are the function evaluations at stage $i$ which are determined by coefficients $a_{ij}$, $b_i$ and $c_i$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All of the critical coefficients for any explicit RK scheme can be arranged in a ``Butcher tableau'' \n",
    "\n",
    "$$\\begin{array} {c|ccccc} 0 & 0\\\\ c_2 & a_{21}\\\\ c_3 & a_{31} & a_{32} \\\\ \\vdots & \\vdots & & \\ddots \\\\ c_s& a_{s1}& a_{s2} & \\cdots & a_{s,s-1}\\\\ \\hline & b_1 & b_2 & \\cdots & b_{s-1} & b_s \\end{array}\n",
    "$$\n",
    "\n",
    "where $c_i$ gives the fraction of the time step for each stage.  $b_i$ gives the weighting for each function evaluation $K_i$ in the general quadrature and the $a_{ij}$ determine how to use previous stages to advance the value of $u$ at stage $i$.\n",
    "\n",
    "$$\n",
    "K_{i} = \\Delta t f\\left(t_n + c_i \\Delta t, u_n + \\sum_{j=1}^{i-1} a_{ij} K_{j}\\right)\n",
    "$$\n",
    "\n",
    "for more information see [here](https://www.johndcook.com/blog/2020/02/13/runge-kutta-methods/#:~:text=For%20explicit%20Runge%2DKutta%20methods,Butcher%20who%20classified%20RK%20methods)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "All of the schemes we have discussed can be written as Butcher Tableaus (although they can be a bit hard to parse).  But for example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Euler's method\n",
    "\n",
    "$$\n",
    "\\begin{array} {c|c} 0 & 0\\\\ \\hline & 1\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Backwards Euler's method\n",
    "\n",
    "$$\n",
    "\\begin{array} {c|c} 1 & 1\\\\ \\hline & 1\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Classical RK4\n",
    "\n",
    "$$\n",
    "\\begin{array} {c|cccc} 0\\\\ 1/2 & 1/2\\\\ 1/2 &0 &1/2 \\\\ 1& 0& 0& 1\\\\ \\hline & 1/6 & 1/3 & 1/3 &1/6 \\end{array}\n",
    "$$\n",
    "\n",
    "Wikipedia provides a convenient [list](https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods)  of a large range of explicit and implicit RK schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Adaptive Time Stepping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Why should we care about all of these schemes and their errors?\n",
    "\n",
    "* Even though we know the formal error.  It is with respect to a true solution we don't know.\n",
    "* In itself, the error estimates don't tell us how to choose a time step $\\Delta t$ to keep the solution within a given tolerance\n",
    "* However,  in combination, we can use multiple methods to control the error and provide **Adaptive** time stepping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:  Compare 1 step of Euler to one step of RK2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "# RK2 step\n",
    "dt = 1e3\n",
    "U0 = 1.0\n",
    "K1 = dt * f(0., U0)\n",
    "Y1 = U0 + K1/2\n",
    "K2 = dt * f(dt/2., Y1)\n",
    "U1 = U0 + K2\n",
    "\n",
    "t = numpy.linspace(U0, 1600.)\n",
    "u_true = U0 * numpy.exp(decay_constant * t)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t, u_true, label=\"True Solution\")\n",
    "axes.plot(0., U0, 'ro')\n",
    "axes.text(0., U0+.01, '$U_0$', fontsize=16)\n",
    "axes.plot((0.0, dt), (U0, U0 + K1), 'k')\n",
    "#axes.plot((0.0, dt/2.), (U0, U0 + K1/2.), 'k')\n",
    "\n",
    "#Euler step\n",
    "axes.plot(dt, U0 + K1, 'ro')\n",
    "#axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "axes.text(dt + 10., U0 + K1, '$U_{euler}$', fontsize=18)\n",
    "\n",
    "axes.plot(dt, U1, 'go')\n",
    "#axes.plot((0.0, 0.0), (U0, Y1), 'k--')\n",
    "#axes.text(10., 0.85, '$\\\\frac{K_1}{2}$', fontsize=18)\n",
    "#axes.plot((0.0, dt/2), (Y1, Y1), 'k--')\n",
    "#axes.text(250, Y1 - 0.05, '$\\\\frac{\\Delta t}{2}$', fontsize=18)\n",
    "\n",
    "# RK2 Step\n",
    "axes.plot(dt, U1, 'go')\n",
    "axes.plot((0., dt), (U0, U1), 'k')\n",
    "axes.text(dt+20, U1, '$U_{RK2}$', fontsize=18)\n",
    "#axes.plot((0.0, 0.0), (U0, U1), 'g--')\n",
    "#axes.plot((0.0, dt), (U1, U1), 'g--')\n",
    "\n",
    "# difference\n",
    "axes.plot((dt, dt), (U1, U0+K1),'k--')\n",
    "axes.text(dt+40, 0.5*(U1 + U0+K1), '$\\Delta\\propto\\Delta t^2$', fontsize=18)\n",
    "\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.05))\n",
    "axes.legend(loc='best')\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Relative Truncation Error \n",
    "\n",
    "If we consider the *Step Error* for each of our schemes, we know that\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    U^{euler}_{n+1}  &= u_{n+1} + O(\\Delta t^2)\\\\\n",
    "    U^{RK2}_{n+1} &= u_{n+1}  + O(\\Delta t^3)\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Therefore we can compute the *relative truncation error* as\n",
    "\n",
    "$$\n",
    "    \\Delta = | U^{euler}_{n+1} - U^{RK2}_{n+1} | = O(\\Delta t^{?})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* $\\Delta$ is Computable!\n",
    "* has a known dependence on step-size $\\Delta t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Adaptive Time Stepping\n",
    "\n",
    "Given the relative truncation error and its scaling with $\\Delta t$, we can now use this to choose a single good time step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example:\n",
    "\n",
    "Suppose we wanted our relative truncation error to be small relative to the solution or zero, we could set \n",
    "\n",
    "$$\n",
    "    \\Delta_{target} = \\mathtt{rtol}\\,U^{RK2}_{n+1} + \\mathtt{atol}\n",
    "$$\n",
    "\n",
    "where $\\mathtt{rtol}$ and $\\mathtt{atol}$ are relative and absolute tolerances (and we assume that $U^{RK2}_{n+1}$ is a reasonably good estimate of the true solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Moreover,  we know how the relative truncation error should scale with time step, i.e.\n",
    "\n",
    "$$\n",
    "    \\Delta_{target} \\propto \\Delta t_{target}^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "\n",
    "But our measured relative truncation error, $\\Delta$ depends on  the step size we just took i.e\n",
    "\n",
    "$$\n",
    "    \\Delta_{measured} \\propto \\Delta t_n^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Adaptive Time Stepping\n",
    "\n",
    "If we take the ratio of both relationships we get\n",
    "\n",
    "$$\n",
    "    \\frac{\\Delta_{target}}{\\Delta_{measured}} = \\left[\\frac{\\Delta t_{target}}{\\Delta t_{n}}\\right]^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or rearranging, we can solve for the target step size\n",
    "\n",
    "$$\n",
    "    \\Delta t_{target} = \\Delta t_{n}\\left[\\frac{\\Delta_{target}}{\\Delta_{measured}}\\right]^{\\frac{1}{2}}\n",
    "$$\n",
    "\n",
    "which tells us how to grow or shrink our time step to maintain accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In general,  if we have two methods with different step errors such that \n",
    "\n",
    "$$\n",
    "    \\Delta \\propto \\Delta t^p\n",
    "$$\n",
    "\n",
    "then our adaptive stepper will look like\n",
    "\n",
    "$$\n",
    "    \\Delta t_{target} = \\Delta t_{n}\\left[\\frac{\\Delta_{target}}{\\Delta_{measured}}\\right]^{1/p}\n",
    "$$\n",
    "\n",
    "This leads to all sorts of adaptive schemes most are included in standard libraries.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Embedded Runge-Kutta Schemes\n",
    "\n",
    "There are in fact a whole family of **Embedded RK** schemes which are $N$ stage schemes but can combine the $N$ function evaluations in two different ways to produce methods with different error estimates.  \n",
    "\n",
    "A popular one is **RK45** (available in `SciPy`) which is based on the [Dormand-Prince 5(4)](https://doi.org/10.1016/0771-050X(80)90013-3) pair which uses 6 function evaluations per step to produce a 4th order and 5th order scheme.  The 4th order scheme controls the time step, and the 5th order scheme actually is the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Dormand-Prince 5(4) Butcher Tableau\n",
    "\n",
    "$$\n",
    "\\begin{array} {c|ccccccc} \n",
    "0\\\\ \n",
    "1/5 & 1/5\\\\ \n",
    "3/10 &3/40 & 9/40 \\\\ \n",
    "4/5 & 44/45& & −56/15\t&32/9\\\\ \n",
    "8/9\t& 19372/6561\t& −25360/2187\t& 64448/6561 &\t−212/729\\\\\n",
    "1 & 9017/3168 &\t−355/33\t& 46732/5247 & 49/176& −5103/18656 \\\\\n",
    "1 & 35/384 & 0\t& 500/1113\t& 125/192\t&−2187/6784\t& 11/84\t\\\\ \\hline\n",
    "& 35/384 &\t0\t& 500/1113 &\t125/192 &\t−2187/6784& \t11/84\t&0\\\\\n",
    "& 5179/57600 &\t0\t& 7571/16695\t& 393/640\t& −92097/339200\t& 187/2100\t& 1/40\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "\n",
    "The first row of $b$ coefficients gives the fifth-order accurate solution and the second row gives the fourth-order accurate solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "def f_vanderpol(t, u, mu=5):\n",
    "    return numpy.array([u[1], mu * (1.0 - u[0]**2) * u[1] - u[0]])\n",
    "\n",
    "t_span = (0., 50.)\n",
    "u0 = [ 1., 0. ]\n",
    "f = lambda t, u : f_vanderpol(t, u, mu=20) \n",
    "sol = solve_ivp(f, t_span, u0, method='BDF',rtol=1.e-3,atol=1.e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "\n",
    "axes.plot(sol.t, sol.y[0],'o-')\n",
    "axes.set_title(\"Solution to Van der Pol Oscillator\", fontsize=18)\n",
    "axes.set_xlabel(\"t\", fontsize=16)\n",
    "axes.set_ylabel(\"y(t)\", fontsize=16)\n",
    "axes.grid()\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "delta_t = sol.t[1:] - sol.t[:-1]\n",
    "axes.plot(sol.t[:-1], delta_t)\n",
    "axes.grid()\n",
    "axes.set_xlabel('$t$', fontsize=16)\n",
    "axes.set_ylabel('$\\Delta t$', fontsize=16)\n",
    "axes.set_title('Timesteps, N = {}'.format(len(sol.t)), fontsize=18)\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 its 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",
    "$$\n",
    "    u^{(p)}(t_n) = f^{(p-1)}(t_n, u(t_n))\n",
    "$$\n",
    "\n",
    "leading to the method\n",
    "\n",
    "$$\n",
    "\\begin{align}\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)) \\\\\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",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "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*}$$\n",
    "\n",
    "With Step error $O(\\Delta t^3)$ and truncation error $T$, $O(\\Delta t^2)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example\n",
    "\n",
    "Let's use our simplest problem $u'(t) = \\lambda u$ with $f=\\lambda u$.  Therefore\n",
    "\n",
    "$$\\begin{align*}\n",
    "    f(t,u) &= \\lambda u\\\\\n",
    "    f'(t,u) &= \\lambda u' = \\lambda f = \\lambda^2 u\\\\\n",
    "    f''(t,u) &= \\lambda^2 u' = \\lambda^2 f = \\lambda^3 u\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "so a third order scheme would look like\n",
    "$$\n",
    "\\begin{align}\n",
    "    U(t_{n+1}) &= U(t_n)\\left[ 1 + \\lambda\\Delta t  + \\frac{(\\lambda\\Delta t)^2}{2} + \\frac{(\\lambda\\Delta t)^3}{6}\\right]+ O(\\Delta t^4)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def Taylor_3_flambda_u(lamda, t_span, u0, N):\n",
    "    \"\"\" implement constant step size  3rd order Taylor Series method for f(t,u) = \\lambda u\"\"\"\n",
    "    \n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    lambda_dt = lamda*(t[1] - t[0])\n",
    "    u = numpy.empty(t.shape)\n",
    "    u[0] = u0 \n",
    "    for (n, t_n) in enumerate(t[:-1]):\n",
    "        u[n+1] = u[n] * ( 1. + lambda_dt  + (lambda_dt**2)/2.  + (lambda_dt**3)/6.)\n",
    "    return t, u\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "lam = -1.\n",
    "t_span = [0., 5.]\n",
    "u0 = 1.\n",
    "\n",
    "f = lambda t,u : -u\n",
    "t_exact = numpy.linspace(t_span[0], t_span[1], 100)\n",
    "u_exact = u0*numpy.exp(-t_exact)\n",
    "\n",
    "N = 20\n",
    "t_taylor, u_taylor = Taylor_3_flambda_u(lam, t_span, u0, N)\n",
    "t_euler, u_euler = euler(f, t_span, u0, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t_exact,u_exact,'k',label='exact')\n",
    "axes.plot(t_euler, u_euler, 'ro', label='euler')\n",
    "axes.plot(t_taylor, u_taylor, 'bo', label='Taylor3')\n",
    "\n",
    "axes.grid()\n",
    "axes.set_xlabel('t', fontsize=16)\n",
    "axes.set_ylabel('u', fontsize=16)\n",
    "axes.legend(loc='best')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Some Drawbacks\n",
    "\n",
    "**Taylor Series methods**\n",
    " - require differentiating the given equation which can be cumbersome and difficult to implement\n",
    " - require a new routine for every $f$\n",
    " \n",
    "**General one-step/multi-stage methods**\n",
    "  - higher order methods often require a large number of evaluations of $f$ per time step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Overview -- so far\n",
    "\n",
    "So far we have discussed 3 basic techniques for integration of ODE IVP's\n",
    "\n",
    "* **Single-Step Multi-Stage** schemes (explicit and implicit)\n",
    "* **Taylor's Series** Methods\n",
    "* **Linear Multi-step** schemes (just started this)\n",
    "\n",
    "as well as \n",
    "* **truncation error** of each method (and it's relation to step-error)\n",
    "* **adaptive stepping** for Single-Step Schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In general Single-Step Multi-Stage methods (e.g. Embedded RK schemes) plus adaptive time stepping make for a very robust family of solvers. However there are some other classical schemes worth mentioning that have some advantages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear Multi-Step Methods\n",
    "\n",
    "**Multi-step methods**  are ODE methods that \n",
    " - require only *one* new function evaluation per time step  to work.  \n",
    " - reuse values and function evaluations at some number of previous time steps\n",
    " \n",
    "**Disadvantages over single step methods**\n",
    "\n",
    " - Methods are not self-starting, i.e. they require other methods to find the initial values\n",
    " - Difficult to adapt. The time step  Δ𝑡  in one-step methods can be changed at any time while multi-step methods this is much more complex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Simplest example:  The leap-frog method\n",
    "\n",
    "The leap-frog method is similar to Euler's method in that it uses the information from two-previous time steps to advance the problem.  We can write the problem as a centered first-derivative centered   at time $t_{n+1}, U_{n+1}$ i.e. \n",
    "\n",
    "$$\\frac{U_{n+2} - U_{n}}{2\\Delta t} = f(t_{n+1}, U_{n+1})$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "    U_{n+2} = U_{n} + 2\\Delta t\\, f(t_{n+1}, U_{n+1})\n",
    "$$ \n",
    "\n",
    "this method is known as the leap-frog method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "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(figsize=(8,6))\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",
    "u1 = c_0 * numpy.exp(decay_constant * dt / 2.0)\n",
    "u_np = c_0 + dt * (decay_constant * u1)\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",
    "axes.plot([0., dt/2, dt], [ c_0, u1, u_np],'ro')\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.grid()\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": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def leap_frog(f, t_span, u0, N, start=RK2):\n",
    "    \"\"\" calculate fixed step with leap-frog iterator with a single step starter\n",
    "    \"\"\"\n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    delta_t = t[1] - t[0]\n",
    "    u = numpy.zeros(t.shape)\n",
    "    u[0] = u0\n",
    "    # use a single-step multi-stage method to start\n",
    "    \n",
    "    t_start, u_start = start(f, (t[0],t[1]), u0, 2)\n",
    "    u[1] = u_start[-1]\n",
    "    for (n, t_np) in enumerate(t[1:-1]):\n",
    "        u[n+2] = u[n] +  2 *delta_t * f(t_np, u[n+1]) \n",
    "    return t, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "u0 = 1.0\n",
    "t_span = (0., 1600.) \n",
    "N = 7\n",
    "\n",
    "# Stable example\n",
    "decay_constant = -numpy.log(2.0) / 1600.0\n",
    "f = lambda t, u: decay_constant * u\n",
    "\n",
    "t_exact = numpy.linspace(t_span[0], t_span[1], N)\n",
    "u_exact = u0 * numpy.exp( decay_constant * t_exact)\n",
    "\n",
    "t_leapfrog, u_leapfrog = leap_frog(f, t_span, u0, N, start=RK4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(t_leapfrog, u_leapfrog, 'or-', label=\"Leap-Frog\")\n",
    "axes.plot(t_exact, u_exact, 'k--', label=\"True Solution\")\n",
    "axes.grid()\n",
    "axes.set_title(\"Leap-Frog\", fontsize=18)\n",
    "axes.set_xlabel(\"t (years)\", fontsize=16)\n",
    "axes.set_xlabel(\"$c(t)$\", fontsize=16)\n",
    "axes.legend(loc='best', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "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 around time $t_{n+1}$ to yield\n",
    "$$\\begin{aligned}\n",
    "    u_{n+2} &= u_{n+1} + \\Delta t f(t_{n+1},u_{n+1}) + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  + \\Delta t^3 \\frac{u'''(t_{n+1})}{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",
    "    U_{n+2} = U_{n} + 2 \\Delta t f(t_{n+1}, U_{n+1}).\n",
    "$$\n",
    "To handle the $U_{n}$ term we need to write this with relation to $u(t_{n+1})$.  Again we use the Taylor series\n",
    "$$\n",
    "    u(t_n) = u_{n+1} - \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  - \\Delta t^3 \\frac{u'''(t_{n+1})}{6} + \\mathcal{O}(\\Delta t^4)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    u(t_{n+2}) &= u_{n+1} + \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  + \\Delta t^3 \\frac{u'''(t_{n+1})}{6} + \\mathcal{O}(\\Delta t^4)\\\\\n",
    "    u(t_{n}) &= u_{n+1} - \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  - \\Delta t^3 \\frac{u'''(t_{n+1})}{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} + 2 \\Delta t f_{n+1}}_{U_{n+2}} - \\underbrace{\\left(u_{n+1} + \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  + \\Delta t^3 \\frac{u'''(t_{n+1})}{6} + \\mathcal{O}(\\Delta t^4) \\right )}_{u(t_{n+2})} \\right ] \\\\\n",
    "    &=\\frac{1}{\\Delta t} \\left [ \\underbrace{ \\left(u_{n+1} - \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  - \\Delta t^3 \\frac{u'''(t_{n+1})}{6} + \\mathcal{O}(\\Delta t^4)\\right)}_{u_{n}} + 2\\Delta t f_n - \\underbrace{\\left(u_{n+1} + \\Delta t f_{n+1} + \\Delta t^2 \\frac{u''(t_{n+1})}{2}  + \\Delta t^3 \\frac{u'''(t_{n+1})}{6} + \\mathcal{O}(\\Delta t^4) \\right )}_{u(t_{n+2})} \\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",
    "\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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between Euler, RK2 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_span = (0.0, 10.0)\n",
    "num_steps = [2**n for n in range(4,11)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "delta_t = numpy.empty(len(num_steps))\n",
    "error_euler = numpy.empty(len(num_steps))\n",
    "error_RK2 = 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(t_span[0], t_span[1], N)\n",
    "    tt, u_euler = euler(f, t_span, u_0, N )\n",
    "    tt, u_rk2 = RK2(f, t_span, u_0, N)\n",
    "    tt, u_leapfrog = leap_frog(f, t_span, u_0, N, start=euler)\n",
    "    \n",
    "    delta_t[i] = t[1] - t[0]\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_RK2[i] = numpy.linalg.norm(delta_t[i] * (u_rk2 - 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(figsize=(8,6))\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=1$')\n",
    "axes.loglog(delta_t, error_RK2, 'ro', label='RK2, $n=2$')\n",
    "axes.loglog(delta_t, error_leapfrog, 'go', label=\"Leap-Frog, $n=2$\")\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_RK2[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.grid()\n",
    "axes.legend(loc=2, fontsize=14)\n",
    "axes.set_title(\"Comparison of Errors\", fontsize=18)\n",
    "axes.set_xlabel(\"$\\Delta t$\",fontsize=16)\n",
    "axes.set_ylabel(\"$|U(t_f) - u(t_f)|$\", fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Look at the errors for Leap-Frog\n",
    "\n",
    "They're actually quite large...If you make a quick plot of u_leapfrog vs $t$ you'll see what's happening (and is a good example of an issue we will need to address in future lectures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N= 50\n",
    "t_leapfrog, u_leapfrog = leap_frog(f, t_span, u_0, N, start=euler)\n",
    "## Your plotting code here\n",
    "plt.figure()\n",
    "plt.plot(t_leapfrog, u_leapfrog)\n",
    "plt.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### General Linear Multi-Step Methods\n",
    "\n",
    "Leap-frog is perhaps the simplest of multi-step methods but 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.\n",
    "\n",
    "For example:  our Leap-frog method can be written using $r=2$, $\\alpha = \\begin{bmatrix} -1 & 0 & 1\\\\\n",
    "\\end{bmatrix}$, $\\beta = \\begin{bmatrix} 0 & 2 & 0 \\\\ \\end{bmatrix}$"
   ]
  },
  {
   "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 eliminate 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": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def AB2(f, t_span, u0, N, start=RK2):\n",
    "    \"\"\" calculate fixed step Adams-Bashforth 2-step method with a single step starter\n",
    "        reuses previous function evaluations\n",
    "    \"\"\"\n",
    "    t = numpy.linspace(t_span[0], t_span[1], N)\n",
    "    delta_t = t[1] - t[0]\n",
    "    u = numpy.zeros(t.shape)\n",
    "    \n",
    "    u[0] = u0\n",
    "    # use a single-step multi-stage method to start\n",
    "    t_start, u_start = start(f, (t[0],t[1]), u0, 2)\n",
    "    u[1] = u_start[-1]\n",
    "    \n",
    "    # set initial function evaluations\n",
    "    fn = f(t[0], u[0])\n",
    "    fnp = f(t[1], u[1])\n",
    "    for (n, t_np) in enumerate(t[2:]):\n",
    "        u[n+2] = u[n + 1] + delta_t / 2.0 * (-fn + 3.0 * fnp)\n",
    "        fn = fnp\n",
    "        fnp = f(t_np, u[n+2])\n",
    "    return t, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Use 2-step Adams-Bashforth to compute solution\n",
    "f = lambda t, u: -u\n",
    "t_span = (0., 10.)\n",
    "u0 = 1.0\n",
    "\n",
    "N = 20\n",
    "t, u_ab2 = AB2(f, t_span, u0, N, start=RK2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t_exact = numpy.linspace(t_span[0], t_span[1], 100)\n",
    "u_exact = numpy.exp(-t_exact)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\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\", fontsize=18)\n",
    "axes.set_xlabel(\"t\", fontsize=16)\n",
    "axes.set_ylabel(\"u(t)\",fontsize=16)\n",
    "axes.legend(loc=1, fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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(figsize=(8,6))\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 ($f=-u$)\", fontsize=18)\n",
    "axes.set_xlabel(\"t\", fontsize=16)\n",
    "axes.set_ylabel(\"u(t)\", fontsize=16)\n",
    "axes.legend(loc=1, fontsize=14)\n",
    "axes.grid()\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": "fragment"
    }
   },
   "source": [
    "Using the general expansion and evaluation 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}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "collecting terms of order $u^{(p)}$\n",
    "\n",
    "$$\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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Compare accuracy between RK-2, AB-2 and AM-2, RK-4\n",
    "f = lambda t, u: -u\n",
    "u_exact = lambda t: numpy.exp(-t)\n",
    "\n",
    "t_f = 10.0\n",
    "t_span = (0.0, t_f)\n",
    "\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_rk4 = 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",
    "    tt, U_rk  = RK2(f, t_span, u0, N)\n",
    "    \n",
    "    # Compute RK4\n",
    "    tt, U_rk4  = RK4(f, t_span, u0, N)\n",
    "        \n",
    "    # Compute Adams-Bashforth 2-stage\n",
    "    tt, U_ab  = AB2(f, t_span, u0, N)\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_rk4[i] = numpy.linalg.norm(delta_t[i] * (U_rk4 - 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(figsize=(8,6))\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",
    "axes.loglog(delta_t, error_rk4, 'co', label='RK-4')\n",
    "\n",
    "\n",
    "order_C = lambda delta_x, error, order: numpy.exp(numpy.log(error) - order * numpy.log(delta_x))\n",
    "\n",
    "colors = [ 'r', 'b', 'g', 'c']\n",
    "start = [ error_ab[0], error_ab[0], error_am[0], error_rk4[0] ]\n",
    "for i,c in enumerate(colors):\n",
    "    order = i+1.\n",
    "    axes.loglog(delta_t, order_C(delta_t[0], start[i], order) * delta_t**order, '--{}'.format(c),\n",
    "                label='$p={}$'.format(order))\n",
    "\n",
    "axes.legend(loc=4, fontsize=14)\n",
    "axes.set_title(\"Comparison of Errors\",fontsize=18)\n",
    "axes.set_xlabel(\"$\\Delta t$\",fontsize=16)\n",
    "axes.set_ylabel(\"$|U(t) - u(t)|$\", fontsize=16)\n",
    "axes.grid()\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 now 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 \\left[f(U_n) + f(\\hat{U}_{n+1}) \\right]\n",
    "\\end{aligned}$$\n",
    "leading to a second order accurate method.  Note this algorithm is identical to \\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\__________?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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 = 50\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",
    "t_ie, u_ieuler = improved_euler(f, (0.0, 10.), 1., N)\n",
    "fig = plt.figure(figsize=(8,6))\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",
    "axes.plot(t_ie, u_ieuler, 'bx', label=\"Improved Euler\")\n",
    "\n",
    "\n",
    "axes.set_title(\"Adams-Bashforth-Moulton P/C Method\", fontsize=18)\n",
    "axes.set_xlabel(\"t\", fontsize=18)\n",
    "axes.set_ylabel(\"u(t)\", fontsize=18)\n",
    "axes.legend(loc='best', fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}