11,642,855 members (69,760 online)

# Fast Numerical Integration

, 28 Jun 2014 Public Domain 88.3K 1.8K 93
 Rate this:
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:

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

```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:

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
{
public:
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.

We implement a function object for our integrand as follows:

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

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

## Share

 Singular Value Consulting 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...

 First Prev Next
 Where is the source code? Member 1168303211-May-15 19:30 Member 11683032 11-May-15 19:30
 Re: Where is the source code? John D. Cook12-May-15 1:54 John D. Cook 12-May-15 1:54
 Re: Where is the source code? Member 116830326-Jun-15 14:34 Member 11683032 6-Jun-15 14:34
 My vote of 5! Volynsky Alex12-Nov-14 4:05 Volynsky Alex 12-Nov-14 4:05
 My vote of 5 CatchExAs9-Jul-14 3:55 CatchExAs 9-Jul-14 3:55
 File invalid Member 342179026-Jun-14 23:58 Member 3421790 26-Jun-14 23:58
 Re: File invalid John D. Cook28-Jun-14 10:05 John D. Cook 28-Jun-14 10:05
 Re: File invalid Member 342179029-Jun-14 21:54 Member 3421790 29-Jun-14 21:54
 Thank you! JoraSteelRat1-Oct-13 0:04 JoraSteelRat 1-Oct-13 0:04
 My vote of 5 Koh Twee Toong29-Jun-13 7:39 Koh Twee Toong 29-Jun-13 7:39
 Great code (and a quick question) Member 992827620-Mar-13 10:43 Member 9928276 20-Mar-13 10:43
 My vote of 5 silleryxu8-Oct-12 22:53 silleryxu 8-Oct-12 22:53
 My vote of 5 burndavid9327-May-12 5:20 burndavid93 27-May-12 5:20
 My vote of 5 MattEFowler12-Apr-12 10:27 MattEFowler 12-Apr-12 10:27
 a function to test dariusz_brzezinski22-Mar-12 23:23 dariusz_brzezinski 22-Mar-12 23:23
 Re: a function to test John D. Cook23-Mar-12 0:10 John D. Cook 23-Mar-12 0:10
 Really good work [modified] V. Thieme8-Nov-10 10:27 V. Thieme 8-Nov-10 10:27
 very good, thank you. jemisha28-Jul-10 6:27 jemisha 28-Jul-10 6:27
 Krylov method Majid Shahabfar6-Jul-09 22:06 Majid Shahabfar 6-Jul-09 22:06
 Re: Krylov method John D. Cook7-Jul-09 0:58 John D. Cook 7-Jul-09 0:58
 Re: Multiple integration routine? John D. Cook18-Feb-09 15:52 John D. Cook 18-Feb-09 15:52
 Very cool! Dr.Luiji16-Jan-09 22:52 Dr.Luiji 16-Jan-09 22:52
 Re: Very cool! John D. Cook17-Jan-09 3:00 John D. Cook 17-Jan-09 3:00
 Why function objects John D. Cook19-Dec-08 6:12 John D. Cook 19-Dec-08 6:12
 code also available for g++ andreas19808-Dec-08 22:26 andreas1980 8-Dec-08 22:26
 Re: code also available for g++ John D. Cook8-Dec-08 23:51 John D. Cook 8-Dec-08 23:51
 Re: code also available for g++ andreas198010-Dec-08 5:56 andreas1980 10-Dec-08 5:56
 Accuracy Member 27859438-Dec-08 14:10 Member 2785943 8-Dec-08 14:10
 Re: Accuracy John D. Cook8-Dec-08 15:38 John D. Cook 8-Dec-08 15:38
 Mathematical details posted [modified] John D. Cook8-Dec-08 3:45 John D. Cook 8-Dec-08 3:45
 Symbolic counterpart? Michael Mogensen8-Dec-08 2:03 Michael Mogensen 8-Dec-08 2:03
 Re: Symbolic counterpart? John D. Cook8-Dec-08 3:39 John D. Cook 8-Dec-08 3:39
 Thanks, but need more details Hatem Mostafa7-Dec-08 21:49 Hatem Mostafa 7-Dec-08 21:49
 Re: Thanks, but need more details John D. Cook8-Dec-08 0:43 John D. Cook 8-Dec-08 0:43
 3D? Johann Gerell7-Dec-08 20:11 Johann Gerell 7-Dec-08 20:11
 Re: 3D? John D. Cook8-Dec-08 0:48 John D. Cook 8-Dec-08 0:48
 Re: 3D? Johann Gerell8-Dec-08 20:39 Johann Gerell 8-Dec-08 20:39
 Brilliant f27-Dec-08 5:22 f2 7-Dec-08 5:22
 Re: Brilliant [modified] John D. Cook7-Dec-08 8:26 John D. Cook 7-Dec-08 8:26
 Re: Brilliant f27-Dec-08 9:39 f2 7-Dec-08 9:39
 great Luc Pattyn6-Dec-08 16:11 Luc Pattyn 6-Dec-08 16:11
 Last Visit: 31-Dec-99 18:00     Last Update: 2-Aug-15 18:09 Refresh 1