17.9 QAWO adaptive integration for oscillatory functions

The QAWO algorithm is designed for integrands with an oscillatory factor, \sin(\omega x) or \cos(\omega x). In order to work efficiently the algorithm requires a table of Chebyshev moments which must be pre-computed with calls to the functions below.

Function: gsl_integration_qawo_table * gsl_integration_qawo_table_alloc (double omega, double L, enum gsl_integration_qawo_enum sine, size_t n)

This function allocates space for a gsl_integration_qawo_table struct and its associated workspace describing a sine or cosine weight function W(x) with the parameters (\omega, L),

W(x) = sin(omega x)
W(x) = cos(omega x)

The parameter L must be the length of the interval over which the function will be integrated L = b - a. The choice of sine or cosine is made with the parameter sine which should be chosen from one of the two following symbolic values:

GSL_INTEG_COSINE
GSL_INTEG_SINE

The gsl_integration_qawo_table is a table of the trigonometric coefficients required in the integration process. The parameter n determines the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that n levels are sufficient for subintervals down to the length L/2^n. The integration routine gsl_integration_qawo returns the error GSL_ETABLE if the number of levels is insufficient for the requested accuracy.

Function: int gsl_integration_qawo_table_set (gsl_integration_qawo_table * t, double omega, double L, enum gsl_integration_qawo_enum sine)

This function changes the parameters omega, L and sine of the existing workspace t.

Function: int gsl_integration_qawo_table_set_length (gsl_integration_qawo_table * t, double L)

This function allows the length parameter L of the workspace t to be changed.

Function: void gsl_integration_qawo_table_free (gsl_integration_qawo_table * t)

This function frees all the memory associated with the workspace t.

Function: int gsl_integration_qawo (gsl_function * f, const double a, const double epsabs, const double epsrel, const size_t limit, gsl_integration_workspace * workspace, gsl_integration_qawo_table * wf, double * result, double * abserr)

This function uses an adaptive algorithm to compute the integral of f over (a,b) with the weight function \sin(\omega x) or \cos(\omega x) defined by the table wf,

I = \int_a^b dx f(x) sin(omega x)
I = \int_a^b dx f(x) cos(omega x)

The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. The subintervals and their results are stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated size of the workspace.

Those subintervals with “large” widths d where d\omega > 4 are computed using a 25-point Clenshaw-Curtis integration rule, which handles the oscillatory behavior. Subintervals with a “small” widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration.