Integrate a function of three variables C++

Question

I spent quiet some time looking on the internet to find a solution to this, maybe it's out there but nothing of what I saw helped me.

I have a function !

double integrand(double r, double phi, double theta) 

That I want to integrate with some given definite bounds over the three dimensions. I found multiple lines of code on the internet that implement single variable definite integrals numerical schemes. I was thinking to myself "well, I'll just integrate along one dimension after the other".

Algorithmically speaking what I wanted to do was :

double firstIntegral(double r, double phi) {
    double result = integrationFunction(integrand,lower_bound,upper_bound);
    return result;
}

And simply do it again two more times. This works easily in languages like Matlab where I can create functions handler anywhere but I don't know how to do it in C++. I would have to first define a function that some r and phi will calculate integrand(r, phi, theta) for any theta and make it in C++ a function of one variable only but I don't know how to do that.

How can I compute the triple integral of my three-variables function in C++ using a one -dimensional integration routine (or anything else really...) ?

Answer 1

If your concern is just about getting the right prototype to pass to the integration function, you can very well use alternative data passing mechanisms, the simpler of which is using global variables.

Assuming that the order of integration is on theta, then phi, then r, write three functions of a single argument:

It(theta) computes the integrand from the argument theta passed explicitly and the global phi and r.

Ip(phi) computes the bounds on theta from the argument phi passed explicitly and the global r; it also copies the phi argument to the global variable and invokes integrationFunction(It, lower_t, upper_t).

Ir(r) computes the bounds on phi from the argument r passed explicitly; it also copies the r argument to the global variable and invokes integrationFunction(Ip, lower_p, upper_p).

Now you are ready to call integrationFunction(Ir, lower_r, upper_r).

It may also be that integrationFunction supports a "context" argument where you can store what you want.

Answer 2

This is a very slow and inexact version for integrals over cartesian coordinates, which should work with C++11.

It is using std::function and lambdas to implement the numerical integration. No steps have been taken to optimize this.

A template based solution could be much faster (by several orders of magnitude) than this, because it may allow the compiler to inline and simplify some of the code.

#include<functional>
#include<iostream>

static double integrand(double /*x*/, double y, double /*z*/)
{
  return y;
}

double integrate_1d(std::function<double(double)> const &func, double lower, double upper)
{
  static const double increment = 0.001;

  double integral = 0.0;
  for(double x = lower; x < upper; x+=increment) {
    integral += func(x) * increment;
  }
  return integral;
}

double integrate_2d(std::function<double(double, double)> const &func, double lower1, double upper1, double lower2, double upper2)
{
  static const double increment = 0.001;

  double integral = 0.0;
  for(double x = lower2; x < upper2; x+=increment) {
    auto func_x = [=](double y){ return func(x, y);};
    integral += integrate_1d(func_x, lower1, upper1) * increment;
  }
  return integral;
}

double integrate_3d(std::function<double(double, double, double)> const &func,
                    double lower1, double upper1,
                    double lower2, double upper2,
                    double lower3, double upper3)
{
  static const double increment = 0.001;

  double integral = 0.0;
  for(double x = lower3; x < upper3; x+=increment) {
    auto func_x = [=](double y, double z){ return func(x, y, z);};
    integral += integrate_2d(func_x, lower1, upper1, lower2, upper2) * increment;
  }
  return integral;
}


int main()
{
  double integral = integrate_3d(integrand, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
  std::cout << "Triple integral: " << integral << std::endl;
  return 0;
}

Answer 3

You can use functors

#include <iostream>

struct MyFunctorMultiply
{
    double m_coeff;

    MyFunctorMultiply(double coeff)
    {
        m_coeff  = coeff;
    }

    double operator()(double value)
    {
        return m_coeff * value;
    }
};

struct MyFunctorAdd
{
    double m_a;

    MyFunctorAdd(double a)
    {
        m_a  = a;
    }

    double operator()(double value)
    {
        return m_a + value;
    }
};

template<class t_functor>
double calculate(t_functor functor, double value, double other_param)
{

    return functor(value) - other_param;
}

int main()
{
    MyFunctorMultiply multiply2(2.);
    MyFunctorAdd      add3(3.);

    double result_a  = calculate(multiply2, 4, 1); // should obtain 4 * 2 - 1 = 7
    double result_b  = calculate(add3, 5, 6);      // should obtain 5 + 3 - 6 = 2

    std::cout << result_a << std::endl;
    std::cout << result_b << std::endl;
}