{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "628f9256-8e52-4780-9d3e-b896e344b00c",
   "metadata": {},
   "source": [
    "\n",
    "# Numerical integration and Richardson extrapolation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c54ed57-45e1-4408-98d5-e0f1e7654873",
   "metadata": {},
   "source": [
    "\n",
    "Trapezoidal rule with the integration step $h$:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\int_a^b f(x) \\, \\mathrm{d}x \n",
    "    & =\n",
    "    \\frac{h}{2} \\left[\n",
    "    f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4)+ \\ldots + 2f(x_{n-1}) + f(x_n)\n",
    "    \\right] + \\alpha h^p  + \\beta h^q + \\gamma h^r + \\ldots \\\\\n",
    "    &= T(h, n) + \\alpha h^p  + \\beta h^q + \\gamma h^r + \\ldots, \n",
    "    \\tag{1}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "where $r > q > p > 0$ and for later use we introduced the notation\n",
    "\n",
    "$$\n",
    "    T(h, n) \\equiv \\frac{h}{2} \\left[\n",
    "    f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4)+ \\ldots + 2f(x_{n-1}) + f(x_n)\n",
    "    \\right] , \n",
    "    \\tag{2}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "    h = \\frac{b-a}{n - 1}, \\quad\n",
    "    x_i = a + (i - 1) \\, h, \\quad i = 1, \\ldots, n.\n",
    "$$\n",
    "\n",
    "Here $\\alpha h^p$ is the *leading term* of *discretization error* of the\n",
    "algorithm. The coefficients $\\alpha$, $\\beta$, $\\gamma$, $\\ldots$ are usually unknown. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd5007fa-a3f9-45a0-b8e4-ea380c4ecb98",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "117777a7-8308-4286-bcef-85b40c127c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "    ans = mytrapezoids(fun, a, b, n)\n",
    "\n",
    "Numerically evaluate the integral int_a^b fun(x) dx using the\n",
    "trapezoidal rule: I = h/2*(f_1 + 2f_2 + ... + 2f_{n-1} + f_n),\n",
    "where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, f_i = fun(x_i).\n",
    "\"\"\"\n",
    "function mytrapezoids(fun, a, b, n)\n",
    "    h = (b - a)/(n - 1)\n",
    "    s1 = fun(a) + fun(b)\n",
    "    s2 = 0.0\n",
    "    for i = 2:(n-1)\n",
    "        s2 += fun(a + (i-1)*h)\n",
    "    end\n",
    "    return (s1 + 2*s2)*h/2\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d46c6b5e-38e8-43df-85a6-70a0ccf76b11",
   "metadata": {},
   "source": [
    "\n",
    "### Testing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea5d85b4-3a7e-4de6-ac24-f4571bf719d9",
   "metadata": {},
   "source": [
    "\n",
    "#### Simple test #1: $\\quad \\displaystyle \\int_1^2 \\frac{\\mathrm{d}x}{x} = \\log(2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8599b3e-f07a-4fc8-9da3-912304bd60b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fun1(x) = 1/x\n",
    "a1 = 1.0\n",
    "b1 = 2.0\n",
    "np1 = 20\n",
    "exact1 = log(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a174a21c-6c73-4196-9470-8f6e89526acb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "round(abs(mytrapezoids(fun1, a1, b1, np1) - exact1), sigdigits=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f22d0653-62c7-49a0-b0b3-84bb0638ed5b",
   "metadata": {},
   "source": [
    "\n",
    "#### Simple test #2: $\\quad \\displaystyle \\int_0^1  e^x \\, \\mathrm{d}x = e - 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3e44f30-e8b1-48e9-b1ea-6b85b3bb4f00",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fun2(x) = exp(x)\n",
    "a2 = 0.0\n",
    "b2 = 1.0\n",
    "np2 = 100\n",
    "exact2 = exp(1) - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a748222-de64-4c10-8388-caf5d9d52977",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "round(abs(mytrapezoids(fun2, a2, b2, np2) - exact2), sigdigits=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b3f98df-c9ed-47fe-b9a9-a218da4d081f",
   "metadata": {},
   "source": [
    "\n",
    "### Error of trapezoidal formula"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5781bf57-ead9-4274-aa2c-ea907fbd38d8",
   "metadata": {},
   "source": [
    "\n",
    "A numerical experiment to determine the error of trapezoidal rule as a function of integration step $h$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef917639-f99f-41c1-a508-d1fede3ab84e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fun3(x) = sin(x)\n",
    "a3 = 0.0\n",
    "b3 = pi\n",
    "exact3 = 2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51bc1072-3ae0-42b3-9646-6b588b215e99",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ndp = 10           # number of tests\n",
    "errt = zeros(ndp)  # absolute errors\n",
    "hh = zeros(ndp);   # integration steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bb63042-e732-4adb-b0c8-079a3eddc792",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for k = 1:ndp \n",
    "    np = 2^k + 1\n",
    "    hh[k] = (b3 - a3)/(np - 1)\n",
    "    errt[k] = abs(mytrapezoids(fun3, a3, b3, np) - exact3)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "737191c8-704f-4f31-8f35-bf0c8b306a59",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "loglog(hh, errt, label=\"numerical experiment\", marker=\".\")\n",
    "p = 2\n",
    "loglog(hh, hh .^ p, linestyle=\"dashed\", label=L\"hh^{%$p}\")\n",
    "grid(true)\n",
    "xlabel(\"integration step\")\n",
    "ylabel(\"absolute error\")\n",
    "title(\"Error of trapezoidal rule\")\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f388c852-0232-42dc-b34b-06a5870bfea4",
   "metadata": {},
   "source": [
    "\n",
    "## Richardson extrapolation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b074d95-a332-468a-a91e-1ca3420792cb",
   "metadata": {},
   "source": [
    "\n",
    "Recall that \n",
    "$$\\int_a^b f(x) \\, \\mathrm{d}x = T(h, n) + \\alpha h^2  + \\beta h^q + \\ldots, \\tag{3}$$ \n",
    "where $T(h, n)$ is given by Eq. (2).\n",
    "\n",
    "Next, we write down the trapezoidal rule with step $2h$. It will involve\n",
    "only the nodes at $x_1$, $x_3$, $x_5$, $\\ldots$, $x_{n-2}$, $x_n$:\n",
    "$$\\int_a^b f(x) \\, \\mathrm{d}x = h [ f(x_1) + 2f(x_3) + 2f(x_5) + \\ldots + 2f(x_{n-2}) + f(x_n) ] + \\alpha (2h)^2  + \\beta (2h)^q + \\ldots \n",
    "= T(2h) + 4\\alpha h^2  + \\beta'h^q + \\ldots . \\tag{4}$$\n",
    "\n",
    "The leading error in Eq. (4) is 4 times the error in Eq. (3), since we doubled the step size. So we multiply both sides of Eq. (3) by 4,\n",
    "$$\n",
    "4\\int_a^b f(x) \\, \\mathrm{d}x = 4T(h) + 4\\alpha h^2 + \\beta''h^q + \\ldots,  \\tag{5}\n",
    "$$\n",
    "and subtract Eq. (4) from Eq. (5), intending to\n",
    "cancel out the leading error term. The result that we get is the\n",
    "*Simpson's formula*:\n",
    "$$\n",
    "\\int_a^b f(x) \\, \\mathrm{d}x  = \\frac{1}{3}\\left(4T(h) - T(2h) \\right) + \\beta'''h^q + \\ldots .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a64dac0c-4520-4427-b9b5-6ab62e41d452",
   "metadata": {},
   "source": [
    "### Simpson's formula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e499f89-d6af-4519-bf32-79242640a90b",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "    ans = mytrapezoids2(f, h)\n",
    "\n",
    "Numerically evaluate the integral int_a^b fun(x) dx using the\n",
    "trapezoidal rule: I = h/2*(f[1] + 2f[2] + ... + 2f[n-1] + f[n]),\n",
    "where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, f[i] = fun(x_i).\n",
    "\"\"\"\n",
    "function mytrapezoids2(fvec, h)\n",
    "    s1 = fvec[1] + fvec[end]\n",
    "    s2 = sum(fvec)\n",
    "    return (2*s2 - s1)*h/2\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c1f4b14-f1c6-4d00-8077-7fae9793eb96",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\"\"\"\n",
    "    ans = mysimpsons2(fun, a, b, n)\n",
    "\n",
    "Numerically evaluate the integral int_a^b fun(x) dx using the\n",
    "Simpson's rule, where h = (b - a)/(n - 1), x_i = a + (i - 1)*h, fvec[i] = fun(x_i).\n",
    "\"\"\"\n",
    "function mysimpsons2(fun, a, b, n)\n",
    "    x = range(a, b, n)\n",
    "    fvec = fun.(x)\n",
    "    h = x[2] - x[1]\n",
    "    t1 = mytrapezoids2(fvec, h)\n",
    "    t2 = mytrapezoids2(fvec[1:2:n], 2*h)\n",
    "    return (4*t1 - t2)/3\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b563b843-74d6-40b9-9a06-4882144c3c4a",
   "metadata": {},
   "source": [
    "\n",
    "#### Simple test #1: $\\quad \\displaystyle \\int_1^2 \\frac{\\mathrm{d}x}{x} = \\log(2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c5ac4ee-da7a-4f79-b2de-8670f37ec35d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "np1 = 21\n",
    "round(abs(mysimpsons2(fun1, a1, b1, np1) - exact1), sigdigits=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c8eb05-e873-4a28-97ad-231d1c0e890e",
   "metadata": {},
   "source": [
    "\n",
    "Approximately, how many nodes is necessary for trapezoidal rule to achieve the same precision?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b148128-2cf8-4231-b1b5-4731927d0de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "np1t = 501\n",
    "round(abs(mytrapezoids(fun1, a1, b1, np1t) - exact1), sigdigits=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1524120f-d1ba-4340-9a24-d9e301922fdb",
   "metadata": {},
   "source": [
    "\n",
    "#### Simple test #2: $\\quad \\displaystyle \\int_0^1  e^x \\, \\mathrm{d}x = e - 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fece0b7-3d0b-4240-adda-5005f224638d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "np2 = 101\n",
    "round(abs(mysimpsons2(fun2, a2, b2, np2) - exact2), sigdigits=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1fbecac-33f5-4afd-850e-9f1030ae030e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "np2t = 34001\n",
    "round(abs(mytrapezoids(fun2, a2, b2, np2t) - exact2), sigdigits=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff45c07b-879c-486e-bb84-9f8a2033e515",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.6",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}