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?