# Adaptive

Adaptive integration algorithm using RMS formulae

Controller: **CodeCogs**

## Interface

C++

## Overview

The aim of this module is to implement an algorithm that can automatically estimate the definite integral of any one-variable function over some prescribed closed interval : The algorithm should give proper estimates for the definite integral even if the function to integrate is singular at unknown positions in the integration interval . This module implements an adaptive integration algorithm using RMS formulae, which is an improved version of the algorithm used in the**quadpack**Netlib library. The algorithm is able to automatically handle bad behavior in the integrand, like singularities at unknown positions in the integration interval.

Recursive monotone stable (RMS) formulae are a set of integration formulae , , ..., to be applied sequentially, with the following properties:

This module is based on the - the nodes used by the formula are a subset of the nodes used by , for all , which means that the previously computed integral estimates are not wasted,
- composition of these formulae does not imply wasting of functional values,
- the precision obtained by using formula is higher than the one obtained by using , for all ,
- the distance between the nodes is shorter near the boundary of the integration interval,
- all formulae are generalizations of the simple trapezoidal formula,
- all the used weights are positive, i.e. all formulae are numerically stable.

**quadpack**Netlib library, but replaces its local quadrature module (LQM) from using a pair of Gauss-Kronrod formulae called

**QK21**, to using four RMS formulae with 13, 19, 27 and 41 nodes called , , and correspondingly. These are the best formulae found by the original authors, among the total of 27 families of RMS formulae. The algorithm is numbered 691 in the ACM mathematical software list and its authors are Paola Favati, Grazia Lotti and Francesco Romani. A link to their paper is provided in the reference list below. The main function in this module is

**integrate**, which takes as parameters the integrand and the limits of the integration interval, and returns the estimated value of the integral. The following example uses the

**integrate**function to estimate the integral of over the unit interval [0, 1]. Note that this integrand is singular at both end points of the unit interval.

#include <codecogs/maths/calculus/quadrature/adaptive.h> #include <iostream> #include <iomanip> // define the integrand double f(double x) { return log(x) * sqrt(x / (1 - x)); } int main() { // calculate the integral double integral = Maths::Calculus::Quadrature::Adaptive::integrate(f, 0, 1); // set the display precision to 15 digits std::cout << std::setprecision(15); // display the integral value std::cout << "Integral: " << integral << std::endl; return 0; }

**Output:**

Integral: -0.606789763508705

### References

- Algorithm 691: Improving QUADPACK automatic integration routines, http://portal.acm.org/citation.cfm?id=108556.108580
- The quadpack library, http://www.netlib.org/quadpack/
- The MathWorld page on Wynn's epsilon method, http://mathworld.wolfram.com/WynnsEpsilonMethod.html

## Class Adaptive

## Members of Adaptive

#### Integrate

This function estimates the definite integral using an adaptive integration algorithm based on RMS formulae.static doubleintegrate( double (* `f`)(double)*[function pointer]*double `a`double `b`) f the function to integrate a the lower limit of the integration interval b the upper limit of the integration interval

### Returns

- the estimate for the definite integral of
*f*over [*a*,*b*]

#### Dqxrul

This method computes with error estimate and conditionally computes by using a RMS rule.static voiddqxrul( double (* `f`)(double)*[function pointer]*double `xl`double `xu`double *&*`y`double *&*`ya`double *&*`ym`int *&*`ke`int *&*`k1`double* `fv1`double* `fv2`int *&*`l1`int *&*`l2`) f function defining the integrand xl lower limit of integration xu upper limit of integration y stores an approximation to the integral, computed by applying the requested RMS rule ya if *ke*=*k1*, then this will store an approximation to the integralym if *ke*=*k1*, then this will store an approximation to the integral of overke key for choice of local integral rule; a RMS rule is used with: 13 points if *ke*= 2, 19 points if*ke*= 3, 27 points if*ke*= 4 and 42 points if*ke*= 5k1 value of key for which the additional estimates *ya*and*ym*are to be computedfv1 vector containing *l1*saved functional values of positive abscissasfv2 vector containing *l1*saved functional values of positive abscissasl1 number of elements in *fv1*l2 number of elements in *fv2*

#### Dqxlqm

This method implements the local quadrature module (LQM) used when estimating the integral. Instead of using the pairstatic voiddqxlqm( double (* `f`)(double)*[function pointer]*double `a`double `b`double *&*`result`double *&*`abserr`double *&*`resabs`double *&*`resasc`double* `vr`double* `vs`int *&*`lr`int *&*`ls`) **QK21**of Gauss-Kronrod rules, this LQM implements four RMS formulae of 13, 19, 27 and 41 nodes called , , and respectively. The integrals computed by this method are with its error estimate and:f function defining the integrand a lower limit of integration b upper limit of integration result approximation to the integral abserr estimate of the modulus of the absolute error, which should not exceed resabs approximation to the integral resasc approximation to the integral of over vr vector of length *lr*containing the saved functional values of positive abscissasvs vector of length *ls*containing the saved functional values of positive abscissaslr number of elements in *vr*ls number of elements in *vs*

#### Dqxrrd

This method reorders the computed functional values before the bisection of an interval.static voiddqxrrd( double (* `f`)(double)*[function pointer]*double* `z`int `lz`double *&*`xl`double *&*`xu`double* `r`double* `s`int *&*`lr`int *&*`ls`) f function defining the integrand z vector containing *lz*saved functional valueslz number of elements in *z*xl lower limit of interval xu upper limit of interval r vector containing *lr*saved functional values for the new intervalss vector containing *ls*saved functional values for the new intervalslr number of elements in *r*ls number of elements in *s*

#### Dqpsrt

This method maintains the descending ordering in the list of the local error estimates resulting from the interval subdivision process. At each call, two error estimates are inserted using the sequential search method, top-down for the largest error estimate and bottom-up for the smallest error estimate.static voiddqpsrt( int `last`int *&*`maxerr`double *&*`ermax`double* `elist`int* `iord`int `nrmax`) last number of error estimates currently in the list maxerr points to the *nrmax*-th largest error estimate currently in the listermax *nrmax*-th largest error estimate;*ermax*=*elist*[*maxerr*]elist vector of dimension *last*containing the error estimatesiord vector of dimension *last*, whose first*k*elements contain pointers to the error estimates, such that*elist*[*iord*[1]], ...,*elist*[*iord*[k]] form a decreasing sequence, where*k*=*last*if*last*(*DIV_LIMIT*/ 2 + 2) and*k*=*DIV_LIMIT*+ 1 -*last*otherwisenrmax value with the property that *maxerr*=*iord*[*nrmax*]

#### Dqelg

This method determines the limit of a given sequence of approximation, by means of the algorithm of P. Wynn. An estimate of the absolute error is also given. The condensed table is computed. Only those elements needed for the computation of the next diagonal are preserved.static voiddqelg( int *&*`n`double* `epstab`double *&*`result`double *&*`abserr`double* `res3la`int *&*`nres`) n *epstab*[*n*] contains the new element in the first column of the tableepstab vector of dimension 52 containing the elements of the two lower diagonals of the triangular table; the elements are numbered starting at the right-hand corner of the triangle result resulting approximation to the integral abserr estimate of the absolute error computed from *result*and the three previous resultsres3la vector of dimension 3 containing the last three results nres number of calls to the routine (should be zero at first call)

#### Dqxgse

This method calculates an approximation to the definite integral hopefully satisfying the following claim for accuracy:static voiddqxgse( double (* `f`)(double)*[function pointer]*double `a`double `b`double *&*`result`double *&*`abserr`int *&*`ier`double* `alist`double* `blist`double* `rlist`double* `elist`int* `iord`int *&*`last`double* `valp`double* `valn`int* `lp`int* `ln`) f function defining the integrand a lower limit of integration b upper limit of integration result approximation to the integral abserr estimate of the modulus of the absolute error, which should equal or exceed | - *result|*ier if *ier*= 0, then the termination of this method was normal and reliable (it is assumed that the requested accuracy has been achieved); if*ier*> 0, there was an abnormal termination of this method, meaning that the estimates for the integral and the error are less reliable (it is assumed that the requested accuracy has not been achieved)alist vector of dimension at least *DIV_LIMIT*, the first*last*elements of which are the left end points of the subintervals in the partition of the given integration range [*a*,*b*]blist vector of dimension at least *DIV_LIMIT*, the first*last*elements of which are the right end points of the subintervals in the partition of the given integration range [*a*,*b*]rlist vector of dimension at least *DIV_LIMIT*, the first*last*elements of which are the integral approximations on the subintervalselist vector of dimension at least *DIV_LIMIT*, the first*last*elements of which are the moduli of the absolute error estimates on the subintervalsiord vector of dimension at least *DIV_LIMIT*, the first*k*elements of which are pointers to the error estimates over the subintervals, such that*elist*[*iord*[1]], ...,*elist*[*iord*[k]] form a decreasing sequence, with*k*=*last*if*last*(*DIV_LIMIT*/ 2 + 2) and*k*=*DIV_LIMIT*+ 1 -*last*otherwiselast number of subintervals actually produced in the subdivision process valp array of dimension at least (21, *DIV_LIMIT*) used to save the functional valuesvaln array of dimension at least (21, *DIV_LIMIT*) used to save the functional valueslp vectors of dimension at least *DIV_LIMIT*, used to store the actual number of functional values saved in the corresponding column of*valp*ln vectors of dimension at least *DIV_LIMIT*, used to store the actual number of functional values saved in the corresponding column of*valn*