Performing 2d numerical integration with Boost Cpp

I've been learning to use the numerical quadrature of the Boost library for Cpp. In the documentation, I've found an example for 1D Gauss-Kronrod Quadrature using Boost.

#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <iostream>

int main(int argc, char *argv[])
{

auto f1 = [](double t) { return std::exp(-t*t / 2); };

double error;
double Q = gauss_kronrod<double, 61>::integrate(f1, 0, std::numeric_limits<double>::infinity(), 5, 1e-14, &error);
std::cout << Q << " " << error<<std::endl;
return 0;
}


This calculates the integral

$$Q = \int\limits_{0}^{\infty} e^{-t^2} dt$$

up to some tolerance 1e-14 and with max depth 5. I'm trying to build upon this and calculate a double integral (cubature) like this

$$Q = \int\limits_{0}^{\infty}\int\limits_{0}^{\infty} e^{-(t^2+ s^2+ts)} dt\ ds$$

I've been trying to calculate this double integral using two 1d integrals. I tried to look up if there is any double integtal function, but found none, so this is what I've tried so far

#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <iostream>

int main(int argc, char *argv[])
{

auto f1 = [](double t, double s) { return std::exp(-(t*t+s*s+t*s) / 2); };

double error;
double Q = gauss_kronrod<double, 15>::integrate(gauss_kronrod<double, 61>::integrate(f1, 0, std::numeric_limits<double>::infinity(), 5), 0, std::numeric_limits<double>::infinity(), 5, 1e-9, &error);
std::cout << Q << " " << error<<std::endl;
return 0;
}


Here, I've introduced an integral of 2 variables. Two 1d quadrature functions are nested. But this gives compilation error. What am I doing wrong and is there any other cleaner way of doing 2d-numerical integrations?

The reason you didn't find 2D quadrature is that we haven't implemented it yet. As to your compilation failure, this does the trick:

#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <iostream>

int main(int argc, char *argv[])
{

auto f1 = [](double t, double s) { return std::exp(-(t*t+s*s+t*s) / 2); };
auto f = [&](double t) {
auto g = [&](double s) {
return f1(t, s);
};
return gauss_kronrod<double, 61>::integrate(g, 0, std::numeric_limits<double>::infinity(), 5);
};

double error;
double Q = gauss_kronrod<double, 15>::integrate(f, 0, std::numeric_limits<double>::infinity(), 5, 1e-9, &error);
std::cout << Q << ", error estimated at " << error <<std::endl;
return 0;
}


That spits out 1.2092 for me.

• This works fine and should work out for any convergent integral, Thanks for your time. Will Boost implement cubatures in near future, do you have any information regarding this? – Galilean May 11 at 4:23
• It's on the roadmap, but it's behind multivariate interpolation. Check out the Padua points which looks pretty easy to implement. – user14717 May 11 at 12:31