17.14 Examples

The integrator QAGS will handle a large class of definite integrals. For example, consider the following integral, which has an algebraic-logarithmic singularity at the origin,

\int_0^1 x^{-1/2} log(x) dx = -4


The program below computes this integral to a relative accuracy bound of 1e-7.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>

double f (double x, void * params) {
double alpha = *(double *) params;
double f = log(alpha*x) / sqrt(x);
return f;
}

int
main (void)
{
gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);

double result, error;
double expected = -4.0;
double alpha = 1.0;

gsl_function F;
F.function = &f;
F.params = &alpha;

gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
w, &result, &error);

printf ("result          = % .18f\n", result);
printf ("exact result    = % .18f\n", expected);
printf ("estimated error = % .18f\n", error);
printf ("actual error    = % .18f\n", result - expected);
printf ("intervals       = %zu\n", w->size);

gsl_integration_workspace_free (w);

return 0;
}


The results below show that the desired accuracy is achieved after 8 subdivisions.

\$ ./a.out

result          = -4.000000000000085265
exact result    = -4.000000000000000000
estimated error =  0.000000000000135447
actual error    = -0.000000000000085265
intervals       = 8


In fact, the extrapolation procedure used by QAGS produces an accuracy of almost twice as many digits. The error estimate returned by the extrapolation procedure is larger than the actual error, giving a margin of safety of one order of magnitude.