Numerical integration

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



In [1]:

    
using Symata;

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,5.999999999999999,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 export the Symata variable to Julia using ExportJ (or 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.007753122,2.0000000000000204 \right]  $$



In [9]:

    
Timing(gamma1(3.0))









    Out[9]:





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

Version and date



In [10]:

    
VersionInfo()









    



Symata version     0.4.1-dev.3
Julia version      0.7.0-beta2.1
Python version     2.7.14+
SymPy version      1.0



In [11]:

    
InputForm(Now())









    Out[11]:





2018-07-20T13:50:27.723