SG++-Doxygen-Documentation
Quadrature in Python

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.py (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.

The function, which sgpp::base::OperationQuadratureMC takes, has one parameter, a sequence (C++ provides a tuple) with the coordinates of the grid point $$\in [0,1]^d$$.

1 # import pysgpp library
2 import pysgpp
3
4 # the standard parabola (arbitrary-dimensional)
5 def f(x):
6  res = 1.0
7  for i in range(len(x)):
8  res *= 4.0*x[i]*(1.0-x[i])
9  return res
10
11 # a pyramid-like shape (arbitrary-dimensional)
12 def g(x):
13  res = 1.0
14  for i in range(len(x)):
15  res *= 2.0*min(x[i], 1.0-x[i])
16  return res
17

Create a two-dimensional piecewise bi-linear grid of level 3

1 dim = 2
2 grid = pysgpp.Grid.createLinearGrid(dim)
3 gridStorage = grid.getStorage()
4 print("dimensionality: {}".format(dim))
5
6 # create regular grid, level 3
7 level = 3
8 gridGen = grid.getGenerator()
9 gridGen.regular(level)
10 print("number of grid points: {}".format(gridStorage.getSize()))
11

Calculate the surplus vector alpha for the interpolant of $$f(x)$$. Since the function can be evaluated at any point. Hence. we simply evaluate it at the coordinates of the grid points to obtain the nodal values. Then we use hierarchization to obtain the surplus value.

1 # create coefficient vector
2 alpha = pysgpp.DataVector(gridStorage.getSize())
3 for i in range(gridStorage.getSize()):
4  gp = gridStorage.getPoint(i)
5  p = tuple([gp.getStandardCoordinate(j) for j in range(dim)])
6  alpha[i] = f(p)
7
8 pysgpp.createOperationHierarchisation(grid).doHierarchisation(alpha)
9

Now we compute and compare the quadrature using four different methods available in SG++.

1 # direct quadrature
2 opQ = pysgpp.createOperationQuadrature(grid)
3 res = opQ.doQuadrature(alpha)
4 print("exact integral value: {}".format(res))
5
6 # Monte Carlo quadrature using 100000 paths
7 opMC = pysgpp.OperationQuadratureMC(grid, 100000)
8 res = opMC.doQuadrature(alpha)
9 print("Monte Carlo value: {:.6f}".format(res))
10 res = opMC.doQuadrature(alpha)
11 print("Monte Carlo value: {:.6f}".format(res))
12
13
14 # Monte Carlo quadrature of a standard parabola
15 res = opMC.doQuadratureFunc(f)
16 print("MC value (f): {:.6f}".format(res))
17
18 # Monte Carlo quadrature of error
19 res = opMC.doQuadratureL2Error(f, alpha)
20 print("MC L2-error (f-u) {:.7f}".format(res))
21
22
23 # Monte Carlo quadrature of a piramidal function
24 res = opMC.doQuadratureFunc(g)
25 print( "MC value (g): {:.6f}".format(res))
26
27 # Monte Carlo quadrature of error
28 res = opMC.doQuadratureL2Error(g, alpha)
29 print( "MC L2-error (g-u) {:.7f}".format(res))

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