Performing 2d numerical integration with Boost Cpp

Question

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[])
{

using namespace boost::math::quadrature;

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[])
{

using namespace boost::math::quadrature;

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?

Answer 1

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[])
{
    using namespace boost::math::quadrature;

    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.