examples/4_custom_schemes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This notebook teaches how to:**  \n",
    "create custom integration schemes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Custom Integration Schemes\n",
    "\n",
    "In this tutorial we want to estimate $\\pi$ with the following equation:\n",
    "\n",
    "$\\pi = 4 \\int_0^1 \\sqrt{1-t^2}\\mathrm{d}t$\n",
    "\n",
    "We set up the simulation frame as explained in the previous tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe import Frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = Frame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Adding field for $\\pi$ and integration variable**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addfield(\"pi\", 0., description=\"Approximation of pi\")\n",
    "sim.addintegrationvariable(\"t\", 0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Differentiator**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def f(frame, x, Y):\n",
    "    return 4.*np.sqrt(1-x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.pi.differentiator = f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step size**\n",
    "\n",
    "We set the step size to 0.25, i.e., the integral function is only evaluated four times in the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dt(frame):\n",
    "    return 0.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t.updater = dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not need to write outputs for this model. We only have to tell the integrator when to stop the simulation, i.e., the upper bound of the integral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t.snapshots = [1.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setting the integrator**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe import Integrator\n",
    "from simframe import Instruction\n",
    "from simframe import schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator = Integrator(sim.t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator.instructions = [Instruction(schemes.expl_1_euler, sim.pi)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Running the simulation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Execution time: \u001b[94m0:00:00\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from IPython.display import Markdown as md\n",
    "def print_table(sim):\n",
    "    return md(\"| |$\\pi$|rel. error|\\n|-|-|-|\\n|real|{:10.8f}||\\n|approx.|{:10.8f}|{:9.3e}|\".format(np.pi,sim.pi,np.abs(np.pi-sim.pi)/np.pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "| |$\\pi$|rel. error|\n",
       "|-|-|-|\n",
       "|real|3.14159265||\n",
       "|approx.|3.49570907|1.127e-01|"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_table(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relative error is about 11 %.\n",
    "\n",
    "There are in principle two way how we can improve this estimate. One is to decrease the stepsize. You can rerun the notebook with a smaller step size and compare the results.\n",
    "\n",
    "Another method is to use more sophisticated integration schemes. We are now going to implement the 2nd-order [Midpoint Method](https://en.wikipedia.org/wiki/Midpoint_method) and the 4th-order [Runge-Kutta method](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods).\n",
    "\n",
    "## Writing custom integration schemes\n",
    "\n",
    "All we have to do is write a function that takes the current value of the integration variable **x0**, the current value of the variable to be integrated **Y0**, and the step size **dx**. The function needs to return the change **dY** of our variable **Y0** after the integration. We can access the derivative of the variable with `Y0.derivative(x, Y)`.\n",
    "\n",
    "For the midpoint method this looks as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def midpoint(x0, Y0, dx, *args, **kwargs):\n",
    "    x1 = x0 + 0.5*dx\n",
    "    Y1 = Y0 + 0.5*dx*Y0.derivative(x0, Y0)\n",
    "    return dx*Y0.derivative(x1, Y1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All we have to do now is to create an integration scheme from this function. This can be done by using `AbstractScheme` provided by `simframe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe.integration import Scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "expl_2_midpoint = Scheme(midpoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`expl_2_midpoint` is now our new integration scheme. The naming convention is `<expl/impl>_<order>_<name><_other>`.\n",
    "\n",
    "We can now assign a new instruction set using our new method to the integrator just as with the 1st-order Euler method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator.instructions = [Instruction(expl_2_midpoint, sim.pi)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we restart the simulation we have to reset to the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t = 0\n",
    "sim.pi = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Execution time: \u001b[94m0:00:00\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "| |$\\pi$|rel. error|\n",
       "|-|-|-|\n",
       "|real|3.14159265||\n",
       "|approx.|3.18392922|1.348e-02|"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_table(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error is now reduced to 1.3 % only by using a higher order method.\n",
    "\n",
    "Note: the higher order method needs more operations and is therefore slower. This does not really matter in this case, but might be important for more complex simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4th-order Runge-Kutta**\n",
    "\n",
    "The scheme function for the 4th-order explicit Runge-Kutta method looks as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rk4(x0, Y0, dx, *args, **kwargs):\n",
    "    k1 = Y0.derivative(x0         , Y0            )\n",
    "    k2 = Y0.derivative(x0 + 0.5*dx, Y0 + 0.5*dx*k1)\n",
    "    k3 = Y0.derivative(x0 + 0.5*dx, Y0 + 0.5*dx*k2)\n",
    "    k4 = Y0.derivative(x0 +     dx, Y0 +     dx*k3)\n",
    "    return dx*(1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "expl_4_rungekutta = Scheme(rk4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator.instructions = [Instruction(rk4, sim.pi)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t = 0\n",
    "sim.pi = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Execution time: \u001b[94m0:00:00\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "| |$\\pi$|rel. error|\n",
       "|-|-|-|\n",
       "|real|3.14159265||\n",
       "|approx.|3.12118917|6.495e-03|"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print_table(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this method the error is reduced down to 0.6 %."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Available integration schemes\n",
    "\n",
    "But before you take effort into developing your own integration schemes, take a look at the schemes already provided by `simframe`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "|Scheme|Description|\n",
       "|------|-----------|\n",
       "|expl_1_euler|Explicit adaptive 1st-order Euler method|\n",
       "|expl_2_fehlberg_adptv|Explicit adaptive 2nd-order Fehlberg's method|\n",
       "|expl_2_heun|Explicit 2nd-order Heun's method|\n",
       "|expl_2_heun_euler_adptv|Explicit adaptive 2nd-order Heun-Euler method|\n",
       "|expl_2_midpoint|Explicit 2nd-order midpoint method|\n",
       "|expl_2_ralston|Explicit 2nd-order Ralston's method|\n",
       "|expl_3_bogacki_shampine_adptv|Explicit adaptive 3rd-order Bogacki-Shampine method|\n",
       "|expl_3_gottlieb_shu_adptv|Explicit adaptive 3rd-order Gottlieb-Shu method|\n",
       "|expl_3_heun|Explicit 3rd-order Heun's method|\n",
       "|expl_3_kutta|Explicit 3rd-order Kutta's method|\n",
       "|expl_3_ralston|Explicit 3rd-order Ralston's method|\n",
       "|expl_3_ssprk|Explicit 3rd-order Strong Stability Preserving Runge-Kutta method|\n",
       "|expl_4_38rule|Explicit 4th-order 3/8 rule method|\n",
       "|expl_4_ralston|Explicit 4th-order Ralston's method|\n",
       "|expl_4_runge_kutta|Explicit 4th-order classical Runge-Kutta method|\n",
       "|expl_5_cash_karp_adptv|Explicit adaptive 5th-order Cash-Karp method|\n",
       "|expl_5_dormand_prince_adptv|Explicit adaptive 5th-order Dormand-Prince method|\n",
       "|impl_1_euler_direct|Implicit 1st-order direct Euler method|\n",
       "|impl_1_euler_gmres|Implicit 1st-order Euler method with GMRES solver|\n",
       "|impl_2_midpoint_direct|Implicit 2nd-order direct midpoint method|\n"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_  = \"|Scheme|Description|\\n\"\n",
    "_ += \"|------|-----------|\\n\"\n",
    "for s in schemes.__dir__():\n",
    "    if s == \"update\": continue\n",
    "    if not s.startswith(\"_\"): _ += \"|\"+s+\"|\"+schemes.__dict__[s]().description+\"|\\n\"\n",
    "md(_)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}