Van-der-PolOscillatorEx6Problem21.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "import plotly.graph_objects as go"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RK45; Adapted from https://stackoverflow.com/questions/54494770/how-to-set-fixed-step-size-with-scipy-integrate\n",
    "def DoPri45Step(f,t,x,h):\n",
    "\n",
    "    k1 = f(t,x) \n",
    "    k2 = f(t + 1./5*h, x + h*(1./5*k1) ) \n",
    "    k3 = f(t + 3./10*h, x + h*(3./40*k1 + 9./40*k2) ) \n",
    "    k4 = f(t + 4./5*h, x + h*(44./45*k1 - 56./15*k2 + 32./9*k3) ) \n",
    "    k5 = f(t + 8./9*h, x + h*(19372./6561*k1 - 25360./2187*k2 + 64448./6561*k3 - 212./729*k4) ) \n",
    "    k6 = f(t + h, x + h*(9017./3168*k1 - 355./33*k2 + 46732./5247*k3 + 49./176*k4 - 5103./18656*k5) )\n",
    "\n",
    "    v5 = 35./384*k1 + 500./1113*k3 + 125./192*k4 - 2187./6784*k5 + 11./84*k6\n",
    "    k7 = f(t + h, x + h*v5)\n",
    "    v4 = 5179./57600*k1 + 7571./16695*k3 + 393./640*k4 - 92097./339200*k5 + 187./2100*k6 + 1./40*k7\n",
    "\n",
    "    return v4,v5\n",
    "\n",
    "def DoPri45integrate(f, t, x0):\n",
    "    N=len(t)\n",
    "    x = [x0]\n",
    "    for k in range(N-1):\n",
    "        v4, v5 = DoPri45Step(f,t[k],x[k],t[k+1]-t[k])\n",
    "        x.append(\n",
    "            x[k] + (t[k+1]-t[k])*v5\n",
    "        )\n",
    "    return np.array(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Explicit Euler\n",
    "\n",
    "def ExplicitEulerStep(f, t, x, h):\n",
    "    return f(t, x)\n",
    "\n",
    "def ExplicitEulerIntegrate(f, t, x0):\n",
    "    N = len(t)\n",
    "    x = [x0]\n",
    "    for k in range(N-1):\n",
    "        v1 = ExplicitEulerStep(f, t[k], x[k], t[k+1]-t[k])\n",
    "        x.append(\n",
    "            x[k] + (t[k+1]-t[k])*v1\n",
    "        )\n",
    "    return np.array(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## PROBLEM DEFINITION\n",
    "# x = y[0], y = y[1]\n",
    "# Initial Conditions\n",
    "y0 = [1,0]\n",
    "# PARAMS\n",
    "epsilons = [0.16, 0.01, 0.0025, 0.0009, 0.000625]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SOLVER DEFINITION\n",
    "h = 1e-3\n",
    "ts = np.arange(0,10,h)\n",
    "solvers = [\n",
    "    {'func':ExplicitEulerIntegrate, 'name': 'Explicit Euler'}, \n",
    "    {'func':DoPri45integrate, 'name': 'Dormand-Prince-45 (Explicit RKM of Order 4(5))'}\n",
    "]\n",
    "\n",
    "solutions = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:5: RuntimeWarning: overflow encountered in double_scalars\n",
      "  \"\"\"\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:12: RuntimeWarning: invalid value encountered in add\n",
      "  if sys.path[0] == '':\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in subtract\n",
      "  if __name__ == '__main__':\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in subtract\n",
      "  # This is added back by InteractiveShellApp.init_path()\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in subtract\n",
      "  del sys.path[0]\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in add\n",
      "  del sys.path[0]\n",
      "C:\\Users\\NENNIA\\AppData\\Local\\Continuum\\miniconda3\\envs\\test\\lib\\site-packages\\ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in add\n",
      "  \"\"\"\n"
     ]
    }
   ],
   "source": [
    "for solver in solvers:\n",
    "    for epsilon in epsilons:\n",
    "        # ODE\n",
    "        f = lambda t,y: np.array([\n",
    "            (-y[1] + y[0]*(1 - y[0]**2))/epsilon,\n",
    "            y[0]\n",
    "        ])\n",
    "        solutions.append({\n",
    "            'ys':solver['func'](f, ts, y0).T, \n",
    "            'ts':ts,\n",
    "            'name':solver['name'],\n",
    "            'epsilon':epsilon,\n",
    "        })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epsilon:  0.16\n",
      "The solver successfully reached the end of the integration interval.\n",
      "Epsilon:  0.01\n",
      "The solver successfully reached the end of the integration interval.\n",
      "Epsilon:  0.0025\n",
      "The solver successfully reached the end of the integration interval.\n",
      "Epsilon:  0.0009\n",
      "The solver successfully reached the end of the integration interval.\n",
      "Epsilon:  0.000625\n",
      "The solver successfully reached the end of the integration interval.\n"
     ]
    }
   ],
   "source": [
    "# RKM 45, Wanted to use explicit Jacobian, but its too ill-conditioned\n",
    "for epsilon in epsilons:\n",
    "    print(\"Epsilon: \", epsilon)\n",
    "    # ODE\n",
    "    f = lambda t,y: np.array([\n",
    "        (-y[1] + y[0]*(1 - y[0]**2))/epsilon,\n",
    "        y[0]\n",
    "     ])\n",
    "    # Solution\n",
    "    sol = solve_ivp(\n",
    "        fun=f, \n",
    "        t_span=[0,10],\n",
    "        t_eval=ts,\n",
    "        y0=y0,\n",
    "        method='Radau',\n",
    "        dense_output=False,\n",
    "    )\n",
    "    print(sol.message)\n",
    "    \n",
    "    solutions.append({\n",
    "        'ts': sol.t,\n",
    "        'ys': sol.y,\n",
    "        'name':'Radau(Variable Step Size, Implicit RKM of Order 5)',\n",
    "        'epsilon':epsilon,\n",
    "    })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figs = []\n",
    "for s in solutions:\n",
    "    #print(s)\n",
    "    fig = go.Figure(\n",
    "        data=go.Scatter(\n",
    "            x=s['ys'][0], \n",
    "            y=s['ys'][1], \n",
    "            marker=dict(color=s['ts'], colorscale='Viridis', showscale=True),\n",
    "            mode='lines+markers',\n",
    "        )\n",
    "    )\n",
    "    fig.update_layout(\n",
    "        title=\"Solver: {0}, ε={1}\".format(s['name'], s['epsilon'])\n",
    "    )\n",
    "    \n",
    "    key = (s['name'], s['epsilon'])\n",
    "    fig.write_image(\"images/van-der-pol-%s-%s.pdf\"%key)\n",
    "    figs.append(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (test)",
   "language": "python",
   "name": "test"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}