NIntegrate
performs numerical integration. The integrand can be a Symata expression or a compiled function
using Symata;
Compute the gamma function for a=3
. The result and estimated error are returned.
a = 3.0
NIntegrate(Exp(-t)*t^(a-1), [t,0,Infinity])
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
.
gamma(a_, x_:Infinity) := NIntegrate(Exp(-t)*t^(a-1), [t,0,x])[1]
Test the function on a List
of numbers.
Map(gamma, [3.0, 4.0, 5.0])
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.
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
gamma1(a1_, x_:Infinity) := (SetJ(a,a1), NIntegrate(gammaint, [0,x]))[1]
gamma1(3.0)
Timing shows that gamma1
is about 50 times faster than gamma
.
Timing(gamma(3.0))
Timing(gamma1(3.0))
VersionInfo()
Symata version 0.4.6 Julia version 1.6.0-DEV.58 Python version 3.8.3 SymPy version 1.5.1
InputForm(Now())
2020-05-29T22:43:09.423