{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f85212f7-7727-4c56-814a-ffc502db98d5",
   "metadata": {},
   "source": [
    "\n",
    "# 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 for a first-order ordinary differential equation\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": [
    "\n",
    "## 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": [
    "\n",
    "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": "markdown",
   "id": "9392b08b-fe54-476a-afb0-e2dd7f3fdaee",
   "metadata": {},
   "source": [
    "\n",
    "## Global error of Euler's method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "022c4d42-74af-4d52-bf99-c3561fc9827c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndp = 9\n",
    "err = zeros(ndp)\n",
    "hs = zeros(ndp)\n",
    "ycheck = yexact(b);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef335953-3037-4216-9f9c-4724699621bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for l = 1:ndp\n",
    "    n = 2^(l+3)\n",
    "    _, y = myeulers(fun, a, b, n, y1)\n",
    "    err[l] = abs(ycheck - y[end])\n",
    "    hs[l] = (b - a)/(n - 1)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2cf154d-35fe-4fcc-aaed-8361903f106b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "loglog(hs, err, \".-\", label=\"experiment\")\n",
    "loglog(hs, hs .^ 1, \"--\", label=L\"h^1\")\n",
    "loglog(hs, hs .^ 2, \"--\", label=L\"h^2\")\n",
    "loglog(hs, hs .^ 3, \"--\", label=L\"h^3\")\n",
    "grid(true)\n",
    "legend()\n",
    "xlabel(\"integration step\")\n",
    "ylabel(\"absolute error\")\n",
    "title(\"Global error of Euler's method\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adefef10-f259-4ac8-9f7b-9c33d8c844b7",
   "metadata": {},
   "source": [
    "## Converting higher order ODEs to systems of first order ordinary differential equations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57b26f90-50d7-4b12-be20-07799b48fa1e",
   "metadata": {},
   "source": [
    "\n",
    "As an introductory example, let's consider an initial value problem for a second order ODE:\n",
    "\n",
    "$$\\frac{\\mathrm{d}^2 y}{\\mathrm{d} t^2} = f(t, y, \\frac{\\mathrm{d} y}{\\mathrm{d} t}), \\quad a \\le t \\le b, \\quad y(a) = y_a, \\quad\n",
    "\\frac{\\mathrm{d} y}{\\mathrm{d} t} (a) = y'_a.$$ \n",
    "\n",
    "We change the dependent variables:\n",
    "\n",
    "$$y_1 = y, \\quad y_2 = \\frac{\\mathrm{d} y}{\\mathrm{d} t}$$\n",
    "\n",
    "Then,\n",
    "\n",
    "$$\\frac{\\mathrm{d} y_1}{\\mathrm{d} t} = y_2,$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\frac{\\mathrm{d} y_2}{\\mathrm{d} t} = \\frac{\\mathrm{d}^2 y_1}{\\mathrm{d} t^2} = f(t, y, \\frac{\\mathrm{d} y}{\\mathrm{d} t}) = f(t, y_1, y_2).$$\n",
    "\n",
    "The last two equations for a system of two first order differential equations:\n",
    "$$\\frac{\\mathrm{d} y_1}{\\mathrm{d} t} = y_2,$$\n",
    "$$\\frac{\\mathrm{d} y_2}{\\mathrm{d} t} = f(t, y_1, y_2).$$\n",
    "\n",
    "The initial conditions for the system are as follows:\n",
    "\n",
    "$$y_1(a) = y_a, \\quad y_2(a) = y'_a.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ac247ea-ef90-46d0-b62a-86648cf1b536",
   "metadata": {},
   "source": [
    "## Euler's method for a system of first order ODEs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d679766c-412e-49cc-b1f0-cceacf87bdee",
   "metadata": {},
   "source": [
    "Let's introduce vector notations:\n",
    "\n",
    "$$Y = \\begin{pmatrix} y_1 \\\\ y_2 \\end{pmatrix}, \\quad F(t, Y) = \\begin{pmatrix} y_2 \\\\ f(t, y_1, y_2) \\end{pmatrix}.$$\n",
    "\n",
    "Then the system of equations can be written in a compact vector form:\n",
    "\n",
    "$$\\frac{\\mathrm{d} Y}{\\mathrm{d} t} = F(t, Y).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de832b8c-1464-41a5-adc2-e24eef121004",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67277d5f-e2fc-4dc6-8859-34b16797857b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "    t, y = myeulersv(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 matrix of correspondinig values\n",
    "of the solution, y(t_i)\n",
    "\"\"\"\n",
    "function myeulersv(fun, a, b, n, y1)\n",
    "    t = range(a, b, n)\n",
    "    neqs = length(y1)\n",
    "    y = zeros(neqs, 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": "markdown",
   "id": "94d294b2-893d-46ed-999e-7ff7e8ce8c7e",
   "metadata": {},
   "source": [
    "Consider the motion of harmonic oscillator with unit frequency:\n",
    "\n",
    "$$\\frac{\\mathrm{d}^2 y}{\\mathrm{d} t^2} = - y.$$\n",
    "\n",
    "In vector notations,\n",
    "\n",
    "$$Y = \\begin{pmatrix} y_1 \\\\ y_2 \\end{pmatrix}, \\quad F(t, Y) = \\begin{pmatrix} y_2 \\\\ -y_1 \\end{pmatrix}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b19efed-e3ef-4128-b5b0-3fe21ee0b9cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "a = 0.0\n",
    "b = 10.0\n",
    "n = 10000\n",
    "y1 = [1.0, 0.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e404180-94ab-4509-948a-ff9de21174f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fun(t, y) = [y[2], -y[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82bf5e06-2c76-40bb-a9bf-7d673a4b1bc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "t, y = myeulers(fun, a, b, n, y1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8cbeb0-609e-417d-9efe-326335d215ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plot(t, y[1, :], label=\"position\")\n",
    "plot(t, y[2, :], label=\"velocity\")\n",
    "grid(true)\n",
    "legend()\n",
    "xlabel(\"time\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d98f6cf-faae-4411-9847-367014ad4165",
   "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
}