{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runge-Kutta 4th Order Method\n",
    "\n",
    "Our goal is to implement the fourth order Runge-Kutta method commonly refered to as RK4. We will compare its accuracy to those of the forward Euler and explicit trapezoidal methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import numpy.linalg\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def degradation_model( t, y, params ):\n",
    "    y1 = y[0]\n",
    "    y2 = y[1]\n",
    "    dydt = np.zeros( (2,) )\n",
    "    dydt[0] = -params[0] * y1\n",
    "    dydt[1] = -params[1] * y2\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "def forwardEuler( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler( fn, t, y0, params )\n",
    "    Given a function handle fn for dydt, find the value of\n",
    "    y at timepoints t, given y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[:,i-1],params);\n",
    "        y[:,i] = y[:,i-1] + h*yprime;\n",
    "    \n",
    "    return (t,y)\n",
    "\n",
    "def explicitTrapezoidal( fn, t, y0, params ):\n",
    "    \"\"\"[t,y] = explicitTrapezoidal( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (y0.size,t.size) )\n",
    "    y[:,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[:,i-1],params);\n",
    "        yhatFE = y[:,i-1] + h*yprime;\n",
    "        yprimehatFE = fn(t[i],yhatFE,params);\n",
    "        y[:,i] = y[:,i-1] + h/2*(yprime+yprimehatFE);\n",
    "\n",
    "    return (t,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implement the method\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "f_1 &=& f\\left(t,\\widehat{y}(t)\\right)\\\\\n",
    "f_2 &=& f\\left(t+\\frac{\\Delta t}{2}, \\widehat{y}(t) + \\frac{\\Delta t}{2} f_1\\right)\\\\\n",
    "f_3 &=& f\\left(t+\\frac{\\Delta t}{2}, \\widehat{y}(t) + \\frac{\\Delta t}{2} f_2\\right)\\\\\n",
    "f_4 &=& f\\left(t+\\Delta t, \\widehat{y}(t) + \\Delta t f_3 \\right)\\\\\n",
    "\\widehat{y}(t + \\Delta t)&=& \\widehat{y}(t) + \\frac{\\Delta t}{6} \\cdot (f_1 + 2f_2 + 2f_3 + f_4)\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK4( fn, t, y0, params ):\n",
    "    \"\"\"\n",
    "    [t,y] = RK4( fn, t, y0, params )\n",
    "    Given a function handle fn for y', find the value of\n",
    "    y at timepoints t, assuming y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    Uses a 4th order explicit Runge-Kutta method\n",
    "    \"\"\"\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run and plot result from RK4.\n",
    "\n",
    "Run the RK4 method on the degradation model with degradation rates (2,5) and initial conditions (1,1) with time beginning at 0 and ending at 2.\n",
    "\n",
    "Then plot the true solution (solid lines) and the forward Euler's solution (dashed lines), explicit Trapezoidal's solution (dotted lines), and RK4's solution (dash-dotted lines).\n",
    "\n",
    "Then plot the difference over time between the true solution and each of the other methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the error vs timestep size\n",
    "\n",
    "For 10 time step sizes ranging from $10^{-4}$ to 0.2, determine the error for forward Euler, explicit trapezoidal, and RK4. To compute the error, use the mean over time of the 2-norm at each time step. \n",
    "\n",
    "Plot the simulations in one figure and the error vs step size in another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Code\n",
    "### Examine substeps\n",
    "\n",
    "Examine the substeps taken to get the final slope estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fn = lambda t, y, params: -params[0]*y\n",
    "params = (5,)\n",
    "y0 = 1\n",
    "t0 = 0\n",
    "h = 0.2\n",
    "tru_t = np.linspace( t0, t0+h, 10 )\n",
    "tru_y = y0*np.exp(-params[0]*tru_t)\n",
    "\n",
    "f1 = fn(t0,y0,params);\n",
    "k1_t = t0+h/2\n",
    "k1_y = y0+0.5*h*f1\n",
    "f2 = fn(k1_t,k1_y,params);\n",
    "k2_t = t0+h/2\n",
    "k2_y = y0+0.5*h*f2\n",
    "f3 = fn(k2_t,k2_y,params);\n",
    "k3_t = t0+h\n",
    "k3_y = y0+h*f3\n",
    "f4 = fn(k3_t,k3_y,params);\n",
    "kf_y = y0 + h/6*(f1+2*f2+2*f3+f4);\n",
    "kf_t = t0 + h\n",
    "eps = h/10;\n",
    "\n",
    "plt.plot( tru_t, tru_y.T, 'C0' );\n",
    "plt.plot( [t0,k1_t], [y0,k1_y], 'C1--' );\n",
    "plt.plot( [k1_t-eps,k1_t+eps], [k1_y-eps*f2,k1_y+eps*f2], 'C2')\n",
    "plt.plot( [t0,k2_t], [y0,k2_y], 'C2--' );\n",
    "plt.plot( [k2_t-eps,k2_t+eps], [k2_y-eps*f3,k2_y+eps*f3], 'C3')\n",
    "plt.plot( [t0,k3_t], [y0,k3_y], 'C3--' );\n",
    "plt.plot( [k3_t-eps,k3_t+eps], [k3_y-eps*f4,k3_y+eps*f4], 'C5')\n",
    "plt.plot( [t0,kf_t], [y0,kf_y], 'C4--' );\n",
    "plt.plot( k1_t, k1_y, 'C1.')\n",
    "plt.plot( k2_t, k2_y, 'C2.')\n",
    "plt.plot( k3_t, k3_y, 'C3.')\n",
    "plt.plot( kf_t, kf_y, 'C4.')\n",
    "plt.legend( bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0., \n",
    "           labels=['True Y', 'Step 1','Slope for Step 2','Step 2','Slope for Step 3','Step 3','Slope to be used with others for Step 4','Step 4']);\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( 'Y');"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}