Integration of 1D real functions

The class and methods defined here allow the integration of a real function. Each user-defined function to be integrated should be defined as a method in a class which inherits from the class defining the integration methods.

Integrator class, Real_integrable

1. Define a subclass which inherits from this class for the particular function to be integrated.
2. Define a method in the class to calculate the function:

      double operator()(const double x) const
      

3. Call the appropriate integration method:
   • Trapezium rule with accuracy checking: Refines number of steps by doubling until change between successive tries is within user - specified accuracy limit. Based on idea of calculating only intermediate points, so don't make extra calls to function routine. Use actual accuracy rather than fractional for robustness when f(x) = 0 Minimum number of steps supplied to help routine if function oscillatory. double trapezium( const double x1, const double x2, // lower and upper limits of integration const double acc = 1.0e-6, // accuracy required const int maxrep = 10, // max times allowed to half step (0 => no limit) const int minsteps = 10, // min number of steps expected to be necessary const double eps = 1.0e-10 // numerical bandwidth used to detect stepsize underflow )
   • Simpson's rule with accuracy checking: Refines number of steps by doubling until change between successive tries is within user - specified accuracy limit. Based on idea of calculating only intermediate points, so don't make extra calls to function routine, but with more aesthetic initialisation. Use actual accuracy rather than fractional for robustness when f(x) = 0 Minimum number of steps supplied to help routine if function oscillatory.

        double simpson(
           const double x1, const double x2, // lower and upper limits of integration
           const double acc = 1.0e-6, // accuracy required
           const int maxrep = 10, // max times allowed to half step (0 => no limit)
           const int minsteps = 10, // min number of steps expected to be necessary
           const double eps = 1.0e-10 // numerical bandwidth used to detect stepsize underflow
        )
     
        
     Romberg integration:
        Romberg's method of arbitrary order and accuracy.
        Performs successive Trapezoidal integrations, storing the
        results (stepsize against integral) until polynomial extrapolation
        to zero stepsize (using Neville's algorithm) converges.
        double romberg(
           const double x1, const double x2, // lower and upper limits of integration
           const double acc = 1.0e-6, // accuracy required
           const int iorder = 4, // order of integration, min 4, must be even
           const int maxrep = 10, // max times allowed to half step (0 => no limit)
           const int minsteps = 10, // min number of steps expected to be necessary
           const double eps = 1.0e-10 // numerical bandwidth used to detect stepsize underflow
        )