SG++

To be able to quickly start with a toolkit, it is often advantageous (not only for the impatient users), to look at some code examples first.
In this tutorial, we give a short example program which interpolates a bivariate function on a regular sparse grid. Identical versions of the example are given in all languages currently supported by SG++: C++, Python, Java, and MATLAB.
For instructions on how to compile and run the example, please see Installation and Usage.
In the example, we create a twodimensional regular sparse grid with level 3 (with grid points \(\vec{x}_j \in [0, 1]^2\)) using piecewise bilinear basis functions \(\varphi_j\colon [0, 1]^2 \to \mathbb{R}\). We then interpolate 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 \]
with
\[ u\colon [0, 1]^2 \to \mathbb{R},\quad u(x_0, x_1) := \sum_{j=0}^{N1} \alpha_j \varphi_j(x_0, x_1) \]
by calculating the coefficients \(\alpha_j\) such that \(u(\vec{x}_j) = f(\vec{x}_j)\) for all \(j\). This process is called hierarchization in sparse grid contexts; the \(\alpha_j\) are called (hierarchical) surpluses. Note that \(f\) vanishes at the boundary of the domain \([0, 1]^2\); therefore, we don't have to spend sparse grid points on the boundary. Finally, we evaluate the sparse grid function \(u\) at a point \(\vec{p} = (0.52, 0.73)\).
Before starting, the function \(f\), which we want to interpolate, is defined.
First, we create a twodimensional grid (type sgpp::base::Grid) with piecewise bilinear basis functions with the help of the factory method sgpp::base::Grid.createLinearGrid().
Then we obtain a reference to the grid's sgpp::base::GridStorage object which allows us, e.g., to access grid points, to obtain the dimensionality (which we print) and the number of grid points.
Now, we use a sgpp::base::GridGenerator to create a regular sparse grid of level 3. Thus, gridStorage.getSize()
returns 17, the number of grid points of a twodimensional regular sparse grid of level 3.
We create an object of type sgpp::base::DataVector which is essentially a wrapper around a double
array. The DataVector
is initialized with as many entries as there are grid points. It serves as a coefficient vector for the sparse grid interpolant we want to construct. As the entries of a freshly created DataVector
are not initialized, we set them to 0.0. (This is superfluous here as we initialize them in the next few lines anyway.)
The for
loop iterates over all grid points: For each grid point gp
, the corresponding coefficient \(\alpha_j\) is set to the function value at the grid point's coordinates which are obtained by getStandardCoordinate(dim)
. The current coefficient vector is then printed.
An object of sgpp::base::OperationHierarchisation is created and used to hierarchize the coefficient vector, which we print.
Finally, a second DataVector is created which is used as a point to evaluate the sparse grid function at. An object is obtained which provides an evaluation operation (of type sgpp::base::OperationEvaluation), and the sparse grid interpolant is evaluated at \(\vec{p}\), which is close to (but not exactly at) a grid point.
The example results in the following output:
dimensionality: 2 number of grid points: 17 length of alpha vector: 17 alpha before hierarchization: [1, 0.75, 0.75, 0.4375, 0.9375, 0.9375, 0.4375, 0.75, 0.75, 0.4375, 0.9375, 0.9375, 0.4375, 0.5625, 0.5625, 0.5625, 0.5625] alpha after hierarchization: [1, 0.25, 0.25, 0.0625, 0.0625, 0.0625, 0.0625, 0.25, 0.25, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625] u(0.52, 0.73) = 0.7696
It can be clearly seen that the surpluses decay with a factor of 1/4: On the first level, we obtain 1, on the second 1/4, and on the third 1/16 as surpluses.