{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PY 502 -  Numerical Integration \n",
    "### 09/22/2022\n",
    "#### Gabe Schumm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Integral of $sin(x)$ from $0$ to $\\pi$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function f(x)\n",
    "    return sin(x)\n",
    "end\n",
    "\n",
    "function F(x)\n",
    "    return -cos(x)\n",
    "end\n",
    "\n",
    "x0 = 0\n",
    "xf = pi\n",
    "\n",
    "I_exact = F(xf) - F(x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function trap(f, x0, xf, n)\n",
    "    h = (xf-x0)/n\n",
    "    \n",
    "    I = 0.5 * (f(x0) + f(xf)) #add end points\n",
    "    \n",
    "    for i=1:n-1\n",
    "        I += f(x0 + i*h)\n",
    "    end\n",
    "    return I*h\n",
    "end\n",
    "\n",
    "\n",
    "function simp(f, x0, xf, n)\n",
    "    h = (xf-x0)/n\n",
    "    \n",
    "    I = (f(x0) + f(xf)) # add end points\n",
    "    \n",
    "    for i=1:n-1 # n must be even\n",
    "        if mod(i, 2) == 1\n",
    "            I += 4 * f(x0 + i*h)\n",
    "        else\n",
    "            I += 2 * f(x0 + i*h)\n",
    "        end\n",
    "    end\n",
    "    return I*h/3\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loop over h values and compute trap/simps integrals, h = (xf - x0) / n\n",
    "I_trap = []\n",
    "I_simp = []\n",
    "hs = []\n",
    "\n",
    "for n=4:2:30 #only even n\n",
    "    push!(I_trap, trap(f, x0, xf, n))\n",
    "    push!(I_simp, simp(f, x0, xf, n))\n",
    "    push!(hs, (xf-x0)/n)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling of Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A log-log plot shows us power law scaling:\n",
    "\n",
    "$$\\left|I_{trap} - I_{exact}\\right| = \\Delta I_{trap} \\sim h^2  \\Rightarrow \\log\\left(\\Delta I_{trap}\\right)\\sim 2\\log\\left(h\\right)$$\n",
    "$$\\left|I_{simp} - I_{exact}\\right| = \\Delta I_{simp} \\sim h^4 \\Rightarrow \\log\\left(\\Delta I_{simp}\\right)\\sim 4\\log\\left(h\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(log.(hs), log.(abs.(I_exact.-I_trap) ),\n",
    "   markershape = :circle, markersize = 2,  label = \"Trapezoid\")\n",
    "plot!(log.(hs), log.(abs.(I_exact.-I_simp) ), xlabel = \"log(h)\", ylabel = \"log(∆I)\",\n",
    "   markershape = :circle, markersize = 2,  label = \"Simpson's\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate slope of log-log plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "∆x = log(hs[end]) - log(hs[1])\n",
    "\n",
    "trap_∆y = log(abs(I_exact-I_trap[end])) - log(abs(I_exact-I_trap[1]) )\n",
    "simp_∆y = log(abs(I_exact-I_simp[end])) - log(abs(I_exact-I_simp[1]) )\n",
    "\n",
    "\n",
    "println(\"Trapezoid slope = \", trap_∆y/∆x)\n",
    "println(\"Simpson's slope = \", simp_∆y/∆x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Romberg Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"h_0 = \",hs[1], \", I_0 = \",I_trap[1])\n",
    "println(\"h_1 = \",hs[3], \", I_1 = \",I_trap[3])\n",
    "println(\"h_1/h_0 = \", hs[3]/hs[1])\n",
    "romberg_simple = (4/3)*I_trap[3] - (1/3)*I_trap[1]\n",
    "println()\n",
    "println(\"Simple Romberg approximation (m = 1): I = \", romberg_simple)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = plot((hs), I_trap, xlabel = \"h\", ylabel = \"I\",\n",
    "   markershape = :circle, markersize = 2,  label = \"Trapezoid\")\n",
    "plot!((hs), I_simp,\n",
    "   markershape = :circle, markersize = 2,  label = \"Simpson's\")\n",
    "\n",
    "plot!([hs[1], hs[end]], [romberg_simple, romberg_simple], linestyle = :dash, c= :black, label = \"Simple Romberg\",\n",
    "    legend=:right)\n",
    "\n",
    "\n",
    "plot!(\n",
    "    (hs)[5:end], I_trap[5:end], markershape = :circle, markersize = 2,  label=\"\",\n",
    "    inset = (1, bbox(0.1, 0.15, 0.4, 0.4, :bottom, :left)),\n",
    "    xticks = round.(hs[5:3:end], digits = 2),\n",
    "    yticks = [1.997, 2],\n",
    "    subplot = 2\n",
    ")\n",
    "\n",
    "plot!(p[2], (hs)[5:end], I_simp[5:end], markershape = :circle, markersize = 2,  label=\"\")\n",
    "plot!(p[2], [hs[5], hs[end]], [romberg_simple, romberg_simple], linestyle = :dash, c= :black, label=\"\")\n",
    "#scatter!(p[2], [hs[end]], [2], c= :black, label=\"\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration using QuadGK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import Pkg; Pkg.add(\"QuadGK\")\n",
    "using QuadGK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "quadgk(f, x0, xf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function g(x)\n",
    "    return log(cos(exp(x))^2)  #1/(x-pi/2)^2\n",
    "end\n",
    "\n",
    "plot(x0:0.01:xf, g.(x0:.01:xf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_QGK, σ_QGK = quadgk(g, x0, xf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loop over h values and compute trap/simps integrals, h = (xf - x0) / n\n",
    "I_trap = []\n",
    "I_simp = []\n",
    "hs = []\n",
    "\n",
    "for n=400:200:5000 #only even n\n",
    "    push!(I_trap, trap(g, x0, xf, n))\n",
    "    push!(I_simp, simp(g, x0, xf, n))\n",
    "    push!(hs, (xf-x0)/n)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i, j = 4, 9\n",
    "println(\"h_0 = \",hs[i], \", I_0 = \",I_trap[i])\n",
    "println(\"h_1 = \",hs[j], \", I_1 = \",I_trap[j])\n",
    "println(\"h_1/h_0 = \", hs[j]/hs[i])\n",
    "romberg_simple = (4/3)*I_trap[j] - (1/3)*I_trap[i] ;\n",
    "println()\n",
    "println(\"Simple Romberg approximation (m = 1): I = \", romberg_simple)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = plot((hs), I_trap, xlabel = \"h\", ylabel = \"I\",\n",
    "   markershape = :circle, markersize = 2,  label = \"Trapezoid\")\n",
    "plot!((hs), I_simp,\n",
    "   markershape = :circle, markersize = 2,  label = \"Simpson's\")\n",
    "\n",
    "plot!([hs[1], hs[end]], [romberg_simple, romberg_simple], linestyle = :dash, c= :black, label = \"Simple Romberg\")\n",
    "\n",
    "plot!([hs[1], hs[end]], [I_QGK, I_QGK], linestyle = :dash, c= :red, label = \"QuadGK\",\n",
    "    legend=:right)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}