Click here to Skip to main content
15,867,568 members
Articles / Programming Languages / C++

Fast Numerical Integration

Rate me:
Please Sign up or sign in to vote.
4.89/5 (66 votes)
28 Jun 2014BSD7 min read 161.3K   2.7K   100   49
Numerical integration of smooth functions over a finite interval using an optimal algorithm.

Introduction

Most integrals that come up in real applications have to be evaluated numerically. Unfortunately, there is no way to write code that will efficiently and accurately evaluate any integral you throw at it. However, it is possible to write integration routines that work well on particular classes of problems.

The method presented here integrates smooth functions over finite intervals. In a sense, the method is optimally efficient, achieving the most accuracy for a given number of function evaluations.

To be more precise, the method presented here assumes that the functions being integrated are analytic. This means the functions are required to be very smooth, no discontinuities in the function or any of its derivatives. It rules out functions patched together using "if" statements, such as "return x3 if x is less than 0, and return x2 if x is positive". It also rules out functions that become infinite in the middle of the integration interval. On the other hand, it does allow functions that become infinite at the ends of the integration interval.

The method presented here is the double exponential transformation. For an explanation of the mathematical details, see the article "The double-exponential transformation in numerical analysis" by Masatake Mori and Masaaki Sugihara in the Journal of Computational and Applied Mathematics, volume 127 (2001), pages 287-296. The mathematics behind the method is quite sophisticated, and will not be presented here. However, the code that implements the method is simple and easy to use.

Using the Code

The integrator is implemented in a templated class DEIntegrator. (The "DE" in the name stands for the "double exponential" algorithm at the heart of the code.) The template parameter is the class of the function to be integrated. The code uses function classes rather than function pointers because the former are more flexible. Functions often have parameters that tag along, and function classes are the way to handle that in C++. (In functional languages, you might use closures.)

To use DEIntegrator in your code, you need to add two files to your project, DEIntegrator.h and DEIntegratorConstants.h, and you need to add #include "DEIntegrator.h" in the file where you want to call the integrator.

DEIntegrator has two overloaded versions of the Integrate method. The first is the easiest to use:

C++
static double Integrate
(
    const TFunctionObject& f,       // [in] integrand
    double a,                       // [in] left limit of integration
    double b,                       // [in] right limit of integration
    double targetAbsoluteError      // [in] desired bound on error
)

You simply pass in the function to integrate, the interval you want to integrate over, and how accurately you want the result calculated. The return value is the value of the integral.

The other overloaded version allows you to get more information after the integral is calculated.

C++
static double Integrate
(
    const TFunctionObject& f,       // [in] integrand
    double a,                       // [in] left limit of integration
    double b,                       // [in] right limit of integration
    double targetAbsoluteError,     // [in] desired bound on error
    int& numFunctionEvaluations,    // [out] number of function evaluations used
    double& errorEstimate           // [out] estimated error in integration
)

The two additional arguments let you know how many function evaluations were required in computing the integral as well as the algorithm's estimate of the error. Note that targetAbsoluteError specifies how small you want the error to be, but errorEstimate estimates how small the error actually is. Often, the latter will be much smaller than the former, i.e., the method will do a better job than you asked for. However, if the method has difficulty with your integrand (say, due to a strong singularity at one of the end points), then errorEstimate may be larger than targetAbsoluteError, letting you know there is a problem.

Example 1

Suppose we want to integrate f(x) = exp(-x/5) (2 + sin(2x)) from 0 to 10, plotted below:

Plot of exp(-x/5) (2 + sin(2x))

It turns out that this integral can be calculated exactly. It equals 9.1082396073230. This is a somewhat artificial example: usually, you use numerical integration for integrals that cannot be calculated exactly. But, exact integrals are useful for illustrations and testing. Suppose we want to numerically evaluate the integral to six decimal places.

First, we write our function object:

