File: Numerical-integration-examples.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (146 lines) | stat: -rw-r--r-- 5,938 bytes parent folder | download
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146  GNU Scientific Library – Reference Manual: Numerical integration examples

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.