Four dimensional integration in cpp using gsl

[CAF123]

I am trying to perform a four dimensional integration in cpp using gsl. I do this by nesting together the one-dimensional integration routines from gsl library. What I wrote seems to work for a few test integrands but I am having trouble with the integrand that I actually want a value for. See below for the full code I am using.

#include <iostream>
#include <gsl/gsl_integration.h>
#include <math.h>
#include <complex>

double b = 7.2;
double c = 0.130477;
double d = 0.549;

struct Params_t {
    double z1, z2, bx1, by1, bx, by;
};

double f (double z1, void * params) {
    
    Params_t& Params = *(Params_t *)(params);
    double z2 = Params.z2;
    double bx1 = Params.bx1;
    double by1 = Params.by1;
    double bx = Params.bx;
    double by = Params.by;
    
    double r1 = sqrt( z1*z1 + bx1*bx1 + by1*by1 );
    double r2 = sqrt( z2*z2 + bx*bx-bx1*bx1 + by*by-by1*by1 );
    double exp1 = c/(1. + exp( (r1-b)/d) );
    double exp2 = c/(1. + exp( (r2-b)/d) );
    
    return exp1*exp2;
    
}

double firstIntegrand(double z2, void * params) {   //does the integration over z1 and sets up function ready for integration over z2 (with bx1,by1 kept fixed)
    
    Params_t& Params = *(Params_t *)(params);
    
    double epsabs = 1e-5;
    double epsrel = 1e-5;
    double lower = 0.;
    double upper = 10;
    
    double result, error;
    
    Params_t custom_Params;
    custom_Params.z2=z2;
    
    double bx1 = Params.bx1;
    double by1 = Params.by1;
    double bx = Params.bx;
    double by = Params.by;
    
    custom_Params.bx1=bx1;
    custom_Params.by1=by1;
    custom_Params.bx=bx;
    custom_Params.by=by;
    
    
    gsl_integration_workspace * ws = gsl_integration_workspace_alloc (3e8);

      gsl_function F;
      F.function = &f;
      F.params = &custom_Params;

      gsl_integration_qags (&F, lower, upper, epsabs, epsrel, 3e8,
                            ws, &result, &error);

      gsl_integration_workspace_free (ws);
    
    return result;
}

double secondIntegrand(double bx1, void * params) { //does the integration over z2 and sets up function ready for integration over bx1 (with by1 kept fixed)
    
    Params_t& Params = *(Params_t *)(params);
    
    double epsabs = 0;
    double epsrel = 1e-7;
    double lower = 0;
    double upper = 10;
    
    double result, error;
    
    Params_t custom_Params;
    custom_Params.bx1=bx1;
    
    
    double by1 = Params.by1;
    double bx = Params.bx;
    double by = Params.by;
    
    custom_Params.by1=by1;
    custom_Params.bx=bx;
    custom_Params.by=by;
    
    
    gsl_integration_workspace * ws = gsl_integration_workspace_alloc (3e8);

      gsl_function F;
      F.function = &firstIntegrand;
      F.params = &custom_Params;

      gsl_integration_qags (&F, lower, upper, epsabs, epsrel, 3e8,
                            ws, &result, &error);

      gsl_integration_workspace_free (ws);
    
    return result;
}

double thirdIntegrand(double by1, void * params) { //does the integration over bx1 and sets up function ready for integration over by1
    
    Params_t& Params = *(Params_t *)(params);
    
    double epsabs = 0;
    double epsrel = 1e-7;
    double lower = 0.;
    double upper = 10;
    
    double result, error;
    
    Params_t custom_Params;
    custom_Params.by1 = by1;
    
    double bx = Params.bx;
    double by = Params.by;
    
    custom_Params.bx=bx;
    custom_Params.by=by;
    
    gsl_integration_workspace * ws = gsl_integration_workspace_alloc (3e8);

      gsl_function F;
      F.function = &secondIntegrand;
      F.params = &custom_Params;

      gsl_integration_qags (&F, lower, upper, epsabs, epsrel, 3e8,
                            ws, &result, &error);

      gsl_integration_workspace_free (ws);
    
    return result;
}

