SG++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
quadrature.cpp

The following example shows how to integrate in SG++, using both direct integration of a sparse grid function and the use of Monte Carlo integration.

As in the tutorial.cpp (Start Here) example, we deal with the function

\[ f\colon [0, 1]^2 \to \mathbb{R},\quad f(x_0, x_1) := 16 (x_0 - 1) x_0 (x_1 - 1) x_1 \]

which we first interpolate. We then integrate the interpolant, then the function itself using 100000 Monte Carlo points, and we then compute the L2-error.

For instructions on how to compile and run the example, please see Installation and Usage.

The function, which sgpp::base::OperationQuadratureMC takes, has three parameters. First, the dimensionality (int), then a double* with the coordinates of the grid point \(\in[0,1]^d\), and finally a void* with clientdata for the function, see sgpp::base::FUNC.

// include all SG++ base headers
#include <sgpp_base.hpp>
#include <iostream>
// function to interpolate
double f(int dim, double* x, void* clientdata) {
double res = 1.0;
for (int i = 0; i < dim; i++) {
res *= 4.0 * x[i] * (1.0 - x[i]);
}
return res;
}
int main() {
// create a two-dimensional piecewise bi-linear grid
int dim = 2;
std::unique_ptr<sgpp::base::Grid> grid(sgpp::base::Grid::createLinearGrid(dim));
sgpp::base::GridStorage& gridStorage = grid->getStorage();
std::cout << "dimensionality: " << gridStorage.getDimension() << std::endl;
// create regular grid, level 3
int level = 3;
grid->getGenerator().regular(level);
std::cout << "number of grid points: " << gridStorage.getSize() << std::endl;
// create coefficient vector
sgpp::base::DataVector alpha(gridStorage.getSize());
double p[2];
for (size_t i = 0; i < gridStorage.getSize(); i++) {
sgpp::base::GridPoint& gp = gridStorage.getPoint(i);
p[0] = gp.getStandardCoordinate(0);
p[1] = gp.getStandardCoordinate(1);
alpha[i] = f(2, p, NULL);
}
alpha);
// direct quadrature
std::unique_ptr<sgpp::base::OperationQuadrature> opQ(
double res = opQ->doQuadrature(alpha);
std::cout << "exact integral value: " << res << std::endl;
// Monte Carlo quadrature using 100000 paths
sgpp::base::OperationQuadratureMC opMC(*grid, 100000);
res = opMC.doQuadrature(alpha);
std::cout << "Monte Carlo value: " << res << std::endl;
res = opMC.doQuadrature(alpha);
std::cout << "Monte Carlo value: " << res << std::endl;
// Monte Carlo quadrature of a function
res = opMC.doQuadratureFunc(f, NULL);
std::cout << "MC value: " << res << std::endl;
// Monte Carlo quadrature of error
res = opMC.doQuadratureL2Error(f, NULL, alpha);
std::cout << "MC L2-error: " << res << std::endl;
} // end of main

This results in an output similar to:

dimensionality:        2
number of grid points: 17
exact integral value:  0.421875
Monte Carlo value:     0.421298
Monte Carlo value:     0.421971
MC value:              0.444917
MC L2-error:           0.0242639