Is it possible to optimise this integration code so that it runs faster?

Question

double trap(double func(double), double b, double a, double N) {
  double j;
  double s;
  double h = (b-a)/(N-1.0); //Width of trapezia

  double func1 = func(a);
  double func2;

  for (s=0,j=a;j<b;j+=h){
    func2 = func(j+h);
    s = s + 0.5*(func1+func2)*h;
    func1 = func2;
  }

  return s;
}

The above is my C++ code for a 1D numerical integration (using the extended trapezium rule) of func() between limits $[a,b]$ using $N-1$ trapezia.

I am actually doing a 3D integration, where this code is called recursively. I work with $N = 50$ giving me decent results.

Other than reducing $N$ further, is anybody able to suggest how to optimise the code above so that it runs faster? Or, even, can suggest a faster integration method?

Answer 1

Mathematically, your expression is equivalent to:

$$I = h \left(\frac{1}{2}f_1 + f_2 + f_3 +...+f_{n-1} + \frac{1}{2}f_n \right) + O\left(\frac{(b-a)^3 f''}{n^2} \right)$$

So you could implement that. As it was said, the time is probably dominated by the function evaluation, so to get the same accuracy, you can use a better integration method that requires less function evaluations.

Gaussian quadrature is, in modern days, bit more than a toy; only useful if you require very few evaluations. If you want something easy to implement, you can use Simpson's rule, but I wouldn't go further than order $1/N^3$ without a good reason.

If the function's curvature changes a lot, you could use an adaptative step routine, that would select a larger step when the function is flat, and a smaller more accurate one when the curvature is higher.

Answer 2

Chances are that the evaluation of the functions is the most time consuming part of this computation. If that's the case, then you should focus on improving the speed of func() rather than trying to speed up the integration routine itself.

Depending on the properties of func(), it's also likely that you could get a more precise evaluation of the integral with fewer function evaluations by using a more sophisticated integration formula.

Answer 3

Possible? Yes. Useful? No. The optimizations I'm going to list here are unlikely to make more than a tiny fraction of a percent difference in the runtime. A good compiler may already do these for you.

Anyway, looking at your inner loop:

    for (s=0,j=a;j<b;j+=h){
        func2 = func(j+h);
        s = s + 0.5*(func1+func2)*h;
        func1 = func2;
    }

