{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f85212f7-7727-4c56-814a-ffc502db98d5",
   "metadata": {},
   "source": [
    "# Euler's method for ODEs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf38ae6-3db1-445c-a1e0-29b8afc0e577",
   "metadata": {},
   "source": [
    "\n",
    "We are interested in the numerical solution of the initial value problem (IVP) for an ordinary differential equation:\n",
    "\n",
    "$$ \\frac{\\mathrm{d}y}{\\mathrm{d}t} = f(t, y), \\quad a \\le t \\le b, \\quad y(a) = y_1.$$\n",
    "\n",
    "The idea is to start from $t = a$ (since we know $y(a)$), increment $t$\n",
    "by sufficiently small integration step $h$, and use the differential equation\n",
    "to determine $y(t + h)$. The process is then repeated until we reach\n",
    "$t = b$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c26a7ad3-5aa2-4159-966c-a9f93bb54db4",
   "metadata": {},
   "source": [
    "### Notations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "720b1f03-fa65-4917-88a9-19d9115cc75a",
   "metadata": {},
   "source": [
    "\n",
    "We denote the value of independent variable at the $i$th integration\n",
    "step by $t_{i+1}$, $i = 1, 2, \\ldots$, $t_1 = a$; the computed\n",
    "solution at the $i$th step by $y_{i+1}$,\n",
    "\n",
    "$$y_{i+1} \\equiv y(t_{i+1}) , \\quad i = 1, \\ldots, n-1;$$\n",
    "\n",
    "the value of the right hand side of the differential equation at the $i$th\n",
    "integration step by $f_{i+1}$,\n",
    "\n",
    "$$f_{i+1} \\equiv f(t_{i+1}, y_{i+1}) .$$\n",
    "\n",
    "The step size $h$ (assumed to be a constant for the sake of\n",
    "simplicity) is:\n",
    "\n",
    "$$h = t_i - t_{i-1} = \\frac{b - a}{n - 1} .$$\n",
    "\n",
    "Here $n-1$ is the total number of integration steps (corresponding to\n",
    "$n$ function evaluations of the right hand side of the ode."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cad8bdea-1697-40ac-a38f-a84d70d15a2e",
   "metadata": {},
   "source": [
    "\n",
    "## The method\n",
    "\n",
    "The Taylor series expansion of $y(t_{j+1})$ about $t_j$ correct up to\n",
    "the $h^2$ term is as following,\n",
    "\n",
    "$$y(t_{j+1}) = y(t_j + h) =\n",
    "  y(t_j) + h \\left.\\frac{\\mathrm{d}y}{\\mathrm{d}t}\\right|_{t_j} \n",
    "  + \\frac{h^2}{2} \\left.\\frac{\\mathrm{d}^2y}{\\mathrm{d}t^2}\\right|_{t_j} \n",
    "  + O(h^3).$$\n",
    "\n",
    "Using the differential equation for $\\frac{\\mathrm{dy}}{\\mathrm{d}t}$,\n",
    "\n",
    "$$y(t_{j+1}) =   y(t_j) + h \\, f(t_j, y_j) + \\alpha \\frac{h^2}{2} + O(h^3),$$\n",
    "\n",
    "or\n",
    "\n",
    "$$y_{j+1} = y_j + h f_j +  \\alpha \\frac{h^2}{2} + O(h^3),$$\n",
    "\n",
    "where $\\alpha$ is an unknown constant.\n",
    "\n",
    "Ignoring the quadratic and higher order terms, we obtain the\n",
    "expression for Euler's integration step:\n",
    "\n",
    "$$y_{j+1} = y_j + h f_j .$$\n",
    "\n",
    "In addition to deriving the expression for Euler's integration step, we learned that the leading\n",
    "in $h$ error term is quadratic in $h$, therefore Euler's method is a first order method."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79a04b95-ab27-45cc-8517-8b7ec4b1782e",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4116d1d-596f-408c-bbf8-ad4fbc0c8b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "    t, y = myeulers(fun, a, b, n, y1)\n",
    "\n",
    "Solve IVP y' = fun(t, y), a <= t <= b, y(a) = y1 using Euler's method.\n",
    "Use the integration step h = (b - a)/(n - 1). Return a vector of values\n",
    "of the independent variable t_i, and a vector of correspondinig values\n",
    "of the solution, y(t_i)\n",
    "\"\"\"\n",
    "function myeulers(fun, a, b, n, y1)\n",
    "    t = range(a, b, n)\n",
    "    y = zeros(n)\n",
    "    h = t[2] - t[1]\n",
    "    y[1] = y1\n",
    "    for i = 1:n-1\n",
    "        k1 = h*fun(t[i], y[i])\n",
    "        y[i+1] = y[i] + k1\n",
    "    end\n",
    "    return t, y\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec24c719-71c2-4f89-8e93-dadcad421127",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "a = 0.0\n",
    "b = 5.0\n",
    "fun(t, y) = exp(-sin(t)) - y * cos(t)\n",
    "y1 = 0.0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "143f8c42-c572-48ea-bc62-90ecd7f42cfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "yexact(t) = t * exp(-sin(t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb67e480-d11d-4378-8fbe-4e3117c109cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3cea235-73e8-4f2e-a652-b41eb26ede2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "n = 64\n",
    "t, y = myeulers(fun, a, b, n, y1)\n",
    "plot(t, y, label=\"Euler's\", marker=\".\")\n",
    "plot(t, yexact.(t), label=\"exact\")\n",
    "grid(true)\n",
    "legend()\n",
    "xlabel(\"t\")\n",
    "ylabel(\"y(t)\")\n",
    "title(\"Euler's method for IVP\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7cefcfe-d2e7-4b78-a0af-562f5536c60a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.5",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}