{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:TITLE: On Schemes for Exponential Decay -->\n",
    "# On 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: **Sep 24, 2015**\n",
    "\n",
    "Copyright 2015, Hans Petter Langtangen. Released under CC Attribution 4.0 license\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/CN_logo.png, width=300 frac=0.3] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/CN_logo.png\" width=300>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "<!-- Each slide starts with !split and a title inside 5 = on each side -->\n",
    "<!-- (i.e., DocOnce subsections are used to specify slide titles, -->\n",
    "<!-- sections are used for parts/sections of the talk to appear in a -->\n",
    "<!-- table of contents) -->\n",
    "\n",
    "## Goal\n",
    "\n",
    "The primary goal of this demo talk is to demonstrate how to write\n",
    "talks with [DocOnce](https://github.com/hplgit/doconce)\n",
    "and get them rendered in numerous HTML formats.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Problem setting and methods\n",
    "<!-- Short title: Problem -->\n",
    "\n",
    "<!-- dom:FIGURE: [https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/method.png, width=600 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/method.png\" width=600>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "\n",
    "## We aim to solve the (almost) simplest possible differential equation 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)\n",
    "\\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\n",
    "\\label{initial:value} \\tag{2}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here,\n",
    "\n",
    " * $t\\in (0,T]$\n",
    "\n",
    " * $a$, $I$, and $T$ are prescribed parameters\n",
    "\n",
    " * $u(t)$ is the unknown function\n",
    "\n",
    " * The ODE [(1)](#ode) has the initial condition [(2)](#initial:value)\n",
    "\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/teacher2.jpg, width=250 frac=0.5] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/teacher2.jpg\" width=250>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "## The ODE problem is solved by a finite difference scheme\n",
    "\n",
    "\n",
    " * Mesh in time: $0= t_0< t_1 \\cdots < t_N=T$\n",
    "\n",
    " * Assume constant $\\Delta t = t_{n}-t_{n-1}$\n",
    "\n",
    " * $u^n$: numerical approx to the exact solution at $t_n$\n",
    "\n",
    "\n",
    "The $\\theta$ rule,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^{n+1} = \\frac{1 - (1-\\theta) a\\Delta t}{1 + \\theta a\\Delta t}u^n,\n",
    "\\quad n=0,1,\\ldots,N-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "contains the [Forward Euler](http://en.wikipedia.org/wiki/Forward_Euler_method) ($\\theta=0$),\n",
    "the [Backward Euler](http://en.wikipedia.org/wiki/Backward_Euler_method) ($\\theta=1$),\n",
    "and the [Crank-Nicolson](http://en.wikipedia.org/wiki/Crank-Nicolson) ($\\theta=0.5$)\n",
    "schemes.\n",
    "\n",
    "\n",
    "\n",
    "## The Forward Euler scheme explained\n",
    "\n",
    "\n",
    "<!-- dom:MOVIE: [http://youtu.be/PtJrPEIHNJw, width=640 height=480] -->\n",
    "<!-- begin movie -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "_s = \"\"\"\n",
    "<iframe width=\"640\" height=\"480\" src=\"http://www.youtube.com/embed/PtJrPEIHNJw\" frameborder=\"0\" allowfullscreen></iframe>\n",
    "\"\"\"\n",
    "HTML(_s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- end movie -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Implementation\n",
    "\n",
    "**Implementation in a Python function:**"
   ]
  },
  {
   "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]; step: dt.\"\"\"\n",
    "    dt = float(dt)           # avoid integer division\n",
    "    N = int(round(T/dt))     # no of time intervals\n",
    "    T = N*dt                 # adjust T to fit time step dt\n",
    "    u = zeros(N+1)           # array of u[n] values\n",
    "    t = linspace(0, T, N+1)  # time mesh\n",
    "\n",
    "    u[0] = I                 # assign initial condition\n",
    "    for n in range(0, N):    # n=0,1,...,N-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": [
    "## How to use the solver function\n",
    "\n",
    "\n",
    "**A complete main program.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Set problem parameters\n",
    "I = 1.2\n",
    "a = 0.2\n",
    "T = 8\n",
    "dt = 0.25\n",
    "theta = 0.5\n",
    "\n",
    "from solver import solver, exact_solution\n",
    "u, t = solver(I, a, T, dt, theta)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(t, u, t, exact_solution)\n",
    "plt.legend(['numerical', 'exact'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results\n",
    "\n",
    "<!-- dom:FIGURE: [https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/results.jpg, width=600 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/results.jpg\" width=600>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "## The Crank-Nicolson method shows oscillatory behavior for not sufficiently small time steps, while the solution should be monotone\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/CN.png, width=600 frac=0.9] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"https://raw.githubusercontent.com/hplgit/doconce/master/doc/pub/slides/demo/fig/CN.png\" width=600>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "## The artifacts can be explained by some theory\n",
    "\n",
    "Exact solution of the scheme:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u^n = A^n,\\quad A = \\frac{1 - (1-\\theta) a\\Delta t}{1 + \\theta a\\Delta t}\\thinspace .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Key results:\n",
    "\n",
    " * Stability: $|A| < 1$\n",
    "\n",
    " * No oscillations: $A>0$\n",
    "\n",
    " * $\\Delta t < 1/a$ for Forward Euler ($\\theta=0$)\n",
    "\n",
    " * $\\Delta t < 2/a$ for Crank-Nicolson ($\\theta=1/2$)\n",
    "\n",
    "\n",
    "**Concluding remarks:**\n",
    "\n",
    "Only the Backward Euler scheme is guaranteed to always give\n",
    "qualitatively correct results."
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 0
}