vanderpol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vanderpol model\n",
    "\n",
    "We consider the following problem :\n",
    "\n",
    "\\begin{equation} \n",
    "\\left\\{ \n",
    "\\begin{aligned} \n",
    "{\\mathrm d}_t y_1 & = y_2\\\\ \n",
    "{\\mathrm d}_t y_2 & = \\varepsilon \\, (1 - y_1^2) \\, y_2  - y_1 \\quad \\text{avec } \\varepsilon > 0\n",
    "\\end{aligned} \n",
    "\\right. \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "from bokeh.io import push_notebook, show, output_notebook\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.layouts import column\n",
    "from bokeh.models import PrintfTickFormatter\n",
    "\n",
    "from ipywidgets import interact, FloatSlider\n",
    "\n",
    "from mylib.model import vanderpol_model\n",
    "import mylib.integration as integration\n",
    "\n",
    "import time\n",
    "\n",
    "output_notebook(hide_banner=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quasi exact solution\n",
    "\n",
    "The quasi-exact solution is obtained by using an implicit Runge-Kutta method of Radau IIA family of order 5 with stepsize control and fine tolerances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_quasi_exact_sol():\n",
    "\n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn_radau\n",
    "\n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "\n",
    "    yini = (0.5, 0)\n",
    "\n",
    "    sol_exa = integration.radau5(tini, tend, yini, fcn, 2, rtol=1.e-12, atol=1.e-12)\n",
    "\n",
    "    fig_sol_exa_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Quasi exact solution : y1\")\n",
    "    plt_y1 = fig_sol_exa_y1.line(sol_exa.t, sol_exa.y[0], line_width=2) \n",
    "\n",
    "    fig_sol_exa_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Quasi exact solution : y2\")\n",
    "    plt_y2 = fig_sol_exa_y2.line(sol_exa.t, sol_exa.y[1], line_width=2) \n",
    "\n",
    "    show(column(fig_sol_exa_y1, fig_sol_exa_y2), notebook_handle=True)\n",
    "\n",
    "    def update(eps) :\n",
    "        vdpm = vanderpol_model(eps)\n",
    "        fcn = vdpm.fcn_radau  \n",
    "        sol_exa = integration.radau5(tini, tend, yini, fcn, 2, rtol=1.e-12, atol=1.e-12)\n",
    "        plt_y1.data_source.data = dict(x=sol_exa.t, y=sol_exa.y[0])\n",
    "        plt_y2.data_source.data = dict(x=sol_exa.t, y=sol_exa.y[1])\n",
    "        #plt_pha.data_source.data = dict(x=sol_exa.y[0], y=sol_exa.y[1])\n",
    "        push_notebook()\n",
    "        \n",
    "    interact(update, eps=FloatSlider(min=1.,max=20.,step=1., value=1., continuous_update=False))\n",
    "    \n",
    "plot_quasi_exact_sol()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radau5 integration (python implementation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_radau5_sol():\n",
    "    \n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn\n",
    "    fcn_radau = vdpm.fcn_radau\n",
    "    jac = vdpm.jac\n",
    "    \n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "    \n",
    "    yini = (0.5, 0)\n",
    "    \n",
    "    tol = 1.e-6\n",
    "    t0 = time.time()\n",
    "    sol_radau = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=tol, atol=tol, jac=jac)\n",
    "    t1 = time.time()\n",
    "\n",
    "    sol_exa = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=1.e-12, atol=1.e-12, t_eval=sol_radau.t, jac=jac)\n",
    "    y1_err = np.abs(sol_exa.y[0] - sol_radau.y[0])\n",
    "    y2_err = np.abs(sol_exa.y[1] - sol_radau.y[1])\n",
    "\n",
    "    fig_sol_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y1\")\n",
    "    fig_sol_y1.x(sol_radau.t, sol_radau.y[0], line_width=2)  \n",
    "    \n",
    "    fig_sol_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y2\")\n",
    "    fig_sol_y2.x(sol_radau.t, sol_radau.y[1], line_width=2)  \n",
    "        \n",
    "    fig_err = figure(x_range=(tini, tend), y_axis_type=\"log\", plot_height=300, plot_width=950, title=\"Global error\")\n",
    "    fig_err.yaxis[0].formatter = PrintfTickFormatter(format=\"%8.1e\")\n",
    "    fig_err.x(sol_radau.t, y1_err, line_width=2, legend=\"y1\")\n",
    "    fig_err.x(sol_radau.t, y2_err, line_width=2, legend=\"y2\", color=\"crimson\")\n",
    "    fig_err.legend.location = \"bottom_right\"\n",
    "    \n",
    "    show(column(fig_sol_y1, fig_sol_y2, fig_err))\n",
    "    \n",
    "    print(f\"Time to integrate : {(t1-t0):.5f} s\")\n",
    "    print(f\"Number of accepted steps : {(sol_radau.t.size-1):d}\")\n",
    "    print(f\"Number of function evaluations : {sol_radau.nfev:d}\")\n",
    "    print(f\"Number of evaluations of the jacobian : {sol_radau.njev:d}\")\n",
    "    print(f\"Number of LU decompositions {sol_radau.nlu:d}\")\n",
    "\n",
    "plot_radau5_sol()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runge-Kutta \"3/8 method\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_rk38_sol():\n",
    "    \n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn\n",
    "    jac = vdpm.jac\n",
    "    \n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "    \n",
    "    yini = (0.5, 0)\n",
    "    \n",
    "    nt = 10000\n",
    "    t0 = time.time()\n",
    "    sol_rk38 = integration.rk38(tini, tend, nt, yini, fcn)\n",
    "    t1 = time.time()\n",
    "    \n",
    "    sol_exa = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=1.e-12, atol=1.e-12, t_eval=sol_rk38.t)\n",
    "    y1_err = np.abs(sol_exa.y[0] - sol_rk38.y[0])\n",
    "    y2_err = np.abs(sol_exa.y[1] - sol_rk38.y[1])\n",
    "\n",
    "    fig_sol_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y1\")\n",
    "    fig_sol_y1.x(sol_rk38.t, sol_rk38.y[0], legend=\"y1\", line_width=2)    \n",
    "    \n",
    "    fig_sol_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y2\")\n",
    "    fig_sol_y2.x(sol_rk38.t, sol_rk38.y[1], line_width=2)  \n",
    "\n",
    "    fig_err = figure(x_range=(tini, tend), y_axis_type=\"log\", plot_height=300, plot_width=950, title=\"Global error\")\n",
    "    fig_err.yaxis[0].formatter = PrintfTickFormatter(format=\"%8.1e\")\n",
    "    fig_err.x(sol_rk38.t, y1_err, line_width=2, legend=\"y1\")\n",
    "    fig_err.x(sol_rk38.t, y2_err, line_width=2, legend=\"y2\", color=\"crimson\")\n",
    "    fig_err.legend.location = \"bottom_right\"\n",
    "    \n",
    "    show(column(fig_sol_y1, fig_sol_y2, fig_err))\n",
    "\n",
    "    print(f\"Time to integrate : {(t1-t0):.5f} s\")\n",
    "    print(f\"Number of time steps : {(nt-1):d}\")\n",
    "    print(f\"Number of function evaluations : {(4*(nt-1)):d}\")\n",
    "\n",
    "plot_rk38_sol()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedded method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_rk_emb_sol():\n",
    "    \n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn\n",
    "    jac = vdpm.jac\n",
    "    \n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "    \n",
    "    yini = (0.5, 0)\n",
    "\n",
    "    tol = 1.e-8\n",
    "    t0 = time.time()\n",
    "    sol_rk_emb = integration.rk_embedded(tini, tend, yini, fcn, tol)\n",
    "    t1 = time.time()\n",
    "\n",
    "    sol_exa = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=1.e-12, atol=1.e-12, t_eval=sol_rk_emb.t)\n",
    "    y1_err = np.abs(sol_exa.y[0] - sol_rk_emb.y[0])\n",
    "    y2_err = np.abs(sol_exa.y[1] - sol_rk_emb.y[1])\n",
    "    \n",
    "    fig_sol_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y1\")\n",
    "    fig_sol_y1.x(sol_rk_emb.t, sol_rk_emb.y[0], legend=\"y1\", line_width=2)    \n",
    "    \n",
    "    fig_sol_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y2\")\n",
    "    fig_sol_y2.x(sol_rk_emb.t, sol_rk_emb.y[1], line_width=2)  \n",
    "\n",
    "    fig_err = figure(x_range=(tini, tend), y_axis_type=\"log\", plot_height=300, plot_width=950, title=\"Global error\")\n",
    "    fig_err.yaxis[0].formatter = PrintfTickFormatter(format=\"%8.1e\")\n",
    "    fig_err.x(sol_rk_emb.t, y1_err, line_width=2, legend=\"y1\")\n",
    "    fig_err.x(sol_rk_emb.t, y2_err, line_width=2, legend=\"y2\", color=\"crimson\")\n",
    "    fig_err.legend.location = \"bottom_right\"\n",
    "    \n",
    "    show(column(fig_sol_y1, fig_sol_y2, fig_err))\n",
    "\n",
    "    print(f\"Time to integrate : {(t1-t0):.5f} s\")\n",
    "    print(f\"Number of time steps : {(sol_rk_emb.t.size-1):d}\")\n",
    "    print(f\"Number of function evaluations : {sol_rk_emb.nfev:d}\")\n",
    "\n",
    "plot_rk_emb_sol()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dopri5 method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_dopri5_sol():\n",
    "    \n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn\n",
    "    jac = vdpm.jac\n",
    "    \n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "    \n",
    "    yini = (0.5, 0)\n",
    "    \n",
    "    tol = 1.e-8\n",
    "    t0 = time.time()\n",
    "    sol_dopri5 = solve_ivp(fcn, (tini, tend), yini, method=\"RK45\", rtol=tol, atol=tol)\n",
    "    t1 = time.time()\n",
    "    \n",
    "    sol_exa = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=1.e-12, atol=1.e-12, t_eval=sol_dopri5.t, jac=jac)\n",
    "    y1_err = np.abs(sol_exa.y[0] - sol_dopri5.y[0])\n",
    "    y2_err = np.abs(sol_exa.y[1] - sol_dopri5.y[1])\n",
    "        \n",
    "    fig_sol_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y1\")\n",
    "    fig_sol_y1.x(sol_dopri5.t, sol_dopri5.y[0], legend=\"y1\", line_width=2)    \n",
    "    \n",
    "    fig_sol_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y2\")\n",
    "    fig_sol_y2.x(sol_dopri5.t, sol_dopri5.y[1], line_width=2)  \n",
    "\n",
    "    fig_err = figure(x_range=(tini, tend), y_axis_type=\"log\", plot_height=300, plot_width=950, title=\"Global error\")\n",
    "    fig_err.yaxis[0].formatter = PrintfTickFormatter(format=\"%8.1e\")\n",
    "    fig_err.x(sol_dopri5.t, y1_err, line_width=2, legend=\"y1\")\n",
    "    fig_err.x(sol_dopri5.t, y2_err, line_width=2, legend=\"y2\", color=\"crimson\")\n",
    "    fig_err.legend.location = \"bottom_right\"\n",
    "\n",
    "    show(column(fig_sol_y1, fig_sol_y2, fig_err))\n",
    "    \n",
    "    print(f\"Time to integrate : {(t1-t0):.5f} s\")\n",
    "    print(f\"Number of time steps : {(sol_dopri5.t.size-1):d}\")\n",
    "    print(f\"Number of function evaluations : {sol_dopri5.nfev:d}\")\n",
    "    \n",
    "plot_dopri5_sol()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radau method (fortran implementation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_radau5_fortran_sol():\n",
    "    \n",
    "    vdpm = vanderpol_model(eps=1)\n",
    "    fcn = vdpm.fcn\n",
    "    fcn_radau = vdpm.fcn_radau\n",
    "    jac = vdpm.jac\n",
    "    \n",
    "    tini = 0. \n",
    "    tend = 100.\n",
    "    \n",
    "    yini = (0.5, 0)\n",
    "    \n",
    "    tol = 1.e-6\n",
    "    t0 = time.time()\n",
    "    sol_radau = integration.radau5(tini, tend, yini, fcn_radau, 2, rtol=tol, atol=tol)\n",
    "    t1 = time.time()\n",
    "    \n",
    "    sol_exa = solve_ivp(fcn, (tini, tend), yini, method=\"Radau\", rtol=1.e-12, atol=1.e-12, t_eval=sol_radau.t, jac=jac)\n",
    "    y1_err = np.abs(sol_exa.y[0] - sol_radau.y[0])\n",
    "    y2_err = np.abs(sol_exa.y[1] - sol_radau.y[1])\n",
    "\n",
    "    fig_sol = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution\")\n",
    "    fig_sol.x(sol_radau.t, sol_radau.y[0], legend=\"y1\", line_width=2)    \n",
    "    fig_sol.legend.location = \"top_left\"\n",
    "    \n",
    "    fig_sol_y1 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y1\")\n",
    "    fig_sol_y1.x(sol_radau.t, sol_radau.y[0], legend=\"y1\", line_width=2)    \n",
    "    \n",
    "    fig_sol_y2 = figure(x_range=(tini, tend), width=950, height=300, title=\"Solution y2\")\n",
    "    fig_sol_y2.x(sol_radau.t, sol_radau.y[1], line_width=2)  \n",
    "\n",
    "    fig_err = figure(x_range=(tini, tend), y_axis_type=\"log\", plot_height=300, plot_width=950, title=\"Global error\")\n",
    "    fig_err.yaxis[0].formatter = PrintfTickFormatter(format=\"%8.1e\")\n",
    "    fig_err.x(sol_radau.t, y1_err, line_width=2, legend=\"y1\")\n",
    "    fig_err.x(sol_radau.t, y2_err, line_width=2, legend=\"y2\", color=\"crimson\")\n",
    "    fig_err.legend.location = \"bottom_right\"\n",
    "    \n",
    "    show(column(fig_sol_y1, fig_sol_y2, fig_err))\n",
    "    \n",
    "    print(f\"Time to integrate : {(t1-t0):.5f} s\")    \n",
    "    print(f\"Number of function evalusations : {sol_radau.nfev:d}\")\n",
    "    print(f\"Number of jacobian evaluations : {sol_radau.njev:d}\")\n",
    "    print(f\"Number of computed steps : {sol_radau.nstep:d}\")\n",
    "    print(f\"Number of accepted steps : {sol_radau.naccpt:d}\")\n",
    "    print(f\"Number of rejected steps : {sol_radau.nrejct:d}\")\n",
    "    print(f\"Number of LU decompositions : {sol_radau.ndec:d}\")\n",
    "    print(f\"Number of forward-backward substitutions : {sol_radau.nsol:d}\")\n",
    "    \n",
    "plot_radau5_fortran_sol()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}