The integrator QAGS
will handle a large class of definite
integrals. For example, consider the following integral, which has a
algebraic-logarithmic singularity at the origin,
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 = α 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 = %d\n", w->size); return 0; }
The results below show that the desired accuracy is achieved after 8 subdivisions.
bash$ ./a.out result = -3.999999999999973799 exact result = -4.000000000000000000 estimated error = 0.000000000000246025 actual error = 0.000000000000026201 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.