At every loop iteration you perform three math operations that can be brought outside: adding j + h, multiplication by 0.5, and multiplication by h. The first you can fix by starting your iterator variable at a + h, and the others by factoring out the multiplications:

    for (s=0, j=a+h; j<=b; j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

Though I would point out that by doing this, due to floating point roundoff error it is possible to miss the last iteration of the loop. (This was also an issue in your original implementation.) To get around that, use an unsigned int or size_t counter:

    size_t n;
    for (s=0, n=0, j=a+h; n<N; n++, j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

As Brian's answer says, your time is better spent optimizing the evaluation of the function func. If the accuracy of this method is sufficient, I doubt you'll find anything faster for the same N. (Though you could run some tests to see if e.g. Runge-Kutta lets you lower N enough that the overall integration takes less time without sacrificing accuracy.)

Answer 4

There are several changes I would recommend to improve the computation:

• For performance and accuracy, use std::fma(), which performs a fused multiply-add.
• For performance, defer multiplying the area of each trapezoid by 0.5 — you can do it once at the end.
• Avoid repeated addition of h, which could accumulate round-off errors.

In addition, I would make several changes for clarity:

• Give the function a more descriptive name.
• Swap the order of a and b in the function signature.
• Rename N → n, h → dx, j → x2, s → accumulator.
• Change n to an int.
• Declare variables in a tighter scope.

#include <cmath>

double trapezoidal_integration(double func(double), double a, double b, int n) {
    double dx = (b - a) / (n - 1);   // Width of trapezoids

    double func_x1 = func(a);
    double accumulator = 0;

    for (int i = 1; i <= n; i++) {
        double x2 = a + i * dx;      // Avoid repeated floating-point addition
        double func_x2 = func(x2);
        accumulator = std::fma(func_x1 + func_x2, dx, accumulator); // Fused multiply-add
        func_x1 = func_x2;
    }

    return 0.5 * accumulator;
}

Answer 5

If your function is a polynomial, possibly weighted by some function (e.g. a gaussian), you can do an exact integration in 3d directly with a cubature formula (e.g. http://people.sc.fsu.edu/~jburkardt/c_src/stroud/stroud.html ) or with a sparse grid (e.g. http://tasmanian.ornl.gov/ ). These methods simply specify a set of points and weights to multiply the function value by, so they are very fast. If your function is smooth enough to be approximated by polynomials, then these methods can still give a very good answer. The formulas are specialized to the type of function you're integrating, so it may take some digging to find the right one.

Answer 6

When you try to calculate an integral numerically, you try to get the precision that you want with the smallest possible effort, or alternatively, try to get the highest possible precision with a fixed effort. You seem to ask how to make the code for one particular algorithm run as fast as possible.

That may give you some little gain, but it will be little. There are much more efficient methods for numerical integration. Google for "Simpson's rule", "Runge-Kutta", and "Fehlberg". They all work quite similar by evaluating some values of the function and cleverly adding multiples of those value, producing much smaller errors with the same number of function evaluations, or the same error with a much smaller number of evaluations.

Answer 7

There are lots of ways to do integration, of which the trapezoidal rule is about the simplest.

If you know anything at all about the actual function you're integrating, you can do better if you exploit that. The idea is to minimize the number of grid points within acceptable levels of error.

For example, trapezoidal is making a linear fit to consecutive points. You could make a quadratic fit, which if the curve is smooth would fit better, which could allow you to use a coarser grid.

Orbital simulations are sometimes done using conics, because orbits are very much like conic sections.

In my work, we are integrating shapes that approximate bell-shaped curves, so it is effective to model them as that (adaptive Gaussian quadrature is considered the "gold standard" in this work).

Answer 8

So, as has been pointed out in other answers, this depends heavily on how expensive your function is. Optimizing your trapz code is only worth it if it is really your bottleneck. If it's not completely obvious, you should check this by profiling your code (tools like Intels V-tune, Valgrind or Visual Studio can do this).

I would however suggest a completely different approach: Monte Carlo integration . Here you simply approximate the integral by sampling your function at random points adding the results. See this pdf in addition to the wiki page for details.

This is works extremely well for high dimensional data, typically much better than the quadrature methods used in 1-d integration.

The simple case is very easy to implement (see the pdf), just be careful that the standard random function in c++98 is quite bad both performance and quality wise. In c++11, you can use the Mersenne Twister in .

If your function has a lot of variation in some areas and less in others, consider using stratified sampling. I would recommend using the GNU scientific library, rather than writing your own though.

Answer 9

I am actually doing a 3D integration, where this code is called recursively.

"recursively" is the key. You are either going through a large data set and considering a lot of data more than once, or you are actually generating your data set yourself from (piecewise?) functions.

Recursively evaluated integrations will be ridiculously expensive, and ridiculously imprecise as the powers increase in recursion.

Create a model for interpolating your data set and do a piecewise symbolic integration. Since a lot of data is then collapsing into coefficients of base functions, the complexity for deeper recursion grows polynomially (and usually rather low powers) rather than exponentially. And you get "exact" results (you still need to figure out good evaluation schemes to get reasonable numeric performance, but it should still be rather feasible to get better than trapezoidal integration).

If you take a look at the error estimates for trapezoidal rules, you'll find that they are related to some derivative of the involved functions, and if the integration/definition is done recursively, the functions will not tend to have well-behaved derivatives.

If your only tool is a hammer, every problem looks like a nail. While you barely touch upon the problem in your description, I have the suspicion that applying the trapezoidal rule recursively is a bad match: you get an explosion of both inaccuracy and computational requirements.

Answer 10

the original code evaluates the function at each N points, then adds the values up, and multiplies the sum by the step size. the only trick is that the values at the beginning and the end are added with weight $1/2$, while all points inside are added with full weight. actually, they are also added with weight $1/2$ but twice. instead of adding them twice, add them only once with full weight. factor out the multiplication by the step size outside of the loop. that's all that can be done to speed this up, really.

    double trap(double func(double), double b, double a, double N){
double j, s;
double h = (b-a)/(N-1.0); //Width of trapezia

double s = 0;
j = a;
for(i=1; i<N-1; i++){
  j += h;
  s += func(j);
}
s += (func(a)+func(b))/2;

return s*h;
}