Click here to Skip to main content
12,750,351 members (37,118 online)
Click here to Skip to main content
Add your own
alternative version


96 bookmarked
Posted 6 Dec 2008

Fast Numerical Integration

, 28 Jun 2014 BSD
Rate this:
Please Sign up or sign in to vote.
Numerical integration of smooth functions over a finite interval using an optimal algorithm.


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:

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.

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:

class DemoFunction
    double operator()(double x) const
        return exp(-x/5.0)*(2.0 + sin(2.0*x));

Then, we evaluate its integral as follows:

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:

class DemoFunction2
    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.

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.


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.


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


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


About the Author

John D. Cook
Singular Value Consulting
United States United States
I am an independent consultant in software development and applied mathematics. I help companies learn from their data to make better decisions.

Check out my blog or send me a note.


You may also be interested in...

Comments and Discussions

QuestionWhat is the formula for the weights? Pin
Member 1287499311-Dec-16 17:15
memberMember 1287499311-Dec-16 17:15 
QuestionWhere is the source code? Pin
Member 1168303211-May-15 20:30
memberMember 1168303211-May-15 20:30 
AnswerRe: Where is the source code? Pin
John D. Cook12-May-15 2:54
memberJohn D. Cook12-May-15 2:54 
GeneralRe: Where is the source code? Pin
Member 116830326-Jun-15 15:34
memberMember 116830326-Jun-15 15:34 
QuestionMy vote of 5! Pin
Volynsky Alex12-Nov-14 5:05
professionalVolynsky Alex12-Nov-14 5:05 
GeneralMy vote of 5 Pin
CatchExAs9-Jul-14 4:55
professionalCatchExAs9-Jul-14 4:55 
QuestionFile invalid Pin
Member 342179027-Jun-14 0:58
memberMember 342179027-Jun-14 0:58 
AnswerRe: File invalid Pin
John D. Cook28-Jun-14 11:05
memberJohn D. Cook28-Jun-14 11:05 
GeneralRe: File invalid Pin
Member 342179029-Jun-14 22:54
memberMember 342179029-Jun-14 22:54 
GeneralThank you! Pin
JoraSteelRat1-Oct-13 1:04
memberJoraSteelRat1-Oct-13 1:04 
GeneralMy vote of 5 Pin
Koh Twee Toong29-Jun-13 8:39
memberKoh Twee Toong29-Jun-13 8:39 
QuestionGreat code (and a quick question) Pin
Member 992827620-Mar-13 11:43
memberMember 992827620-Mar-13 11:43 
GeneralMy vote of 5 Pin
silleryxu8-Oct-12 23:53
membersilleryxu8-Oct-12 23:53 
GeneralMy vote of 5 Pin
burndavid9327-May-12 6:20
memberburndavid9327-May-12 6:20 
GeneralMy vote of 5 Pin
MattEFowler12-Apr-12 11:27
memberMattEFowler12-Apr-12 11:27 
Questiona function to test Pin
dariusz_brzezinski23-Mar-12 0:23
memberdariusz_brzezinski23-Mar-12 0:23 
Maybe you could test your code on following integrand: (1-x)^0.99 in the range (0,1)? I tested it and noticed, that in this case very good accuracy offered by this method decrese dramatically. Even if we eliminate the singularity at the end point.
I would like to hear your comment, opinion, suggestions.
Best regards
AnswerRe: a function to test Pin
John D. Cook23-Mar-12 1:10
memberJohn D. Cook23-Mar-12 1:10 
GeneralReally good work [modified] Pin
V. Thieme8-Nov-10 11:27
memberV. Thieme8-Nov-10 11:27 
Generalvery good, thank you. Pin
jemisha28-Jul-10 7:27
memberjemisha28-Jul-10 7:27 
GeneralKrylov method Pin
Majid Shahabfar6-Jul-09 23:06
memberMajid Shahabfar6-Jul-09 23:06 
GeneralRe: Krylov method Pin
John D. Cook7-Jul-09 1:58
mvpJohn D. Cook7-Jul-09 1:58 
QuestionMultiple integration routine? Pin
Graeme Addison18-Feb-09 12:30
memberGraeme Addison18-Feb-09 12:30 
AnswerRe: Multiple integration routine? Pin
John D. Cook18-Feb-09 16:52
mvpJohn D. Cook18-Feb-09 16:52 
GeneralVery cool! Pin
Dr.Luiji16-Jan-09 23:52
memberDr.Luiji16-Jan-09 23:52 
GeneralRe: Very cool! Pin
John D. Cook17-Jan-09 4:00
mvpJohn D. Cook17-Jan-09 4:00 

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.

Permalink | Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.170215.1 | Last Updated 28 Jun 2014
Article Copyright 2008 by John D. Cook
Everything else Copyright © CodeProject, 1999-2017
Layout: fixed | fluid