# SPARSE_GRID_OPEN Sparse Grids Based on Open Quadrature Rules

SPARSE_GRID_OPEN is a FORTRAN90 library which computes the location of points on a sparse grid based on an open quadrature rule.

One way of looking at the construction of sparse grids is to assume that we start out by constructing a (very dense) product grid. We will assume for now that the order, that is, the number of points used in each component grid, is the same for all dimensions. Moreover, we will assume that the order is a power of 2 minus one, so that we have a natural relationship between the order and the logarithm base 2 of the order plus 1:

```
order = 2( level + 1 ) - 1
```

Thus, if we allow level to grow, the order roughly doubles, as follows:
LevelOrder
0 1
1 3
2 7
3 15
4 31
5 63
6 127
7 255
8 511
9 1023
10 2047

To keep things simple, let us begin by supposing we are selecting points for a grid to be used in an interpolation or quadrature rule. If you successively compute the locations of the points of each level, you will probably see that the points of a level are all included in the grid associated with the next level. (This is not guaranteed for all rules; it's simply a property of the way most such grids are defined!).

This nesting property is very useful. For one thing, it means that when if we've computed a grid of one level, and now proceed to the next, then all the information associated with the current level (point location, the value of functions at those points) is still useful for the next level, and will save us some computation time as well. This also means that, when we have reached a particular level, all the previous levels are still available to us, with no extra storage. These considerations make it possible, for instance, to do efficient and convenient error estimation.

When we move to a problem whose geometry is two-dimensional or more, we can still take the same approach. However, when working in multidimensional geometry, it is usually not a good idea to form a grid using the product of 1D grids, especially when we are determining the order using the idea of levels. Especially in this case, if we go to the next level in each dimension, the total number of points would increase by a factor of roughly 2 to the spatial dimension. Just a few such steps in, say, 6 dimensions, and we would be far beyond our computational capacity.

Instead, in multidimensions, the idea is to construct a sparse grid, which can be thought of in one of two ways:

• the sparse gird is a logical sum of low order product grids; each product grid has a total level (sum of the levels of the 1d rules) that is less than or equal to level_max;
• the sparse grid is a very sparse selection of points from the very high order product grid formed by using rules of level level_max in each dimension.

(There is still a lot of explaining to do to get from the one-dimensional levels to the N-dimensional levels and the selection of the low-level product grids that sum up to the sparse grid...)

Once the grid indices of the sparse grid points have been selected, there are a variety of schemes for distributing the points. We concentrate on two sets of schemes:

Uniform (equally spaced) open rules of order N over [0,1]:

• Newton Cotes Open: [ 1, 2, 3,..., N ] / ( N + 1 );
does not include boundaries. Point spacing is 1/(N+1)
First and last points are distance 1/N to the boundary.

The uniform schemes are easy to understand. However, it has been observed that greater accuracy and stability can be achieved by arranging the points in a nonuniform way that tends to move points towards the boundary and away from the center. A common scheme uses the cosine function to do this, and can be naturally derived from the uniform schemes.

Nonuniform open rules of order N over [-1,1]:

• Fejer Type 2:
Theta = pi * [ 1, 2, 3,..., N ] / ( N + 1 );
Points = cos ( theta );
does not include boundaries.
Analogue to the Newton Cotes Open Rule.
• Gauss-Patterson:
Starting with midpoint rule, then the Gauss rule of order 3,
incremented at each new level by intermediate points,
does not include boundaries.

Note that a standard Gauss-Legendre quadrature rule will not be suitable for use in constructing sparse grids, because rules of successively greater levels are not naturally nested.

### Languages:

SPARSE_GRID_OPEN is available in a C++ version and a FORTRAN90 version and a MATLAB version.

### Related Data and Programs:

CC_DISPLAY, a MATLAB library which can compute and display Clenshaw Curtis grids in two dimensions, as well as sparse grids formed from sums of Clenshaw Curtis grids.

QUADRATURE_RULES, a dataset directory which defines quadrature rules; a number of examples of sparse grid quadrature rules are included.

QUADRULE, a FORTRAN90 library which defines quadrature rules for various intervals and weight functions.

SGMGA, a FORTRAN90 library which creates sparse grids based on a mixture of 1D quadrature rules, allowing anisotropic weights for each dimension.

SMOLPACK, a C library which implements Novak and Ritter's method for estimating the integral of a function over a multidimensional hypercube using sparse grids.

SPARSE_GRID_CC, a dataset directory which contains the abscissas of sparse grids based on a Clenshaw Curtis rule.

SPARSE_GRID_F2, a dataset directory which conains the abscissas of sparse grids based on a Fejer Type 2 rule.

SPARSE_GRID_GL, a FORTRAN90 library which creates sparse grids based on Gauss-Legendre rules.

SPARSE_GRID_GP, a dataset directory which contains the abscissas of sparse grids based on a Gauss Patterson rule.

SPARSE_GRID_HERMITE, a FORTRAN90 library which creates sparse grids based on Gauss-Hermite rules.

SPARSE_GRID_HW, a FORTRAN90 library which creates sparse grids based on Gauss-Legendre, Gauss-Hermite, Gauss-Patterson, or a nested variation of Gauss-Hermite rules, by Florian Heiss and Viktor Winschel.

SPARSE_GRID_LAGUERRE, a FORTRAN90 library which creates sparse grids based on Gauss-Laguerre rules.

SPARSE_GRID_MIXED, a FORTRAN90 library which constructs a sparse grid using different rules in each spatial dimension.

SPARSE_GRID_NCC, a dataset directory which contains files of the abscissas of sparse grids based on a Newton Cotes closed rule.

SPARSE_GRID_NCO, a dataset directory which contains the abscissas of sparse grids based on a Newton Cotes open rule.

SPINTERP, a MATLAB library which carries out piecewise multilinear hierarchical sparse grid interpolation; an earlier version of this software is ACM TOMS Algorithm 847, by Andreas Klimke;

TOMS847, a MATLAB library which carries out piecewise multilinear hierarchical sparse grid interpolation; this library is commonly called SPINTERP (version 2.1); this is ACM TOMS Algorithm 847, by Andreas Klimke;

### Reference:

1. Volker Barthelmann, Erich Novak, Klaus Ritter,
High Dimensional Polynomial Interpolation on Sparse Grids,
Volume 12, Number 4, 2000, pages 273-288.
2. Thomas Gerstner, Michael Griebel,
Numerical Integration Using Sparse Grids,
Numerical Algorithms,
Volume 18, Number 3-4, 1998, pages 209-232.
3. Albert Nijenhuis, Herbert Wilf,
Combinatorial Algorithms for Computers and Calculators,
Second Edition,
ISBN: 0-12-519260-6,
LC: QA164.N54.
4. Fabio Nobile, Raul Tempone, Clayton Webster,
A Sparse Grid Stochastic Collocation Method for Partial Differential Equations with Random Input Data,
SIAM Journal on Numerical Analysis,
Volume 46, Number 5, 2008, pages 2309-2345.
5. Sergey Smolyak,
Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions,
Volume 4, 1963, pages 240-243.
6. Dennis Stanton, Dennis White,
Constructive Combinatorics,
Springer, 1986,
ISBN: 0387963472,
LC: QA164.S79.

### List of Routines:

• ABSCISSA_LEVEL_OPEN_ND: first level at which given abscissa is generated.
• COMP_NEXT computes the compositions of the integer N into K parts.
• F2_ABSCISSA returns the I-th abscissa for the Fejer type 2 rule.
• GET_UNIT returns a free FORTRAN unit number.
• GL_ABSCISSA sets abscissas for "nested" Gauss-Legendre quadrature.
• GP_ABSCISSA returns the I-th abscissa for a Gauss Patterson rule.
• I4_CHOOSE computes the binomial coefficient C(N,K).
• I4_MODP returns the nonnegative remainder of I4 division.
• INDEX_TO_LEVEL_OPEN determines the level of a point given its index.
• LEVEL_TO_ORDER_OPEN converts a level to an order for open rules.
• LEVELS_OPEN_INDEX computes open grids with 0 <= LEVEL <= LEVEL_MAX.
• MULTIGRID_INDEX1 returns an indexed multidimensional grid.
• MULTIGRID_SCALE_OPEN renumbers a grid as a subgrid on a higher level.
• NCO_ABSCISSA returns the I-th abscissa for the Newton Cotes Open rule.
• R8MAT_WRITE writes an R8MAT file.
• S_BLANK_DELETE removes blanks from a string, left justifying the remainder.
• SPARSE_GRID_F2S_SIZE sizes a sparse grid using Fejer Type 2 Slow rules.
• SPARSE_GRID_GPS_SIZE sizes a sparse grid using Gauss-Patterson-Slow rules.
• SPARSE_GRID_OFN_SIZE sizes a sparse grid using Open Fully Nested rules.
• SPARSE_GRID_ONN_SIZE sizes a sparse grid using Open Non-Nested rules.
• SPARSE_GRID_OWN_SIZE sizes a sparse grid using Open Weakly Nested rules.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• VEC_COLEX_NEXT2 generates vectors in colex order.

You can go up one level to the FORTRAN90 source codes.

Last revised on 12 January 2010.