NAME¶
PDL::GSL::INTEG - PDL interface to numerical integration routines in GSL
DESCRIPTION¶
This is an interface to the numerical integration package present in the GNU
Scientific Library, which is an implementation of QUADPACK.
Functions are named
gslinteg_{algorithm} where {algorithm} is the
QUADPACK naming convention. The available functions are:
- gslinteg_qng: Non-adaptive Gauss-Kronrod integration
- gslinteg_qag: Adaptive integration
- gslinteg_qags: Adaptive integration with singularities
- gslinteg_qagp: Adaptive integration with known singular points
- gslinteg_qagi: Adaptive integration on infinite interval of the form
(-\infty,\infty)
- gslinteg_qagiu: Adaptive integration on infinite interval of the form
(a,\infty)
- gslinteg_qagil: Adaptive integration on infinite interval of the form
(-\infty,b)
- gslinteg_qawc: Adaptive integration for Cauchy principal values
- gslinteg_qaws: Adaptive integration for singular functions
- gslinteg_qawo: Adaptive integration for oscillatory functions
- gslinteg_qawf: Adaptive integration for Fourier integrals
Each algorithm computes an approximation to the integral, I, of the function
f(x)w(x), where w(x) is a weight function (for general integrands w(x)=1). The
user provides absolute and relative error bounds (epsabs,epsrel) which specify
the following accuracy requirement:
|RESULT - I| <= max(epsabs, epsrel |I|)
The routines will fail to converge if the error bounds are too stringent, but
always return the best approximation obtained up to that stage
All functions return the result, and estimate of the absolute error and an error
flag (which is zero if there were no problems). You are responsible for
checking for any errors, no warnings are issued unless the option {Warn =>
'y'} is specified in which case the reason of failure will be printed.
You can nest integrals up to 20 levels. If you find yourself in the unlikely
situation that you need more, you can change the value of
'max_nested_integrals' in the first line of the file 'FUNC.c' and recompile.
Please check the GSL documentation for more information.
SYNOPSIS¶
use PDL;
use PDL::GSL::INTEG;
my $a = 1.2;
my $b = 3.7;
my $epsrel = 0;
my $epsabs = 1e-6;
# Non adaptive integration
my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$a,$b,$epsrel,$epsabs);
# Warnings on
my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$a,$b,$epsrel,$epsabs,{Warn=>'y'});
# Adaptive integration with warnings on
my $limit = 1000;
my $key = 5;
my ($res,$abserr,$ierr) = gslinteg_qag(\&myf,$a,$b,$epsrel,
$epsabs,$limit,$key,{Warn=>'y'});
sub myf{
my ($x) = @_;
return exp(-$x**2);
}
FUNCTIONS¶
gslinteg_qng() -- Non-adaptive Gauss-Kronrod integration¶
This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and
87-point integration rules in succession until an estimate of the integral of
f over ($a,$b) is achieved within the desired absolute and relative error
limits, $epsabs and $epsrel. It is meant for fast integration of smooth
functions. It returns an array with the result, an estimate of the absolute
error, an error flag and the number of function evaluations performed.
Usage:
($res,$abserr,$ierr,$neval) = gslinteg_qng($function_ref,$a,$b,
$epsrel,$epsabs,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9);
# with warnings on
my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9,{Warn => 'y'});
sub f{
my ($x) = @_;
return ($x**2.6)*log(1.0/$x);
}
gslinteg_qag() -- Adaptive integration¶
This function applies an integration rule adaptively until an estimate of the
integral of f over ($a,$b) is achieved within the desired absolute and
relative error limits, $epsabs and $epsrel. On each iteration the adaptive
integration strategy bisects the interval with the largest error estimate; the
maximum number of allowed subdivisions is given by the parameter $limit. The
integration rule is determined by the value of $key, which has to be one of
(1,2,3,4,5,6) and correspond to the 15, 21, 31, 41, 51 and 61 point
Gauss-Kronrod rules respectively. It returns an array with the result, an
estimate of the absolute error and an error flag.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qag($function_ref,$a,$b,$epsrel,
$epsabs,$limit,$key,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1);
# with warnings on
my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1,{Warn => 'y'});
sub f{
my ($x) = @_;
return ($x**2.6)*log(1.0/$x);
}
gslinteg_qags() -- Adaptive integration with singularities¶
This function applies the Gauss-Kronrod 21-point integration rule adaptively
until an estimate of the integral of f over ($a,$b) is achieved within the
desired absolute and relative error limits, $epsabs and $epsrel. The algorithm
is such that it accelerates the convergence of the integral in the presence of
discontinuities and integrable singularities. The maximum number of allowed
subdivisions done by the adaptive algorithm must be supplied in the parameter
$limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qags($function_ref,$a,$b,$epsrel,
$epsabs,$limit,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
return ($x)*log(1.0/$x);
}
gslinteg_qagp() -- Adaptive integration with known singular points¶
This function applies the adaptive integration algorithm used by gslinteg_qags
taking into account the location of singular points until an estimate of the
integral of f over ($a,$b) is achieved within the desired absolute and
relative error limits, $epsabs and $epsrel. Singular points are supplied in
the piddle $points, whose endpoints determine the integration range. So, for
example, if the function has singular points at x_1 and x_2 and the integral
is desired from a to b (a < x_1 < x_2 < b), $points =
pdl(a,x_1,x_2,b). The maximum number of allowed subdivisions done by the
adaptive algorithm must be supplied in the parameter $limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qagp($function_ref,$points,$epsabs,
$epsrel,$limit,[{Warn => $warn}])
Example:
my $points = pdl(0,1,sqrt(2),3);
my ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
my $x2 = $x**2;
my $x3 = $x**3;
return $x3 * log(abs(($x2-1.0)*($x2-2.0)));
}
gslinteg_qagi() -- Adaptive integration on infinite interval¶
This function estimates the integral of the function f over the infinite
interval (-\infty,+\infty) within the desired absolute and relative error
limits, $epsabs and $epsrel. After a transformation, the algorithm of
gslinteg_qags with a 15-point Gauss-Kronrod rule is used. The maximum number
of allowed subdivisions done by the adaptive algorithm must be supplied in the
parameter $limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qagi($function_ref,$epsabs,
$epsrel,$limit,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000,{Warn => 'y'});
sub myfn{
my ($x) = @_;
return exp(-$x - $x*$x) ;
}
gslinteg_qagiu() -- Adaptive integration on infinite interval¶
This function estimates the integral of the function f over the infinite
interval (a,+\infty) within the desired absolute and relative error limits,
$epsabs and $epsrel. After a transformation, the algorithm of gslinteg_qags
with a 15-point Gauss-Kronrod rule is used. The maximum number of allowed
subdivisions done by the adaptive algorithm must be supplied in the parameter
$limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qagiu($function_ref,$a,$epsabs,
$epsrel,$limit,[{Warn => $warn}]);
Example:
my $alfa = 1;
my ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
if (($x==0) && ($alfa == 1)) {return 1;}
if (($x==0) && ($alfa > 1)) {return 0;}
return ($x**($alfa-1))/((1+10*$x)**2);
}
gslinteg_qagil() -- Adaptive integration on infinite interval¶
This function estimates the integral of the function f over the infinite
interval (-\infty,b) within the desired absolute and relative error limits,
$epsabs and $epsrel. After a transformation, the algorithm of gslinteg_qags
with a 15-point Gauss-Kronrod rule is used. The maximum number of allowed
subdivisions done by the adaptive algorithm must be supplied in the parameter
$limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qagl($function_ref,$b,$epsabs,
$epsrel,$limit,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000,{Warn => 'y'});
sub myfn{
my ($x) = @_;
return exp($x);
}
gslinteg_qawc() -- Adaptive integration for Cauchy principal values¶
This function computes the Cauchy principal value of the integral of f over
(a,b), with a singularity at c, I = \int_a^b dx f(x)/(x - c). The integral is
estimated within the desired absolute and relative error limits, $epsabs and
$epsrel. The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qawc($function_ref,$a,$b,$c,$epsabs,$epsrel,$limit)
Example:
my ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
return 1.0 / (5.0 * $x * $x * $x + 6.0) ;
}
gslinteg_qaws() -- Adaptive integration for singular functions¶
The algorithm in gslinteg_qaws is designed for integrands with
algebraic-logarithmic singularities at the end-points of an integration
region. Specifically, this function computes the integral given by I =
\int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x). The
integral is estimated within the desired absolute and relative error limits,
$epsabs and $epsrel. The maximum number of allowed subdivisions done by the
adaptive algorithm must be supplied in the parameter $limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) =
gslinteg_qawc($function_ref,$alpha,$beta,$mu,$nu,$a,$b,
$epsabs,$epsrel,$limit,[{Warn => $warn}]);
Example:
my ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
if($x==0){return 0;}
else{
my $u = log($x);
my $v = 1 + $u*$u;
return 1.0/($v*$v);
}
}
gslinteg_qawo() -- Adaptive integration for oscillatory functions¶
This function uses an adaptive algorithm to compute the integral of f over (a,b)
with the weight function sin(omega*x) or cos(omega*x) -- which of sine or
cosine is used is determined by the parameter $opt ('cos' or 'sin'). The
integral is estimated within the desired absolute and relative error limits,
$epsabs and $epsrel. The maximum number of allowed subdivisions done by the
adaptive algorithm must be supplied in the parameter $limit.
Please check the GSL documentation for more information.
Usage:
($res,$abserr,$ierr) = gslinteg_qawo($function_ref,$omega,$sin_or_cos,
$a,$b,$epsabs,$epsrel,$limit,[opt])
Example:
my $PI = 3.14159265358979323846264338328;
my ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000);
# with warnings on
($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
if($x==0){return 0;}
else{ return log($x);}
}
gslinteg_qawf() -- Adaptive integration for Fourier integrals¶
This function attempts to compute a Fourier integral of the function f over the
semi-infinite interval [a,+\infty). Specifically, it attempts tp compute I =
\int_a^{+\infty} dx f(x)w(x), where w(x) is sin(omega*x) or cos(omega*x) --
which of sine or cosine is used is determined by the parameter $opt ('cos' or
'sin'). The integral is estimated within the desired absolute error limit
$epsabs. The maximum number of allowed subdivisions done by the adaptive
algorithm must be supplied in the parameter $limit.
Please check the GSL documentation for more information.
Usage:
gslinteg_qawf($function_ref,$omega,$sin_or_cos,$a,$epsabs,$limit,[opt])
Example:
my ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000);
# with warnings on
($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000,{Warn => 'y'});
sub f{
my ($x) = @_;
if ($x == 0){return 0;}
return 1.0/sqrt($x)
}
BUGS¶
Feedback is welcome. Log bugs in the PDL bug database (the database is always
linked from <
http://pdl.perl.org>).
SEE ALSO¶
PDL
The GSL documentation is online at
http://www.gnu.org/software/gsl/manual/
AUTHOR¶
This file copyright (C) 2003,2005 Andres Jordan <ajordan@eso.org> All
rights reserved. There is no warranty. You are allowed to redistribute this
software documentation under certain conditions. For details, see the file
COPYING in the PDL distribution. If this file is separated from the PDL
distribution, the copyright notice should be included in the file.
The GSL integration routines were written by Brian Gough. QUADPACK was written
by Piessens, Doncker-Kapenga, Uberhuber and Kahaner.
FUNCTIONS¶
qng_meat¶
Signature: (double a(); double b(); double epsabs();
double epsrel(); double [o] result(); double [o] abserr();
int [o] neval(); int [o] ierr(); int warn(); SV* funcion)
info not available
qng_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qag_meat¶
Signature: (double a(); double b(); double epsabs();double epsrel(); int limit();
int key(); double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qag_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qags_meat¶
Signature: (double a(); double b(); double epsabs();double epsrel(); int limit();
double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qags_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qagp_meat¶
Signature: (double pts(l); double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qagp_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qagi_meat¶
Signature: (double epsabs();double epsrel(); int limit();
double [o] result(); double [o] abserr(); int n(); int [o] ierr();int warn();; SV* funcion)
info not available
qagi_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qagiu_meat¶
Signature: (double a(); double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qagiu_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qagil_meat¶
Signature: (double b(); double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qagil_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qawc_meat¶
Signature: (double a(); double b(); double c(); double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();int [o] ierr();int warn();; SV* funcion)
info not available
qawc_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qaws_meat¶
Signature: (double a(); double b();double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();
double alpha(); double beta(); int mu(); int nu();int [o] ierr();int warn();; SV* funcion)
info not available
qaws_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qawo_meat¶
Signature: (double a(); double b();double epsabs();double epsrel();int limit();
double [o] result(); double [o] abserr();int n();
int sincosopt(); double omega(); double L(); int nlevels();int [o] ierr();int warn();; SV* funcion)
info not available
qawo_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
qawf_meat¶
Signature: (double a(); double epsabs();int limit();
double [o] result(); double [o] abserr();int n();
int sincosopt(); double omega(); int nlevels();int [o] ierr();int warn();; SV* funcion)
info not available
qawf_meat does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.