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.0000000000000204 $$

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

In [8]:
Timing(gamma(3.0))
Out[8]:
$$ \left[ 0.003522827,2.0000000000000204 \right] $$
In [9]:
Timing(gamma1(3.0))
Out[9]:
$$ \left[ 0.000407314,2.0000000000000204 \right] $$

Version and date

In [10]:
VersionInfo()
Symata version     0.4.6
Julia version      1.6.0-DEV.58
Python version     3.8.3
SymPy version      1.5.1
In [11]:
InputForm(Now())
Out[11]:
2020-05-29T22:43:09.423