{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical integration\n",
    "\n",
    "`NIntegrate` performs numerical integration. The integrand can be a Symata expression or a compiled function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using Symata;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the gamma function for `a=3`.  The result and estimated error are returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$  \\left[ 2.0000000000000204,4.460650935347041e-9 \\right]  $$"
      ],
      "text/plain": [
       "L\"$$  \\left[ 2.0000000000000204,4.460650935347041e-9 \\right]  $$\""
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " a = 3.0\n",
    "NIntegrate(Exp(-t)*t^(a-1), [t,0,Infinity])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function `gamma`, where `gamma(a)` is the gamma function at `a` and `gamma(a,x)` is the lower incomplete gamma function. The default value of the second argument is `Infinity`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gamma(a_, x_:Infinity) := NIntegrate(Exp(-t)*t^(a-1), [t,0,x])[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the function on a `List` of numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$  \\left[ 2.0000000000000204,5.999999999999999,23.999999999999996 \\right]  $$"
      ],
      "text/plain": [
       "L\"$$  \\left[ 2.0000000000000204,5.999999999999999,23.999999999999996 \\right]  $$\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Map(gamma, [3.0, 4.0, 5.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integrand in the function `gamma` is a Symata expression that is re-evaluated every time the integration routine asks for a value. Using a compiled function in the integrand is more efficient.\n",
    "\n",
    "Create the compiled function like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " gammaint = Compile([t] , Exp(-t)*t^(a-1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`gammaint` refers to a Julia function of one variable, `x`. The variable `a` is free. To use `gammaint`, we have to export the Symata variable to Julia using `ExportJ` (or `SetJ`).\n",
    "\n",
    "The alternative gamma function is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gamma1(a1_, x_:Infinity) := (SetJ(a,a1), NIntegrate(gammaint, [0,x]))[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$ 2.0000000000000213 $$"
      ],
      "text/plain": [
       "L\"$$ 2.0000000000000213 $$\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma1(3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Timing shows that `gamma1` is about 50 times faster than `gamma`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$  \\left[ 0.007753122,2.0000000000000204 \\right]  $$"
      ],
      "text/plain": [
       "L\"$$  \\left[ 0.007753122,2.0000000000000204 \\right]  $$\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Timing(gamma(3.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$  \\left[ 0.000642127,2.0000000000000213 \\right]  $$"
      ],
      "text/plain": [
       "L\"$$  \\left[ 0.000642127,2.0000000000000213 \\right]  $$\""
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Timing(gamma1(3.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Version and date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symata version     0.4.1-dev.3\n",
      "Julia version      0.7.0-beta2.1\n",
      "Python version     2.7.14+\n",
      "SymPy version      1.0\n"
     ]
    }
   ],
   "source": [
    "VersionInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2018-07-20T13:50:27.723"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "InputForm(Now())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.7.0-beta2",
   "language": "julia",
   "name": "julia-0.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}