eulers_method_error.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the DE\n",
    "def f(t,y):\n",
    "    return -2*t*y**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Euler's Approx\n",
    "def euler_approx(t_0, t_end, num_steps, y_0):\n",
    "    \n",
    "    # calculate the time step\n",
    "    dt = (t_end-t_0)/num_steps\n",
    "    \n",
    "    \n",
    "    # calculate the euler approx at a given time step\n",
    "    def y_approx(t, y, dt):\n",
    "        return y + f(t, y) * dt\n",
    "    \n",
    "    # upper bound of second derivative \n",
    "    # along the graph of the solution\n",
    "    def M1(t_k, y_k):\n",
    "\n",
    "        # declare symbols\n",
    "        t,y = sympy.symbols('t y')\n",
    "\n",
    "        # partial with respect to t\n",
    "        dfdt = sympy.diff(f(t,y),t)\n",
    "\n",
    "        # partial with respect to y\n",
    "        dfdy = sympy.diff(f(t,y),y)\n",
    "\n",
    "        output = dfdt + dfdy*f(t,y)\n",
    "        return abs(output.evalf(subs = {t:t_k, y:y_k}))\n",
    "    \n",
    "    # bound for the parial derivative df/dy\n",
    "    def M2(t_k, y_k):\n",
    "        \n",
    "        # declare symbols\n",
    "        t,y = sympy.symbols('t y')\n",
    "\n",
    "        # partial with respect to y\n",
    "        dfdy = sympy.diff(f(t,y),y)\n",
    "        \n",
    "        return abs(dfdy.evalf(subs = {t:t_k, y:y_k}))\n",
    "\n",
    "    \n",
    "    # initialize lists\n",
    "    y = [y_0]\n",
    "    e = [0]\n",
    "    t = [t_0]\n",
    "    k = [0]\n",
    "    \n",
    "    for i in range(0, num_steps):\n",
    "        \n",
    "        # current time step\n",
    "        t_k = t_0+i*dt\n",
    "        \n",
    "        # approx the new value\n",
    "        if i == 0:\n",
    "            # approximate value\n",
    "            y_k = y_approx(t = t_0, y = y_0, dt = dt)\n",
    "            # approximate error\n",
    "            M1_k = M1(t_k, y_k)\n",
    "            M2_k = M2(t_k, y_k)\n",
    "            e_k = 0.5*M1_k*dt**2\n",
    "            \n",
    "        else:\n",
    "            # approximate value\n",
    "            y_k = y_approx(t = t_k, y = y_k, dt = dt)\n",
    "            # approximate error\n",
    "            M1_k = M1(t_k, y_k)\n",
    "            M2_k = M2(t_k, y_k)\n",
    "            e_k = (1 + M2_k*dt)*e_k + 0.5*M1_k*dt**2\n",
    "            \n",
    "        # append values to lists\n",
    "        y.append(y_k)\n",
    "        e.append(e_k)\n",
    "        t.append(t_k)\n",
    "        k.append(i+1)\n",
    "        \n",
    "    # plotting feature\n",
    "    plt.subplot(2, 1, 1)\n",
    "    plt.plot(t, y, '--')\n",
    "    plt.title('Euler Approximation and Error')\n",
    "    plt.ylabel('Approx Value')\n",
    "    plt.xlabel('time [s]')\n",
    "\n",
    "    plt.subplot(2, 1, 2)\n",
    "    plt.plot(k, e, 'r--')\n",
    "    plt.xlabel('Step Number')\n",
    "    plt.ylabel('Error')\n",
    "    plt.show()\n",
    "    \n",
    "    print(\"Estimated Value: \" + str(round(y[-1],6)))\n",
    "    print(\"Estimated Error: \" + str(round(e[-1],6)))\n",
    "    return {'value':y,'error':e,'time':t,'steps':k}\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run example\n",
    "euler_approx(t_0 = 0, t_end = 2, num_steps = 2000, y_0 = 1)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}