7. Integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to integrate functions numerically using the package QuadGK."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have not installed it yet, first add the package QuadGK:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General`\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m git-repo `https://github.com/JuliaRegistries/General.git`\n",
      "\u001b[?25l\u001b[2K\u001b[?25h\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.2/Project.toml`\n",
      "\u001b[90m [no changes]\u001b[39m\n",
      "\u001b[32m\u001b[1m  Updating\u001b[22m\u001b[39m `~/.julia/environments/v1.2/Manifest.toml`\n",
      "\u001b[90m [no changes]\u001b[39m\n"
     ]
    }
   ],
   "source": [
    "using Pkg\n",
    "Pkg.add(\"QuadGK\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrating a function is now relatively straightforward. As an example, I use\n",
    "$$\\int_0^1 x^2 \\,dx.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.3333333333333333, 5.551115123125783e-17)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using QuadGK\n",
    "\n",
    "#define function f\n",
    "f(x)=x^2\n",
    "\n",
    "#integrate f from 0 to 1\n",
    "quadgk(f,0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you see, the integration command \"quadgk\" returns two values. The first is the actual value of the integral that we wanted to calculate $\\int_0^1x^2\\,dx=1/3$. The second is an estimate on the absolute error, i.e. it tells you by how much julia might have gotten it wrong. For usual functions this error is small (as here) and you can ignore it. If the second number is big, then you might not want to trust the numerical result but this will rarely happen."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two more remarks: First, if you do not want to see the error estimate, you simply add \"[1]\" behind the quadgk command like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3333333333333333"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quadgk(f,0,1)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, you do not have to define the function over which you want to integrate in advance but you can use the quick way to define a so called anonymous function right in place:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.3333333333333333, 5.551115123125783e-17)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quadgk(x->x^2,0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced: Functions with discontinuities (or other known irregularities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you know that your function has a discontinuity at certain points, you can tell quadgk about it and this will help quadgk. The way you tell quadgk about it is that you add points of discontinuity between the boundaries of integration in the quadgk command. As an example, I calculate the integral $\\int_1^{-1} g(x)\\,dx$ where \n",
    "$$g(x)=\\begin{cases}1/\\sqrt{x}&\\text{ if }x>0\\\\0 & \\text{ if }x=0\\\\-1/\\sqrt{-x}&\\text{ else}\\end{cases}$$\n",
    "which has a discontinuity at zero. (Note that this is an option. In the most usual cases quadgk will do a fine job even if your function has a discontinuity and you do not tell quadgk about it.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.8554417678952652e-16, 4.011358198067816e-25)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function g(x)\n",
    "    if x>0.0\n",
    "        return 1.0/(sqrt(x))\n",
    "    elseif x==0\n",
    "        return 0.0\n",
    "    else\n",
    "        return -1.0/(sqrt(-x))\n",
    "    end\n",
    "end\n",
    "\n",
    "quadgk(g,-1.0,0.0,1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the symmetry of the function, the true result is, of course, zero but as you see we get pretty close to that (recall that \"e-16\" means \"multiply by $10^{-16}$\"). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced: Double integrals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are julia packages with algorithms specialized in multidimensional integrals but a simple way of evaluating thos is to apply quadgk twice. As an example, let us calculate\n",
    "$$\\int_0^1\\int_2^3 x/y \\,dy\\,dx.$$\n",
    "We will do so by defining a function which evaluates the inner integral for a given $x$, i.e. \n",
    "$$h(x)=\\int_2^3 x/y \\,dy.$$\n",
    "The original double integrall can then be written as $\\int_0^1 h(x)\\,dx$.\n",
    "\n",
    "The only thing one must not forget is to define $h$ such that it does not return the absolute error but only the value of the integral. As written above we do this by adding \"[1]\" at the end of the quadgk command. (By the way, here you can see that it is really handy to define the function $x/y$ \"anonymously\" right where we need it.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.2027325540540822, 0.0)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h(x)=quadgk(y->x/y,2,3)[1]\n",
    "\n",
    "quadgk(h,0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}