1INTEG(3) User Contributed Perl Documentation INTEG(3)
2
6 PDL::GSL::INTEG - PDL interface to numerical integration routines in
7 GSL
8
10 This is an interface to the numerical integration package present in
11 the GNU Scientific Library, which is an implementation of QUADPACK.
12 Functions are named gslinteg_{algorithm} where {algorithm} is the
14 QUADPACK naming convention. The available functions are:
15 gslinteg_qng: Non-adaptive Gauss-Kronrod integration
17 gslinteg_qag: Adaptive integration
18 gslinteg_qags: Adaptive integration with singularities
19 gslinteg_qagp: Adaptive integration with known singular points
20 gslinteg_qagi: Adaptive integration on infinite interval of the form
21 (-\infty,\infty)
22 gslinteg_qagiu: Adaptive integration on infinite interval of the form
23 (la,\infty)
24 gslinteg_qagil: Adaptive integration on infinite interval of the form
25 (-\infty,lb)
26 gslinteg_qawc: Adaptive integration for Cauchy principal values
27 gslinteg_qaws: Adaptive integration for singular functions
28 gslinteg_qawo: Adaptive integration for oscillatory functions
29 gslinteg_qawf: Adaptive integration for Fourier integrals
30 Each algorithm computes an approximation to the integral, I, of the
32 function f(x)w(x), where w(x) is a weight function (for general
33 integrands w(x)=1). The user provides absolute and relative error
34 bounds (epsabs,epsrel) which specify the following accuracy
35 requirement:
36 |RESULT - I| <= max(epsabs, epsrel |I|)
38 The routines will fail to converge if the error bounds are too
40 stringent, but always return the best approximation obtained up to that
41 stage
42 All functions return the result, and estimate of the absolute error and
44 an error flag (which is zero if there were no problems). You are
45 responsible for checking for any errors, no warnings are issued unless
46 the option {Warn => 'y'} is specified in which case the reason of
47 failure will be printed.
48 You can nest integrals up to 20 levels. If you find yourself in the
50 unlikely situation that you need more, you can change the value of
51 'max_nested_integrals' in the first line of the file 'FUNC.c' and
52 recompile.
53
55 Throughout this documentation we strive to use the same variables that
56 are present in the original GSL documentation (see See Also).
57 Oftentimes those variables are called "a" and "b". Since good Perl
58 coding practices discourage the use of Perl variables $a and $b, here
59 we refer to Parameters "a" and "b" as $pa and $pb, respectively, and
60 Limits (of domain or integration) as $la and $lb.
61 Please check the GSL documentation for more information.
63
65 use PDL;
66 use PDL::GSL::INTEG;
67 my $la = 1.2;
69 my $lb = 3.7;
70 my $epsrel = 0;
71 my $epsabs = 1e-6;
72 # Non adaptive integration
74 my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$la,$lb,$epsrel,$epsabs);
75 # Warnings on
76 my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&myf,$la,$lb,$epsrel,$epsabs,{Warn=>'y'});
77 # Adaptive integration with warnings on
79 my $limit = 1000;
80 my $key = 5;
81 my ($res,$abserr,$ierr) = gslinteg_qag(\&myf,$la,$lb,$epsrel,
82 $epsabs,$limit,$key,{Warn=>'y'});
83 sub myf{
85 my ($x) = @_;
86 return exp(-$x**2);
87 }
88
90 qng_meat
91 Signature: (double a(); double b(); double epsabs();
92 double epsrel(); double [o] result(); double [o] abserr();
93 int [o] neval(); int [o] ierr(); int gslwarn(); SV* function)
94 info not available
96 qng_meat does not process bad values. It will set the bad-value flag
98 of all output piddles if the flag is set for any of the input piddles.
99 qag_meat
101 Signature: (double a(); double b(); double epsabs();double epsrel(); int limit();
102 int key(); double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
103 info not available
105 qag_meat does not process bad values. It will set the bad-value flag
107 of all output piddles if the flag is set for any of the input piddles.
108 qags_meat
110 Signature: (double a(); double b(); double epsabs();double epsrel(); int limit();
111 double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
112 info not available
114 qags_meat does not process bad values. It will set the bad-value flag
116 of all output piddles if the flag is set for any of the input piddles.
117 qagp_meat
119 Signature: (double pts(l); double epsabs();double epsrel();int limit();
120 double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
121 info not available
123 qagp_meat does not process bad values. It will set the bad-value flag
125 of all output piddles if the flag is set for any of the input piddles.
126 qagi_meat
128 Signature: (double epsabs();double epsrel(); int limit();
129 double [o] result(); double [o] abserr(); int n(); int [o] ierr();int gslwarn();; SV* function)
130 info not available
132 qagi_meat does not process bad values. It will set the bad-value flag
134 of all output piddles if the flag is set for any of the input piddles.
135 qagiu_meat
137 Signature: (double a(); double epsabs();double epsrel();int limit();
138 double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
139 info not available
141 qagiu_meat does not process bad values. It will set the bad-value flag
143 of all output piddles if the flag is set for any of the input piddles.
144 qagil_meat
146 Signature: (double b(); double epsabs();double epsrel();int limit();
147 double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
148 info not available
150 qagil_meat does not process bad values. It will set the bad-value flag
152 of all output piddles if the flag is set for any of the input piddles.
153 qawc_meat
155 Signature: (double a(); double b(); double c(); double epsabs();double epsrel();int limit();
156 double [o] result(); double [o] abserr();int n();int [o] ierr();int gslwarn();; SV* function)
157 info not available
159 qawc_meat does not process bad values. It will set the bad-value flag
161 of all output piddles if the flag is set for any of the input piddles.
162 qaws_meat
164 Signature: (double a(); double b();double epsabs();double epsrel();int limit();
165 double [o] result(); double [o] abserr();int n();
166 double alpha(); double beta(); int mu(); int nu();int [o] ierr();int gslwarn();; SV* function)
167 info not available
169 qaws_meat does not process bad values. It will set the bad-value flag
171 of all output piddles if the flag is set for any of the input piddles.
172 qawo_meat
174 Signature: (double a(); double b();double epsabs();double epsrel();int limit();
175 double [o] result(); double [o] abserr();int n();
176 int sincosopt(); double omega(); double L(); int nlevels();int [o] ierr();int gslwarn();; SV* function)
177 info not available
179 qawo_meat does not process bad values. It will set the bad-value flag
181 of all output piddles if the flag is set for any of the input piddles.
182 qawf_meat
184 Signature: (double a(); double epsabs();int limit();
185 double [o] result(); double [o] abserr();int n();
186 int sincosopt(); double omega(); int nlevels();int [o] ierr();int gslwarn();; SV* function)
187 info not available
189 qawf_meat does not process bad values. It will set the bad-value flag
191 of all output piddles if the flag is set for any of the input piddles.
192 gslinteg_qng
194 Non-adaptive Gauss-Kronrod integration
195 This function applies the Gauss-Kronrod 10-point, 21-point, 43-point
197 and 87-point integration rules in succession until an estimate of the
198 integral of f over ($la,$lb) is achieved within the desired absolute
199 and relative error limits, $epsabs and $epsrel. It is meant for fast
200 integration of smooth functions. It returns an array with the result,
201 an estimate of the absolute error, an error flag and the number of
202 function evaluations performed.
203 Usage:
205 ($res,$abserr,$ierr,$neval) = gslinteg_qng($function_ref,$la,$lb,
207 $epsrel,$epsabs,[{Warn => $warn}]);
208 Example:
210 my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9);
212 # with warnings on
213 my ($res,$abserr,$ierr,$neval) = gslinteg_qng(\&f,0,1,0,1e-9,{Warn => 'y'});
214 sub f{
216 my ($x) = @_;
217 return ($x**2.6)*log(1.0/$x);
218 }
219 gslinteg_qag
221 Adaptive integration
222 This function applies an integration rule adaptively until an estimate
224 of the integral of f over ($la,$lb) is achieved within the desired
225 absolute and relative error limits, $epsabs and $epsrel. On each
226 iteration the adaptive integration strategy bisects the interval with
227 the largest error estimate; the maximum number of allowed subdivisions
228 is given by the parameter $limit. The integration rule is determined
229 by the value of $key, which has to be one of (1,2,3,4,5,6) and
230 correspond to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod rules
231 respectively. It returns an array with the result, an estimate of the
232 absolute error and an error flag.
233 Please check the GSL documentation for more information.
235 Usage:
237 ($res,$abserr,$ierr) = gslinteg_qag($function_ref,$la,$lb,$epsrel,
239 $epsabs,$limit,$key,[{Warn => $warn}]);
240 Example:
242 my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1);
244 # with warnings on
245 my ($res,$abserr,$ierr) = gslinteg_qag(\&f,0,1,0,1e-10,1000,1,{Warn => 'y'});
246 sub f{
248 my ($x) = @_;
249 return ($x**2.6)*log(1.0/$x);
250 }
251 gslinteg_qags
253 Adaptive integration with singularities
254 This function applies the Gauss-Kronrod 21-point integration rule
256 adaptively until an estimate of the integral of f over ($la,$lb) is
257 achieved within the desired absolute and relative error limits, $epsabs
258 and $epsrel. The algorithm is such that it accelerates the convergence
259 of the integral in the presence of discontinuities and integrable
260 singularities. The maximum number of allowed subdivisions done by the
261 adaptive algorithm must be supplied in the parameter $limit.
262 Please check the GSL documentation for more information.
264 Usage:
266 ($res,$abserr,$ierr) = gslinteg_qags($function_ref,$la,$lb,$epsrel,
268 $epsabs,$limit,[{Warn => $warn}]);
269 Example:
271 my ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000);
273 # with warnings on
274 ($res,$abserr,$ierr) = gslinteg_qags(\&f,0,1,0,1e-10,1000,{Warn => 'y'});
275 sub f{
277 my ($x) = @_;
278 return ($x)*log(1.0/$x);
279 }
280 gslinteg_qagp
282 Adaptive integration with known singular points
283 This function applies the adaptive integration algorithm used by
285 gslinteg_qags taking into account the location of singular points until
286 an estimate of the integral of f over ($la,$lb) is achieved within the
287 desired absolute and relative error limits, $epsabs and $epsrel.
288 Singular points are supplied in the piddle $points, whose endpoints
289 determine the integration range. So, for example, if the function has
290 singular points at x_1 and x_2 and the integral is desired from a to b
291 (a < x_1 < x_2 < b), $points = pdl(a,x_1,x_2,b). The maximum number of
292 allowed subdivisions done by the adaptive algorithm must be supplied in
293 the parameter $limit.
294 Please check the GSL documentation for more information.
296 Usage:
298 ($res,$abserr,$ierr) = gslinteg_qagp($function_ref,$points,$epsabs,
300 $epsrel,$limit,[{Warn => $warn}])
301 Example:
303 my $points = pdl(0,1,sqrt(2),3);
305 my ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000);
306 # with warnings on
307 ($res,$abserr,$ierr) = gslinteg_qagp(\&f,$points,0,1e-3,1000,{Warn => 'y'});
308 sub f{
310 my ($x) = @_;
311 my $x2 = $x**2;
312 my $x3 = $x**3;
313 return $x3 * log(abs(($x2-1.0)*($x2-2.0)));
314 }
315 gslinteg_qagi
317 Adaptive integration on infinite interval
318 This function estimates the integral of the function f over the
320 infinite interval (-\infty,+\infty) within the desired absolute and
321 relative error limits, $epsabs and $epsrel. After a transformation,
322 the algorithm of gslinteg_qags with a 15-point Gauss-Kronrod rule is
323 used. The maximum number of allowed subdivisions done by the adaptive
324 algorithm must be supplied in the parameter $limit.
325 Please check the GSL documentation for more information.
327 Usage:
329 ($res,$abserr,$ierr) = gslinteg_qagi($function_ref,$epsabs,
331 $epsrel,$limit,[{Warn => $warn}]);
332 Example:
334 my ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000);
336 # with warnings on
337 ($res,$abserr,$ierr) = gslinteg_qagi(\&myfn,1e-7,0,1000,{Warn => 'y'});
338 sub myfn{
340 my ($x) = @_;
341 return exp(-$x - $x*$x) ;
342 }
343 gslinteg_qagiu
345 Adaptive integration on infinite interval
346 This function estimates the integral of the function f over the
348 infinite interval (la,+\infty) within the desired absolute and relative
349 error limits, $epsabs and $epsrel. After a transformation, the
350 algorithm of gslinteg_qags with a 15-point Gauss-Kronrod rule is used.
351 The maximum number of allowed subdivisions done by the adaptive
352 algorithm must be supplied in the parameter $limit.
353 Please check the GSL documentation for more information.
355 Usage:
357 ($res,$abserr,$ierr) = gslinteg_qagiu($function_ref,$la,$epsabs,
359 $epsrel,$limit,[{Warn => $warn}]);
360 Example:
362 my $alfa = 1;
364 my ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000);
365 # with warnings on
366 ($res,$abserr,$ierr) = gslinteg_qagiu(\&f,99.9,1e-7,0,1000,{Warn => 'y'});
367 sub f{
369 my ($x) = @_;
370 if (($x==0) && ($alfa == 1)) {return 1;}
371 if (($x==0) && ($alfa > 1)) {return 0;}
372 return ($x**($alfa-1))/((1+10*$x)**2);
373 }
374 gslinteg_qagil
376 Adaptive integration on infinite interval
377 This function estimates the integral of the function f over the
379 infinite interval (-\infty,lb) within the desired absolute and relative
380 error limits, $epsabs and $epsrel. After a transformation, the
381 algorithm of gslinteg_qags with a 15-point Gauss-Kronrod rule is used.
382 The maximum number of allowed subdivisions done by the adaptive
383 algorithm must be supplied in the parameter $limit.
384 Please check the GSL documentation for more information.
386 Usage:
388 ($res,$abserr,$ierr) = gslinteg_qagl($function_ref,$lb,$epsabs,
390 $epsrel,$limit,[{Warn => $warn}]);
391 Example:
393 my ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000);
395 # with warnings on
396 ($res,$abserr,$ierr) = gslinteg_qagil(\&myfn,1.0,1e-7,0,1000,{Warn => 'y'});
397 sub myfn{
399 my ($x) = @_;
400 return exp($x);
401 }
402 gslinteg_qawc
404 Adaptive integration for Cauchy principal values
405 This function computes the Cauchy principal value of the integral of f
407 over (la,lb), with a singularity at c, I = \int_{la}^{lb} dx f(x)/(x -
408 c). The integral is estimated within the desired absolute and relative
409 error limits, $epsabs and $epsrel. The maximum number of allowed
410 subdivisions done by the adaptive algorithm must be supplied in the
411 parameter $limit.
412 Please check the GSL documentation for more information.
414 Usage:
416 ($res,$abserr,$ierr) = gslinteg_qawc($function_ref,$la,$lb,$c,$epsabs,$epsrel,$limit)
418 Example:
420 my ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000);
422 # with warnings on
423 ($res,$abserr,$ierr) = gslinteg_qawc(\&f,-1,5,0,0,1e-3,1000,{Warn => 'y'});
424 sub f{
426 my ($x) = @_;
427 return 1.0 / (5.0 * $x * $x * $x + 6.0) ;
428 }
429 gslinteg_qaws
431 Adaptive integration for singular functions
432 The algorithm in gslinteg_qaws is designed for integrands with
434 algebraic-logarithmic singularities at the end-points of an integration
435 region. Specifically, this function computes the integral given by I =
436 \int_{la}^{lb} dx f(x) (x-la)^alpha (lb-x)^beta log^mu (x-la) log^nu
437 (lb-x). The integral is estimated within the desired absolute and
438 relative error limits, $epsabs and $epsrel. The maximum number of
439 allowed subdivisions done by the adaptive algorithm must be supplied in
440 the parameter $limit.
441 Please check the GSL documentation for more information.
443 Usage:
445 ($res,$abserr,$ierr) =
447 gslinteg_qawc($function_ref,$alpha,$beta,$mu,$nu,$la,$lb,
448 $epsabs,$epsrel,$limit,[{Warn => $warn}]);
449 Example:
451 my ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000);
453 # with warnings on
454 ($res,$abserr,$ierr) = gslinteg_qaws(\&f,0,0,1,0,0,1,0,1e-7,1000,{Warn => 'y'});
455 sub f{
457 my ($x) = @_;
458 if($x==0){return 0;}
459 else{
460 my $u = log($x);
461 my $v = 1 + $u*$u;
462 return 1.0/($v*$v);
463 }
464 }
465 gslinteg_qawo
467 Adaptive integration for oscillatory functions
468 This function uses an adaptive algorithm to compute the integral of f
470 over (la,lb) with the weight function sin(omega*x) or cos(omega*x) --
471 which of sine or cosine is used is determined by the parameter $opt
472 ('cos' or 'sin'). The integral is estimated within the desired
473 absolute and relative error limits, $epsabs and $epsrel. The maximum
474 number of allowed subdivisions done by the adaptive algorithm must be
475 supplied in the parameter $limit.
476 Please check the GSL documentation for more information.
478 Usage:
480 ($res,$abserr,$ierr) = gslinteg_qawo($function_ref,$omega,$sin_or_cos,
482 $la,$lb,$epsabs,$epsrel,$limit,[opt])
483 Example:
485 my $PI = 3.14159265358979323846264338328;
487 my ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000);
488 # with warnings on
489 ($res,$abserr,$ierr) = PDL::GSL::INTEG::gslinteg_qawo(\&f,10*$PI,'sin',0,1,0,1e-7,1000,{Warn => 'y'});
490 sub f{
492 my ($x) = @_;
493 if($x==0){return 0;}
494 else{ return log($x);}
495 }
496 gslinteg_qawf
498 Adaptive integration for Fourier integrals
499 This function attempts to compute a Fourier integral of the function f
501 over the semi-infinite interval [la,+\infty). Specifically, it attempts
502 tp compute I = \int_{la}^{+\infty} dx f(x)w(x), where w(x) is
503 sin(omega*x) or cos(omega*x) -- which of sine or cosine is used is
504 determined by the parameter $opt ('cos' or 'sin'). The integral is
505 estimated within the desired absolute error limit $epsabs. The maximum
506 number of allowed subdivisions done by the adaptive algorithm must be
507 supplied in the parameter $limit.
508 Please check the GSL documentation for more information.
510 Usage:
512 gslinteg_qawf($function_ref,$omega,$sin_or_cos,$la,$epsabs,$limit,[opt])
514 Example:
516 my ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000);
518 # with warnings on
519 ($res,$abserr,$ierr) = gslinteg_qawf(\&f,$PI/2.0,'cos',0,1e-7,1000,{Warn => 'y'});
520 sub f{
522 my ($x) = @_;
523 if ($x == 0){return 0;}
524 return 1.0/sqrt($x)
525 }
526
528 Feedback is welcome. Log bugs in the PDL bug database (the database is
529 always linked from <http://pdl.perl.org>).
530
532 PDL
533 The GSL documentation for numerical integration is online at
535 <https://www.gnu.org/software/gsl/doc/html/integration.html>
536
538 This file copyright (C) 2003,2005 Andres Jordan <ajordan@eso.org> All
539 rights reserved. There is no warranty. You are allowed to redistribute
540 this software documentation under certain conditions. For details, see
541 the file COPYING in the PDL distribution. If this file is separated
542 from the PDL distribution, the copyright notice should be included in
543 the file.
544 The GSL integration routines were written by Brian Gough. QUADPACK was
546 written by Piessens, Doncker-Kapenga, Uberhuber and Kahaner.
547perl v5.32.0 2020-09-17 INTEG(3)