Numerical integration

NIntegrate performs numerical integration. The integrand can be a Symata expression or a compiled function



In [1]:

    
using Symata; isymata()

Compute the gamma function for a=3. The result and estimated error are returned.



In [2]:

    
( a = 3.0,
 NIntegrate(Exp(-t)*t^(a-1), [t,0,Infinity]))









    Out[2]:





$$  \left[ 2.0000000000000204,4.460650935347041e-9 \right]  $$

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.



In [3]:

    
gamma(a_, x_:Infinity) := NIntegrate(Exp(-t)*t^(a-1), [t,0,x])[1]

Test the function on a List of numbers.



In [4]:

    
Map(gamma, [3.0, 4.0, 5.0])









    Out[4]:





$$  \left[ 2.0000000000000204,6.0,23.999999999999996 \right]  $$

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.

Create the compiled function like this.



In [5]:

    
gammaint = Compile([t] , Exp(-t)*t^(a-1));

gammaint refers to a Julia function of one variable, x. The variable a is free. To use gammaint, we have to set the value of the Julia variable a using SetJ.

The alternative gamma function is



In [6]:

    
gamma1(a1_, x_:Infinity) := (SetJ(a,a1), NIntegrate(gammaint, [0,x]))[1]



In [7]:

    
gamma1(3.0)









    Out[7]:





$$ 2.0000000000000213 $$

Timing shows that gamma1 is about 50 times faster than gamma.



In [8]:

    
Timing(gamma(3.0))









    Out[8]:





$$  \left[ 0.010447708,2.0000000000000204 \right]  $$



In [9]:

    
Timing(gamma1(3.0))









    Out[9]:





$$  \left[ 0.000524856,2.0000000000000213 \right]  $$



In [ ]: