I have forgotten

• https://me.yahoo.com

Adaptive integration algorithm using RMS formulae
Controller: CodeCogs

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:
• 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.

The error bounds for RMS formulae are shown to be better than those for Newton-Cotes and comparable to the Clenshaw-Curtis, Gauss and Gauss-Kronrod formulae. Also, as opposed to Gauss and Gauss-Kronrod rules, RMS formulae do not miss the previously computed integral estimates.

The adaptive strategy used by this algorithm is to extrapolate the function values using the algorithm. The purpose of this algorithm is to eliminate the effects of the various kinds of singularities in the function to integrate.

This module is based on the 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

// 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

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

Integrate

 static doubleintegrate( double (*f)(double)[function pointer] double a double b )
This function estimates the definite integral using an adaptive integration algorithm based on RMS formulae.
 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

 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 )
This method computes with error estimate and conditionally computes by using a RMS rule.
 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 integral ym if ke = k1, then this will store an approximation to the integral of over ke 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 = 5 k1 value of key for which the additional estimates ya and ym are to be computed fv1 vector containing l1 saved functional values of positive abscissas fv2 vector containing l1 saved functional values of positive abscissas l1 number of elements in fv1 l2 number of elements in fv2

Dqxlqm

 static 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 )
This method implements the local quadrature module (LQM) used when estimating the integral. Instead of using the pair 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 abscissas vs vector of length ls containing the saved functional values of positive abscissas lr number of elements in vr ls number of elements in vs

Dqxrrd

 static voiddqxrrd( double (*f)(double)[function pointer] double* z int lz double& xl double& xu double* r double* s int& lr int& ls )
This method reorders the computed functional values before the bisection of an interval.
 f function defining the integrand z vector containing lz saved functional values lz 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 intervals s vector containing ls saved functional values for the new intervals lr number of elements in r ls number of elements in s

Dqpsrt

 static voiddqpsrt( int last int& maxerr double& ermax double* elist int* iord int nrmax )
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.
 last number of error estimates currently in the list maxerr points to the nrmax -th largest error estimate currently in the list ermax nrmax -th largest error estimate; ermax = elist [maxerr] elist vector of dimension last containing the error estimates iord 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 otherwise nrmax value with the property that maxerr = iord [nrmax]

Dqelg

 static voiddqelg( int& n double* epstab double& result double& abserr double* res3la int& nres )
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.
 n epstab [n] contains the new element in the first column of the table epstab 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 results res3la vector of dimension 3 containing the last three results nres number of calls to the routine (should be zero at first call)

Dqxgse

 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 )
This method calculates an approximation to the definite integral hopefully satisfying the following claim for accuracy:
 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 subintervals elist vector of dimension at least DIV_LIMIT, the first last elements of which are the moduli of the absolute error estimates on the subintervals iord 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 otherwise last number of subintervals actually produced in the subdivision process valp array of dimension at least (21, DIV_LIMIT) used to save the functional values valn array of dimension at least (21, DIV_LIMIT) used to save the functional values lp 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