C++
class DemoFunction
{
public:
    double operator()(double x) const
    {
        return exp(-x/5.0)*(2.0 + sin(2.0*x));
    }
};

Then, we evaluate its integral as follows:

C++
DemoFunction f;
int evaluations;
double errorEstimate;
double integral = DEIntegrator<DemoFunction>::Integrate
        (f, 0, 10, 1e-6, evaluations, errorEstimate);
std::cout << integral << ", " << errorEstimate 
          << ", " << evaluations << "\n";

The output is:

9.10823960732284, 5.21017118337852e-009, 97

And, so in this example, we requested that the error be less than 10-6, and the method estimates that the error was less than 5 × 10-9, and in fact, the error was more like 10-13. The method under-promised and over-delivered. It used 97 function evaluations in the process.

Example 2

As another example, we look at f(x) = x-1/3(1-x)5, plotted below. We want to integrate this function over the interval from 0 to 1. This example is to show that the integrand does not need to be analytic everywhere, just on the interior of the integration interval. In other words, it's OK if the function blows up at one or both of the ends.

Plot of x^{-1/3>(1-x)^5

We implement a function object for our integrand as follows:

C++
class DemoFunction2
{
public:
    double operator()(double x) const
    {
        return pow(1.0 - x, 5.0)*pow(x, -1.0/3.0);
    }
};

Now, suppose we only want the estimated integral. We're not interested in counting how many function evaluations were required or the estimated error. In this case, the code to evaluate the integral is more compact.

C++
DemoFunction2 f2;
integral = DEIntegrator<DemoFunction2>::Integrate(f2, 0, 1, 1e-6);
std::cout << integral << "\n";

In this case, the integral also happens to be one that can be evaluated exactly. In fact, the integral equals 2187/5236. The software computes the integral as 0.41768525570112. The error is approximately 2 × 10-10.

Error Estimates

The error in elementary numerical integrations goes down like some constant power of the number of integration points. For example, the error in Simpson's rule decreases like N-4 where N is the number of integration points. But, the error in the method presented here decreases exponentially. Specifically, the error is on the order of exp(-cN/log(N)). See this page for the mathematical details.

Internally, the code monitors its own error rate by comparing the differences in consecutive iterations to those the theory would predict. This allows the user to specify an error tolerance, and it allows the code to estimate after the fact how well it did.

Discussion

The method discussed here accurately and efficiently integrates analytic functions over a finite interval. The method is as efficient as possible for the class of functions it integrates, in a technical sense described here.

Unlike most methods presented in college courses, the method presented here does not assume that the function being integrated behaves roughly like a polynomial. Integrands may be singular at an end point, as in Example 2, without harming the accuracy of the method. Also, unlike more elementary methods, the convergence rates are exponential rather than polynomial. This means the method presented here converges very quickly.

There are variations of the double exponential rule for integrals over unbounded regions. Unfortunately, these are difficult to implement in general; the software needs to know more details about the function being integrated. A simpler approach would be to transform the integration problem into one that can be computed with the software presented here. The most obvious approach would be to truncate the unbounded integral to a bounded integral. Another approach would be to use a change of variables to transform the integral into a new integral over a bounded interval.

Truncating an integral over an unbounded region is not recommended. It can be difficult to decide where to truncate. Picking too small an interval can lead to an inaccurate answer. Picking too large an interval can lead to wasting time integrating over regions where the integrand is essentially zero.

Using a change of variables to transform the integral may be more efficient and accurate. For integrals over (0, ∞), a common change of variables is x = t/(1-t). This transforms the integral of f(x) from 0 to ∞ into the integral of f(t/(1-t)) /(1-t)2 from 0 to 1. For integrals over (-∞, ∞), variables x = tan(t) transforms the integral of f(x) from -∞ to ∞ into the integral of f(tan(t)) (1 + sec2(t)) from -1 to 1. These changes of variables often work well. It is possible to get better performance by crafting a transformation specifically for your integration problem, but this is usually not necessary.

History

  • 5th December, 2008: Initial version.
  • 8th December, 2008: Added discussion of error estimates.
  • 29th June, 2009: Added discussion of integrals over unbounded regions.
  • 4th May, 2015: Fixed bug in handling target absolute error.

License

This article, along with any associated source code and files, is licensed under The BSD License


Written By
President John D. Cook Consulting
United States United States
I work in the areas of applied mathematics, data analysis, and data privacy.

Check out my blog or send me a note.

 


Comments and Discussions

 
AnswerRe: File invalid Pin
John D. Cook28-Jun-14 10:05
John D. Cook28-Jun-14 10:05 
GeneralRe: File invalid Pin
David Luca29-Jun-14 21:54
professionalDavid Luca29-Jun-14 21:54 
GeneralThank you! Pin
JoraSteelRat1-Oct-13 0:04
JoraSteelRat1-Oct-13 0:04 
GeneralMy vote of 5 Pin
Koh Twee Toong29-Jun-13 7:39
Koh Twee Toong29-Jun-13 7:39 
QuestionGreat code (and a quick question) Pin
Member 992827620-Mar-13 10:43
Member 992827620-Mar-13 10:43 
GeneralMy vote of 5 Pin
themainsequence8-Oct-12 22:53
themainsequence8-Oct-12 22:53 
GeneralMy vote of 5 Pin
burndavid9327-May-12 5:20
burndavid9327-May-12 5:20 
GeneralMy vote of 5 Pin
MattEFowler12-Apr-12 10:27
MattEFowler12-Apr-12 10:27 
Questiona function to test Pin
dariusz_brzezinski22-Mar-12 23:23
dariusz_brzezinski22-Mar-12 23:23 
AnswerRe: a function to test Pin
John D. Cook23-Mar-12 0:10
John D. Cook23-Mar-12 0:10 
GeneralReally good work [modified] Pin
V. Thieme8-Nov-10 10:27
V. Thieme8-Nov-10 10:27 
Generalvery good, thank you. Pin
jemisha28-Jul-10 6:27
jemisha28-Jul-10 6:27 
GeneralKrylov method Pin
Majid Shahabfar6-Jul-09 22:06
Majid Shahabfar6-Jul-09 22:06 
GeneralRe: Krylov method Pin
John D. Cook7-Jul-09 0:58
John D. Cook7-Jul-09 0:58 
QuestionMultiple integration routine? Pin
Graeme Addison18-Feb-09 11:30
Graeme Addison18-Feb-09 11:30 
AnswerRe: Multiple integration routine? Pin
John D. Cook18-Feb-09 15:52
John D. Cook18-Feb-09 15:52 
GeneralVery cool! Pin
Dr.Luiji16-Jan-09 22:52
professionalDr.Luiji16-Jan-09 22:52 
GeneralRe: Very cool! Pin
John D. Cook17-Jan-09 3:00
John D. Cook17-Jan-09 3:00 
GeneralWhy function objects Pin
John D. Cook19-Dec-08 6:12
John D. Cook19-Dec-08 6:12 
Generalcode also available for g++ Pin
andreas19808-Dec-08 22:26
andreas19808-Dec-08 22:26 
GeneralRe: code also available for g++ Pin
John D. Cook8-Dec-08 23:51
John D. Cook8-Dec-08 23:51 
GeneralRe: code also available for g++ Pin
andreas198010-Dec-08 5:56
andreas198010-Dec-08 5:56 
GeneralAccuracy Pin
Member 27859438-Dec-08 14:10
Member 27859438-Dec-08 14:10 
GeneralRe: Accuracy Pin
John D. Cook8-Dec-08 15:38
John D. Cook8-Dec-08 15:38 
GeneralMathematical details posted [modified] Pin
John D. Cook8-Dec-08 3:45
John D. Cook8-Dec-08 3:45 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.