emher
emher

Reputation: 6024

Multidimensional Integration - Coupled Limits

I need to calculate the value of a high dimensional integral in C++. I have found numerous libraries capable of solving this task for fixed limit integrals,

\int_{0}^{L} \int_{0}^{L} dx dy f(x,y) .

However the integrals which I am looking at have variable limits,

\int_{0}^{L} \int_{x}^{L} dx dy f(x,y) .

To clarify what i mean, here is a naive 2D Riemann sum implementation in 2D, which returns the desired result,

int steps = 100;
double integral = 0;
double dl = L/((double) steps);
double x[2] = {0};

for(int i = 0; i < steps; i ++){
    x[0] = dl*i;
    for(int j = i; j < steps; j ++){
        x[1] = dl*j;
        double val = f(x);
        integral += val*val*dl*dl;
    }
}

where f is some arbitrary function and L the common upper integration limit. While this implementation works, it's slow and thus impractical for higher dimensions.

Effective algorithms for higher dimensions exist, but to my knowledge, library implementations (e.g. Cuba) take a fixed value vector as the limit argument which renders them useless for my problem.

Is there any reason for this and/or is there any trick to circumvent the problem?

Upvotes: 2

Views: 1350

Answers (3)

user3717023
user3717023

Reputation:

When integrating over a triangular region such as 0<=x<=y<=L one can take advantage of symmetry: integrate f(min(x,y),max(x,y)) over the square 0<=x,y<=L and divide the result by 2. This has an advantage over extending f by zero (the method mentioned by LutzL) in that the extended function is continuous, which improves the performance of the integration routine.

I compared these on the example of the integral of 2x+y over 0<=x<=y<=1. The true value of the integral is 2/3. Let's compare the performance; for demonstration purpose I use Matlab routine, but this is not specific to language or library used.

Extending by zero

f = @(x,y) (2*x+y).*(x<=y);
result = integral2(f, 0, 1, 0, 1);
fprintf('%.9f\n',result);

Output:

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

0.666727294

Extending by symmetry

g = @(x,y) (2*min(x,y)+max(x,y));
result2 = integral2(g, 0, 1, 0, 1)/2;
fprintf('%.9f\n',result2);

Output:

0.666666776

The second result is 500 times more accurate than the first.

Unfortunately, this symmetry trick is not available for general domains; but integration over a triangle comes up often enough so it's useful to keep it in mind.

Upvotes: 1

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

Your integration order is wrong, should be dy dx.


You are integrating over the triangle

0 <= x <= y <= L

inside the square [0,L]x[0,L]. This can be simulated by integrating over the full square where the integrand f is defined as 0 outside of the triangle. In many cases, when f is defined on the full square, this can be accomplished by taking the product of f with the indicator function of the triangle as new integrand.

Upvotes: 1

Spektre
Spektre

Reputation: 51873

I was a bit confused by your integral definition but from your code i see it like this:

integral

just did some testing so here is your code:

//---------------------------------------------------------------------------
double f(double *x) { return (x[0]+x[1]); }
void integral0()
    {
    double L=10.0;
    int steps = 10000;
    double integral = 0;
    double dl = L/((double) steps);
    double x[2] = {0};
    for(int i = 0; i < steps; i ++){
        x[0] = dl*i;
        for(int j = i; j < steps; j ++){
            x[1] = dl*j;
            double val = f(x);
            integral += val*val*dl*dl;
        }
    }
    }
//---------------------------------------------------------------------------

Here is optimized code:

//---------------------------------------------------------------------------
void integral1()
    {
    double L=10.0;
    int i0,i1,steps = 10000;
    double x[2]={0.0,0.0};
    double integral,val,dl=L/((double)steps);
    #define f(x) (x[0]+x[1])
    integral=0.0;
    for(x[0]= 0.0,i0= 0;i0<steps;i0++,x[0]+=dl)
    for(x[1]=x[0],i1=i0;i1<steps;i1++,x[1]+=dl)
        {
        val=f(x);
        integral+=val*val;
        }
    integral*=dl*dl;
    #undef f
    }
//---------------------------------------------------------------------------

results:

[ 452.639 ms] integral0
[ 336.268 ms] integral1
  • so the increase in speed is ~ 1.3 times (on 32bit app on WOW64 AMD 3.2GHz)
  • for higher dimensions it will multiply
  • but still I think this approach is slow
  • The only thing to reduce complexity I can think of is algebraically simplify things
    • either by integration tables or by Laplace or Z transforms
    • but for that the f(*x) must be know ...
  • constant time reduction can of course be done
    • by the use of multi-threading
    • and or GPU ussage
    • this can give you N times speed increase
    • because this is all directly parallelisable

Upvotes: 0

Related Questions