Class inte_gauss_cern (o2scl)

O2scl : Class List

template<class func_t = funct, class fp_t = double, class weights_t = inte_gauss_coeffs_double>
class inte_gauss_cern : public o2scl::inte<funct, double>

Gaussian quadrature (CERNLIB)

For any interval \( (a,b) \), we define \( g_8(a,b) \) and \( g_{16}(a,b) \) to be the 8- and 16-point Gaussian quadrature approximations to

\[ I=\int_a^b f(x)~dx \]
and define
\[ r(a,b)=\frac{|g_{16}(a,b)-g_{8}(a,b)|}{1+g_{16}(a,b)} \]
The function integ() returns \( G \) given by
\[ G=\sum_{i=1}^{k} g_{16}(x_{i-1},x_i) \]
where \( x_0=a \) and \( x_k=b \) and the subdivision points \( x_i \) are given by
\[ x_i=x_{i-1}+\lambda(B-x_{i-1}) \]
where \( \lambda \) is the first number in the sequence \( 1,\frac{1}{2},\frac{1}{4},... \) for which
\[ r(x_{i-1},x_i)<\mathrm{eps}. \]
If, at any stage, the ratio
\[ q=\left| \frac{x_i-x_{i-1}}{b-a} \right| \]
is so small so that \( 1+0.005 q \) is indistinguishable from unity, then the accuracy is required is not reachable and the error handler is called.

Unless there is severe cancellation, inte::tol_rel may be considered as specifying a bound on the relative error of the integral in the case that \( |I|>1 \) and an absolute error if \( |I|<1 \). More precisely, if \( k \) is the number of subintervals from above, and if

\[ I_{abs} = \int_a^b |f(x)|~dx \]
then
\[ \frac{|G-I|}{I_{abs}+k}<\mathrm{tol_rel} \]
will nearly always be true when no error is returned. For functions with no singualarities in the interval, the accuracy will usually be higher than this.

If the desired accuracy is not achieved, the integration functions will call the error handler and return the best guess, unless inte::err_nonconv is false, in which case the error handler is not called.

The constructor sets the value of inte::tol_rel to the square root of numeric_limits::epsilon. For double precision numbers, this tolerance is usually about \( 10^{-8} \).

This function is based on the CERNLIB routines GAUSS and DGAUSS which are documented at http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d103/top.html . (3/10/2020: The CERNLIB links are apparently dead and haven’t been supported since 2003.)

Idea for Future:

Allow user to change cst? Is the current value appropriate for multiprecision types?

Public Functions

inline inte_gauss_cern()
inline virtual ~inte_gauss_cern()
inline virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)

Integrate function func from a to b.

Public Members

weights_t wgts

Weights object.

Protected Attributes

const fp_t *w

Pointer to the integration weights.

const fp_t *x

Pointer to the integration abscissae.