double fourthIntegrand(double bx, double by) { //does the integration over by1 finally
    
    double epsabs = 0.;
    double epsrel = 1e-7;
    double lower = 0.;
    double upper = 10;
    
    double result, error;
    
    Params_t custom_Params;
    
    custom_Params.bx=bx;
    custom_Params.by=by;
    
    
    gsl_integration_workspace * ws = gsl_integration_workspace_alloc (3e8);

    gsl_function F;
    F.function = &thirdIntegrand;
    F.params = &custom_Params;

      gsl_integration_qags (&F, lower, upper, epsabs, epsrel, 3e8,
                            ws, &result, &error);

      gsl_integration_workspace_free (ws);
    
    return result;
}

int main(int argc, const char * argv[]) {
     
     std::cout << fourthIntegrand(0.2,0.3) << std::endl; //prints out val of four dim integration for a particular value of bx and by
   
 }

The code works okay for simpler integrands than the exp1*exp2 declared in the opening function declaration. E.g. for bx*by*bx1*by1 or sqrt[z1*z2+bx1*by1]. The problem is in r2, and in particular due to the combination in the sqrt: bx*bx-bx1*bx1 and the fact there is a relative minus (change this to a + and code is OK). I have however been able to my integral in mathematica so I am trying to understand why there is a problem here in gsl. I have played with the relative and absolute error accuracy parameters. Thanks in advance for any comments about this and/or if there is a way to do a multi dimensional integration in gsl more efficiently.

 

[jedishrfu]

Is there any overflow or under flow issues that you can check for?

 

[pbuk]

Was the result not what you expected? Was it close but not accurate enough? Did it fail to (a) compile, (b) execute, or (c) terminate? Does it run too slowly? Did your computer burst into flames?

 

[CAF123]

Is there any overflow or under flow issues that you can check for?

Was the result not what you expected? Was it close but not accurate enough? Did it fail to (a) compile, (b) execute, or (c) terminate? Does it run too slowly? Did your computer burst into flames?

:P yes sorry I should have made that clear. I get the error

gsl: qags.c:553: ERROR: bad integrand behavior found in the integration interval
Default GSL error handler invoked.

originating in firstIntegrand function. So compilation is fine, but execution is not. Thanks.

 

[Mark44]

The code works okay for simpler integrands than the exp1*exp2 declared in the opening function declaration. E.g. for bx*by*bx1*by1 or sqrt[z1*z2+bx1*by1]. The problem is in r2, and in particular due to the combination in the sqrt: bx*bx-bx1*bx1 and the fact there is a relative minus (change this to a + and code is OK).

From post #1:

double r2 = sqrt( z2*z2 + bx*bx-bx1*bx1 + by*by-by1*by1 );

If there is a subinterval for which the expression in the call to sqrt() is negative, that seems to me possibly to be the source of the problem noted in the error message.

gsl: qags.c:553: ERROR: bad integrand behavior found in the integration interval

 

[CAF123]

From post #1:

double r2 = sqrt( z2*z2 + bx*bx-bx1*bx1 + by*by-by1*by1 );

If there is a subinterval for which the expression in the call to sqrt() is negative, that seems to me possibly to be the source of the problem noted in the error message.

Yes I agree, I had thought about this too but then I was thinking
a) why would the error show up in firstIntegrand when the integration there is over z1?
b) mathematica is able to do it and returns real number

 

[pasmith]

Yes I agree, I had thought about this too but then I was thinking
a) why would the error show up in firstIntegrand when the integration there is over z1?
b) mathematica is able to do it and returns real number

You are computing
$$\int_0^{10} \int_0^{10} \int_0^{10} \int_0^{10} f(x,y,x_1,y_1,z_1,z_2)\,dz_1\,dz_2\,dx_1\,dy_1.$$ The integrand $f$ is only evaluated when integrating over the innermost variable, which is done for all values of the outer variables. Thus if $f$ becomes complex for some value of an outer variable, the program will notice this when trying to integrate over the innermost variable at that value of the outer variable.

The GSL routines define their integrand as a gsl_function, which expects a double argument and returns a double value. They are not therefore designed to deal with complex intermediate results, even if the final result is real. Mathematica is designed to deal with complex values.

 

[CAF123]

OK I see, thanks. So the calculation must generate an even number of i factors in the evaluation of the integrals such that the end result is overall real, as given by mathematica. I actually prefer what gsl does because, without specifying what side of the branch cut you are on, likely mathematica can give you the wrong sign of the imaginary part in what it does.