Class_transcript_03_03_2020_ODEstiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In-class transcript from Lecture 16, March 3, 2020\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports and defs for lecture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are the standard imports for CS 111. \n",
    "# This list may change as the quarter goes on.\n",
    "\n",
    "import os\n",
    "import time\n",
    "import math\n",
    "import numpy as np\n",
    "import numpy.linalg as npla\n",
    "import scipy\n",
    "from scipy import linalg as spla\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg\n",
    "from scipy import integrate\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import axes3d\n",
    "%matplotlib inline\n",
    "import cs111\n",
    "np.set_printoptions(precision = 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radius of a match flame: A stiff ODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo of flame example with rk23\n",
    "\n",
    "def flamef(t, y):\n",
    "    \"\"\"ODE describing the radius of a flame\n",
    "    Input:\n",
    "      t is a scalar time\n",
    "      y is a vector of variables\n",
    "    Output:\n",
    "      ydot is the vector dy/dt\n",
    "      \n",
    "    The ODE describes the growth of a flame of radius y,\n",
    "    in terms of its oxygen balance. The flame takes in\n",
    "    oxygen at a rate proportional to its surface area,\n",
    "    and burns oxygen at a rate proportional to its volume:\n",
    "      dy/dt = y**2 - y**3\n",
    "    The flame balances at radius 1 (in these units), \n",
    "    so the exact solution increases monotonically from\n",
    "    an initial radius delta < 1 to the asymptotic value.\n",
    "    \"\"\"\n",
    "    ydot = y**2 - y**3\n",
    "    return ydot\n",
    "\n",
    "delta = .001\n",
    "tspan = (0, 2/delta)\n",
    "yinit = [delta]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = flamef, t_span = tspan, y0 = yinit, method = 'RK23')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('flame: ydot = y**2 - y**3,     %d steps' % len(sol.t))\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Zoom in on the flame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# zoom in on flame example with rk23\n",
    "\n",
    "delta = .001\n",
    "tspan = (0, 2/delta)\n",
    "yinit = [delta]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = flamef, t_span = tspan, y0 = yinit, method = 'RK23')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('flame: ydot = y**2 - y**3,     %d steps' % len(sol.t))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.gca().set_ylim([0, .2])\n",
    "plt.gca().set_xlim([850,1020])\n",
    "plt.title('zoom on initial growth')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.gca().set_ylim([0.99, 1.01])\n",
    "plt.gca().set_xlim([1000,1200])\n",
    "plt.title('zoom on final stage')\n",
    "print('number of steps:', len(sol.t))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Another example: The steepest-descent hiker "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo of steepest-descent hiker with rk23, starting on the river due west of town\n",
    "\n",
    "def valleyf(t, y):\n",
    "    \"\"\"ODE describing the path home in a valley in the fog\n",
    "    Input:\n",
    "      t is a scalar time\n",
    "      y is a vector of variables\n",
    "    Output:\n",
    "      ydot is the vector dy/dt\n",
    "      \n",
    "    You are somewhere in a long valley, whose altitude at\n",
    "    position (y0, y1) is a(y0, y1) =  (y0/1000)**2 + y1**2. \n",
    "\n",
    "    There is a town at the bottom of the valley, at (0,0).\n",
    "\n",
    "    To get to town, you want to walk downhill, \n",
    "    which is in the direction of the negative gradient, \n",
    "\n",
    "      - [da/dy0 , da/dy1]  =  - [y0/500 , 2*y1].\n",
    "\n",
    "    You want to walk at constant speed, \n",
    "    so you scale ydot to have length 1.\n",
    "    \n",
    "    The differential equation is dy/dt = valleyf(t, y).\n",
    "    \"\"\"\n",
    "\n",
    "    gradient = np.array([y[0]/500, 2*y[1]])\n",
    "    ydot = - gradient / npla.norm(gradient)\n",
    "    return ydot\n",
    "\n",
    "tspan = (0, 10)\n",
    "yinit = [-10, 0]  # try [-10, 0] and [-10, 0.1]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = valleyf, t_span = tspan, y0 = yinit, method = 'RK23')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.title('valley walk, %d steps' % len(sol.t))\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo of steepest-descent hiker with rk23, starting just above the river\n",
    "\n",
    "def valleyf(t, y):\n",
    "    \"\"\"ODE describing the path home in a valley in the fog\n",
    "    Input:\n",
    "      t is a scalar time\n",
    "      y is a vector of variables\n",
    "    Output:\n",
    "      ydot is the vector dy/dt\n",
    "      \n",
    "    You are somewhere in a long valley, whose altitude at\n",
    "    position (y0, y1) is a(y0, y1) =  (y0/1000)**2 + y1**2. \n",
    "\n",
    "    There is a town at the bottom of the valley, at (0,0).\n",
    "\n",
    "    To get to town, you want to walk downhill, \n",
    "    which is in the direction of the negative gradient, \n",
    "\n",
    "      - [da/dy0 , da/dy1]  =  - [y0/500 , 2*y1].\n",
    "\n",
    "    You want to walk at constant speed, \n",
    "    so you scale ydot to have length 1.\n",
    "    \n",
    "    The differential equation is dy/dt = valleyf(t, y).\n",
    "    \"\"\"\n",
    "\n",
    "    gradient = np.array([y[0]/500, 2*y[1]])\n",
    "    ydot = - gradient / npla.norm(gradient)\n",
    "    return ydot\n",
    "\n",
    "tspan = (0, 10)\n",
    "yinit = [-10, 0.1]  # try [-10, 0] and [-10, 0.1]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = valleyf, t_span = tspan, y0 = yinit, method = 'RK23')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.title('valley walk, %d steps' % len(sol.t))\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# zoom in on steepest-descent hiker with rk23\n",
    "\n",
    "tspan = (0, 10)\n",
    "yinit = [-10, 0.1]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = valleyf, t_span = tspan, y0 = yinit, method = 'RK23')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.title('valley walk, %d steps' % len(sol.t))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.001, .007])\n",
    "plt.gca().set_xlim([-10,-9.8])\n",
    "plt.title('zoom on turn toward village')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.00002, .00002])\n",
    "plt.gca().set_xlim([-10,-8.])\n",
    "plt.title('zoom on valley floor')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.000005, .000005])\n",
    "plt.gca().set_xlim([-.08,-.073])\n",
    "plt.title('zoom on end of walk')\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An implicit method (the Radau algorithm) used for the hiker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo of steepest-descent hiker with Radau (implicit stiff solver)\n",
    "\n",
    "tspan = (0, 10)\n",
    "yinit = [-10, 0.1]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = valleyf, t_span = tspan, y0 = yinit, method = 'Radau')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.title('valley walk, %d steps' % len(sol.t))\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# zoom in on steepest-descent hiker with Radau \n",
    "\n",
    "tspan = (0, 10)\n",
    "yinit = [-10, 0.1]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = valleyf, t_span = tspan, y0 = yinit, method = 'Radau')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.title('valley walk, %d steps' % len(sol.t))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.001, .007])\n",
    "plt.gca().set_xlim([-10,-9.8])\n",
    "plt.title('zoom on turn toward village')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.00002, .00002])\n",
    "plt.gca().set_xlim([-10,-8.])\n",
    "plt.title('zoom on valley floor')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.y[0], sol.y[1], 'o-', label='ode solution')\n",
    "plt.xlabel('y0')\n",
    "plt.ylabel('y1')\n",
    "plt.gca().set_ylim([-.00002, .00002])\n",
    "plt.gca().set_xlim([-1,0.1])\n",
    "plt.title('zoom on end of walk')\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forward Euler on a simple example: ydot = -y/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo of ode1 with exp(-t/2)\n",
    "def f(t, y):\n",
    "    \"\"\"function to be integrated to solve an ODE or a system of ODEs\n",
    "    Input:\n",
    "      t is a scalar time\n",
    "      y is a vector of variables\n",
    "    Output:\n",
    "      ydot is the vector dy/dt\n",
    "    \"\"\"\n",
    "    ydot = -y/2\n",
    "    return ydot\n",
    "\n",
    "tspan = (0,5)\n",
    "yinit = [1]\n",
    "step_size = 1 # try 1, .5, .1, .01\n",
    "\n",
    "sol_t, sol_y = cs111.ode1(fun = f, t_span = tspan, y0 = yinit, h = step_size)\n",
    "\n",
    "%matplotlib inline\n",
    "plt.plot(sol_t, sol_y[0], 'o', label='ode solution')\n",
    "tt = np.linspace(0, 5, 100)\n",
    "plt.plot(tt, np.exp(-tt/2), label='exp(-t/2)')\n",
    "plt.legend()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('ydot = -y/2, h = %g' % step_size)\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Backward Euler, an implicit method for stiff systems "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simple example ydot = -y/2 , y(0) = 1, with forward and backward Euler\n",
    "# one step of size h from t0 = 0, y0 = 1 to t = h\n",
    "\n",
    "# this is the problem from the sample midterm\n",
    "\n",
    "y0 = 1\n",
    "t0 = 0\n",
    "h = 1\n",
    "print(\"t0 =\", t0, \", y0 =\", y0)\n",
    "print(\"h =\", h)\n",
    "\n",
    "# exact solution: y = e^(-t/2)\n",
    "t = t0 + h\n",
    "print(\"at t =\", t, \"exact y =\", np.e**(-t/2))\n",
    "\n",
    "# forward Euler: y_{n+1} = y_n + h*f(t_n, y_n) = y_n + h*(-y_n / 2) = y_n * (1 - h/2)\n",
    "t = t0 + h\n",
    "y = y0 + h * (-y0/2)\n",
    "print(\"at t =\", t, \"forward Euler y =\", y)\n",
    "\n",
    "# backward Euler: \n",
    "# y_{n+1} = y_n + h*f(t_{n+1},y_{n+1}) = y_n + h*(-y_{n+1}/2),\n",
    "# so (1 + h/2)*y_{n+1} = y_n, \n",
    "# and y_{n+1} = y_n / (1 + h/2)\n",
    "t = t0 + h\n",
    "y = y0 / (1 + h/2)\n",
    "print(\"at t =\", t, \"backward Euler y =\", y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare forward and backward Euler for the simple ODE ydot = -y/2\n",
    "\n",
    "step_size = 1\n",
    "\n",
    "tspan = (0,5)\n",
    "yinit = [1]\n",
    "%matplotlib inline\n",
    "\n",
    "# forward Euler as implemented in cs111.ode1():\n",
    "sol_t, sol_y = cs111.ode1(fun = f, t_span = tspan, y0 = yinit, h = step_size)\n",
    "plt.plot(sol_t, sol_y[0], 'o', label='forward Euler')\n",
    "\n",
    "# backward Euler ONLY for this ODE:\n",
    "# y_{n+1} = y_n + h*f(t_{n+1},y_{n+1}) = y_n + h*(-y_{n+1}/2),\n",
    "# so (1 + h/2)*y_{n+1} = y_n, \n",
    "# and y_{n+1} = y_n / (1 + h/2)\n",
    "t = tspan[0]\n",
    "y = yinit[0]\n",
    "back_y = [y]\n",
    "while t < tspan[1]:\n",
    "    y = y / (1 + step_size/2)\n",
    "    back_y.append(y)\n",
    "    t = t + step_size\n",
    "plt.plot(sol_t, back_y, 'x', label='backward Euler')\n",
    "\n",
    "# the exact solution, which for this particular ODE happens to be easy to find\n",
    "tt = np.linspace(tspan[0], tspan[1], 100)\n",
    "plt.plot(tt, np.exp(-tt/2), label='exp(-t/2)')\n",
    "\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('ydot = -y/2, h = %g' % step_size)\n",
    "plt.legend()\n",
    "\n",
    "print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The flame ODE with the Radau implicit algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flame ODE, with zoom, using Radau\n",
    "\n",
    "delta = .001\n",
    "tspan = (0, 2/delta)\n",
    "yinit = [delta]\n",
    "\n",
    "sol = integrate.solve_ivp(fun = flamef, t_span = tspan, y0 = yinit, method = 'Radau')\n",
    "\n",
    "%matplotlib inline\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('flame: ydot = y**2 - y**3,     %d steps' % len(sol.t))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.gca().set_ylim([0, .2])\n",
    "plt.gca().set_xlim([850,1020])\n",
    "plt.title('zoom on initial growth')\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(sol.t, sol.y[0], '.-', label='ode solution')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.gca().set_ylim([0.99, 1.01])\n",
    "plt.gca().set_xlim([1000,1200])\n",
    "plt.title('zoom on final stage')\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}