Numerical Integration

This chapter describes routines for performing numerical integration (quadrature) of a function in one dimension. There are routines for adaptive and non-adaptive integration of general functions, with specialised routines for specific cases. These include integration over infinite and semi-infinite ranges, singular integrals, including logarithmic singularities, computation of Cauchy principal values and oscillatory integrals. The library reimplements the algorithms used in QUADPACK, a numerical integration package written by Piessens, de Doncker-Kapenga, Ueberhuber and Kahaner. Fortran code for QUADPACK is available on Netlib. Also included are non-adaptive, fixed-order Gauss-Legendre integration routines with high precision coefficients by Pavel Holoborodko.

The functions described in this chapter are declared in the header file gsl_integration.h.

• Numerical Integration Introduction:
• QNG non-adaptive Gauss-Kronrod integration:
• QAG adaptive integration:
• QAGS adaptive integration with singularities:
• QAGP adaptive integration with known singular points:
• QAGI adaptive integration on infinite intervals:
• QAWC adaptive integration for Cauchy principal values:
• QAWS adaptive integration for singular functions:
• QAWO adaptive integration for oscillatory functions:
• QAWF adaptive integration for Fourier integrals:
• CQUAD doubly-adaptive integration:
• Fixed order Gauss-Legendre integration:
• Numerical integration error codes:
• Numerical integration examples:
• Numerical integration References and Further Reading:

Building GSL on Windows

3.22.2014

A few pointers for compiling Windows binaries of the GNU Scientific Library (GSL).

Description

The GNU Scientific Library (GSL) is a fantastic library for scientific computing. Within the C/C++ world I don't think there's anything comparable. My only gripe is the C-style interface, I would prefer C++. Anyways, I tend to use it quite a bit. Personally I use a combination of OS X and Linux, both of which are easy platforms for GSL - the latter moreso, but on OS X you can get it easily through the package managers (I currently use brew).

Windows

Unfortunately many other people in this world still use Windows for scientific computing, so I have to be cross-platform. I've starting using the Visual Studio toolchain, in particular Visual Studio 2013. I can get it for free as a student (through Dreamspark) and it plays nicely with Python.

After some internet sleuthing I found that someone has already made a Visual Studio project for compiling GSL, which can be foundhere. Thanks Brian! He only advertises it as working with Visual Studio 2012. However with a few tweaks I managed to get it to work on 2013. There were a few header files the project looked for in the wrong place. Also it is finicky about order. I found out that you have to build the statically-linked libraries first, and then the dynamically-linked libraries. For the latter I had to copy the *.lib files into wherever it wanted them. Fortunately the error messages encountered tend to be helpful.

As a simplicity for anyone who wants them, I built binaries (lib and dll) using GSL v1.16 and VS2013 for both 32-bit and 64-bit architecture, and they are available here (zip, 12MB). Have fun!

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 =  %d\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          = -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.