{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forward Euler Method Example\n",
    "\n",
    "Here we'll look at a simple implementation of the forward Euler method applied to a 1-dimensional function in order to see some if its properties. Our example will be the initial value problem defined as:\n",
    "\n",
    "$$ f^\\prime(t) = -0.1f(t), \\quad f(0) = 10 $$\n",
    "\n",
    "Of course, we can integrate this function to see that the solution is\n",
    "\n",
    "$$ f(t) = 10e^{-0.1t} $$\n",
    "\n",
    "We'll use this to compare the numerical solutions from the Euler method with the exact solution at $t=50$. The forward Euler method is:\n",
    "\n",
    "$$ f(t + \\Delta t) = f(t) + f^\\prime(t)\\,\\Delta t $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "#f' = -0.1*f\n",
    "def euler_step(f,dt):\n",
    "    return f - 0.1*f*dt\n",
    "\n",
    "##play with the step size number, observe how f changes\n",
    "dt = 1\n",
    "tf = 50.\n",
    "\n",
    "f0 = 10.\n",
    "t = 0.\n",
    "\n",
    "#put first value in array\n",
    "f = [f0]\n",
    "tv = [t]\n",
    "\n",
    "while t + dt <= tf:\n",
    "    f.append(euler_step(f[-1],dt))\n",
    "    t += dt\n",
    "    tv.append(t)\n",
    "\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the exact solution and the computed one\n",
    "\n",
    "tt = np.linspace(0,50,1000)\n",
    "\n",
    "fig,(ax1,ax2) = plt.subplots(1,2,figsize=(15,7))\n",
    "\n",
    "for ax in [ax1,ax2]:\n",
    "    for dt in [5,1,0.2]:\n",
    "        f0 = 10.\n",
    "        t = 0.\n",
    "\n",
    "        #put first value in array\n",
    "        f = [f0]\n",
    "        tv = [t]\n",
    "\n",
    "        while t + dt <= tf:\n",
    "            f.append(euler_step(f[-1],dt))\n",
    "            t += dt\n",
    "            tv.append(t)\n",
    "\n",
    "        error = f[-1] - 10*np.exp(-0.1*50)\n",
    "        ax.plot(tv,f,label=f'dt = {dt:.1f}, error = {error:.2e}',marker='o')\n",
    "\n",
    "    ax.plot(tt,10*np.exp(-0.1*tt),'k-',lw=2,label=r'$f(t) = 10e^{-0.1t}$')\n",
    "    ax.set_ylabel('f(t)')\n",
    "    ax.set_xlabel('t')\n",
    "    ax.legend()\n",
    "ax2.set_yscale('log')\n",
    "#fig.savefig('euler.png',bbox_inches='tight',transparent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the step size is too large, sometimes there can be catastrophic error. Consider this example, where\n",
    "\n",
    "$$ f^\\prime(t) = -f(t) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#f' = -f\n",
    "def euler_step_2(f,dt):\n",
    "    return f -f*dt\n",
    "\n",
    "tt = np.linspace(0,50,1000)\n",
    "\n",
    "fig,ax = plt.subplots(figsize=(12,10))\n",
    "\n",
    "\n",
    "for dt in [0.5,2,2.1]:\n",
    "    f0 = 10.\n",
    "    t = 0.\n",
    "\n",
    "    #put first value in array\n",
    "    f = [f0]\n",
    "    tv = [t]\n",
    "\n",
    "    while t + dt <= tf:\n",
    "        f.append(euler_step_2(f[-1],dt))\n",
    "        t += dt\n",
    "        tv.append(t)\n",
    "    \n",
    "    \n",
    "    error = f[-1] - 10*np.exp(-50)\n",
    "    ax.plot(tv,f,label=f'dt = {dt:.1f}, error = {error:.2e}',marker='o')\n",
    "\n",
    "ax.plot(tt,10*np.exp(-tt),'k-',lw=2,label='Exact')\n",
    "ax.set_ylabel('f(t)')\n",
    "ax.set_xlabel('t')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}