{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: Experiments with Schemes for Exponential Decay -->\n",
    "# Experiments with Schemes for Exponential Decay\n",
    "<!-- dom:AUTHOR: Hans Petter Langtangen Email:hpl@simula.no at Center for Biomedical Computing, Simula Research Laboratory & Department of Informatics, University of Oslo. -->\n",
    "<!-- Author: --> **Hans Petter Langtangen** (email: `hpl@simula.no`), Center for Biomedical Computing, Simula Research Laboratory and Department of Informatics, University of Oslo.\n",
    "\n",
    "Date: **Jul 2, 2016**\n",
    "\n",
    "Copyright 2016, Hans Petter Langtangen. Released under CC Attribution 4.0 license\n",
    "\n",
    "**Summary.** This report investigates the accuracy of three finite difference\n",
    "schemes for the ordinary differential equation $u'=-au$ with the\n",
    "aid of numerical experiments. Numerical artifacts are in particular\n",
    "demonstrated.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Mathematical problem\n",
    "<div id=\"math:problem\"></div>\n",
    "\n",
    "\n",
    "\n",
    "We address the initial-value problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"ode\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u'(t) = -au(t), \\quad t \\in (0,T], \\label{ode} \\tag{1}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"initial:value\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} \n",
    "u(0)  = I,                         \\label{initial:value} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $a$, $I$, and $T$ are prescribed parameters, and $u(t)$ is\n",
    "the unknown function to be estimated. This mathematical model\n",
    "is relevant for physical phenomena featuring exponential decay\n",
    "in time, e.g., vertical pressure variation in the atmosphere,\n",
    "cooling of an object, and radioactive decay.\n",
    "\n",
    "# Numerical solution method\n",
    "<div id=\"numerical:problem\"></div>\n",
    "\n",
    "\n",
    "\n",
    "We introduce a mesh in time with points $0 = t_0 < t_1 \\cdots < t_{N_t}=T$.\n",
    "For simplicity, we assume constant spacing $\\Delta t$ between the\n",
    "mesh points: $\\Delta t = t_{n}-t_{n-1}$, $n=1,\\ldots,N_t$. Let\n",
    "$u^n$ be the numerical approximation to the exact solution at $t_n$.\n",
    "\n",
    "The $\\theta$-rule [[Iserles_2009]](#Iserles_2009)\n",
    "is used to solve [(1)](#ode) numerically:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+1} = \\frac{1 - (1-\\theta) a\\Delta t}{1 + \\theta a\\Delta t}u^n,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $n=0,1,\\ldots,N_t-1$. This scheme corresponds to\n",
    "\n",
    "  * The [Forward Euler](http://en.wikipedia.org/wiki/Forward_Euler_method)\n",
    "    scheme when $\\theta=0$\n",
    "\n",
    "  * The [Backward Euler](http://en.wikipedia.org/wiki/Backward_Euler_method)\n",
    "    scheme when $\\theta=1$\n",
    "\n",
    "  * The [Crank-Nicolson](http://en.wikipedia.org/wiki/Crank-Nicolson)\n",
    "    scheme when $\\theta=1/2$\n",
    "\n",
    "# Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy import zeros, linspace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The numerical method is implemented in a Python function\n",
    "[[Langtangen_2014]](#Langtangen_2014) `solver` (found in the [`model.py`](http://bit.ly/29ayDx3) Python module file):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def solver(I, a, T, dt, theta):\n",
    "    \"\"\"Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt.\"\"\"\n",
    "    dt = float(dt)            # avoid integer division\n",
    "    Nt = int(round(T/dt))     # no of time intervals\n",
    "    T = Nt*dt                 # adjust T to fit time step dt\n",
    "    u = zeros(Nt+1)           # array of u[n] values\n",
    "    t = linspace(0, T, Nt+1)  # time mesh\n",
    "\n",
    "    u[0] = I                  # assign initial condition\n",
    "    for n in range(0, Nt):    # n=0,1,...,Nt-1\n",
    "        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]\n",
    "    return u, t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical experiments\n",
    "\n",
    "\n",
    "\n",
    "A set of numerical experiments has been carried out,\n",
    "where $I$, $a$, and $T$ are fixed, while $\\Delta t$ and\n",
    "$\\theta$ are varied. In particular, $I=1$, $a=2$,\n",
    "$\\Delta t = 1.25, 0.75, 0.5, 0.1$.\n",
    "[Figure](#fig:BE) contains four plots, corresponding to\n",
    "four decreasing $\\Delta t$ values. The red dashed line\n",
    "represent the numerical solution computed by the Backward\n",
    "Euler scheme, while the blue line is the exact solution.\n",
    "The corresponding results for the Crank-Nicolson and\n",
    "Forward Euler methods appear in Figures ref{fig:CN}\n",
    "and ref{fig:FE}.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [BE.png, width=800] The Backward Euler method for decreasing time step values. <div id=\"fig:BE\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fig:BE\"></div>\n",
    "\n",
    "<p>The Backward Euler method for decreasing time step values.</p>\n",
    "<img src=\"BE.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [CN.png, width=800] The Crank-Nicolson method for decreasing time step values. <div id=\"fig:CN\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fig:CN\"></div>\n",
    "\n",
    "<p>The Crank-Nicolson method for decreasing time step values.</p>\n",
    "<img src=\"CN.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [FE.png, width=800] The Forward Euler method for decreasing time step values. <div id=\"fig:FE\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fig:FE\"></div>\n",
    "\n",
    "<p>The Forward Euler method for decreasing time step values.</p>\n",
    "<img src=\"FE.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Error vs $\\Delta t$\n",
    "\n",
    "\n",
    "How the error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "E^n = \\left(\\int_0^T (Ie^{-at} - u^n)^2dt\\right)^{\\frac{1}{2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "varies with $\\Delta t$ for the three numerical methods\n",
    "is shown in [Figure](#fig:error).\n",
    "\n",
    "\n",
    "**Observe:**\n",
    "\n",
    "The data points for the three largest $\\Delta t$ values in the\n",
    "Forward Euler method are not relevant as the solution behaves\n",
    "non-physically.\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [error.png, width=400] Variation of the error with the time step. <div id=\"fig:error\"></div> -->\n",
    "<!-- begin figure -->\n",
    "<div id=\"fig:error\"></div>\n",
    "\n",
    "<p>Variation of the error with the time step.</p>\n",
    "<img src=\"error.png\" width=400>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "The $E$ numbers corresponding to [Figure](#fig:error)\n",
    "are given in the table below.\n",
    "\n",
    "<table border=\"1\">\n",
    "<thead>\n",
    "<tr><th align=\"center\">$\\Delta t$</th> <th align=\"center\">$\\theta=0$</th> <th align=\"center\">$\\theta=0.5$</th> <th align=\"center\">$\\theta=1$</th> </tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr><td align=\"right\">   1.25          </td> <td align=\"right\">   7.4630        </td> <td align=\"right\">   0.2161          </td> <td align=\"right\">   0.2440        </td> </tr>\n",
    "<tr><td align=\"right\">   0.75          </td> <td align=\"right\">   0.6632        </td> <td align=\"right\">   0.0744          </td> <td align=\"right\">   0.1875        </td> </tr>\n",
    "<tr><td align=\"right\">   0.50          </td> <td align=\"right\">   0.2797        </td> <td align=\"right\">   0.0315          </td> <td align=\"right\">   0.1397        </td> </tr>\n",
    "<tr><td align=\"right\">   0.10          </td> <td align=\"right\">   0.0377        </td> <td align=\"right\">   0.0012          </td> <td align=\"right\">   0.0335        </td> </tr>\n",
    "</tbody>\n",
    "</table>\n",
    "\n",
    "**Summary.**\n",
    "\n",
    "1. $\\theta =1$: $E\\sim \\Delta t$ (first-order convergence).\n",
    "\n",
    "2. $\\theta =0.5$: $E\\sim \\Delta t^2$ (second-order convergence).\n",
    "\n",
    "3. $\\theta =1$ is always stable and gives qualitatively corrects results.\n",
    "\n",
    "4. $\\theta =0.5$ never blows up, but may give oscillating solutions\n",
    "   if $\\Delta t$ is not sufficiently small.\n",
    "\n",
    "5. $\\theta =0$ suffers from fast-growing solution if $\\Delta t$ is\n",
    "   not small enough, but even below this limit one can have oscillating\n",
    "   solutions (unless $\\Delta t$ is sufficiently small).\n",
    "\n",
    "\n",
    "\n",
    "# Bibliography\n",
    "\n",
    "\n",
    " 1. <div id=\"Iserles_2009\"></div> **A. Iserles**. \n",
    "    *A First Course in the Numerical Analysis of Differential Equations*,\n",
    "    second edition,\n",
    "    Cambridge University Press,\n",
    "    2009.\n",
    "\n",
    " 2. <div id=\"Langtangen_2014\"></div> **H. P. Langtangen**. \n",
    "    *A Primer on Scientific Programming With Python*,\n",
    "    fifth edition,\n",
    "    Springer,\n",
    "    2016."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 0
}