gsl multivariable numerical integration

Question

I have been using the GSL integration package (gsl_integration.h) in attempt to integrate some multivariable function f(x,y,z) with respect to x between some limit [a,b]. In this example I have as a toy model: f(x,y,z) = xyz

I would like to output the result of this integral using a function

double integral(double a, double b, double y, double z){}

where I can keep y and z arbitrary up until this point. My attempts so far have mostly involved setting y, z equal to a constant in some predefined function and then using the

gsl_integration_qags

function to integrate over that function. However, I want to keep the values of y and z arbitrary until I define them in the above function. The closest I have got so far is as follows

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <gsl/gsl_integration.h>
#include<stdio.h>
#include<math.h>

#define PI 3.1415926535897

double integrand(double x, void * params){
  double y = *(double *) params;
  double z = *(double *) params;            // PROBLEM: this just sets z = y
  double tmp = x*z*y;
  return tmp;
}

double integral(double a, double b, double y, double z){
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);
  gsl_function F;
  F.function = &integrand;               // Set integrand
  F.params = &y, &z;             // Set the parameters you wish to treat as constant in the integration
  double result, error;
  gsl_integration_qags (&F, a, b, 0, 1e-7, 1000,
                        w, &result, &error);
  gsl_integration_workspace_free (w); // Free the memory
  return result;
}




int main(){
std::cout << "Result "<< integral(0.0,1.0,3.0,5.0)<<std::endl;

}

This gives an output of

Result 4.5

The code sets a = 0.0, b = 1.0, y = z = 3.0 -- I want a = 0.0, b = 1.0, y = 3.0, z = 5.0, which would give a result of 7.5.

I would like to stick to GSL integration rather than boost if possible. I have also consulted https://www.gnu.org/software/gsl/doc/html/integration.html but am none-the-wiser. Any advice would be greatly appreciated, thanks.

Answer 1

I'm taking a guess, but it seems to me you want

double params[] = { y, z };
F.params = params;

and

double y = ((double *)params)[0];
double z = ((double *)params)[1];