Cubature formula

A formula for the approximate calculation of multiple integrals of the form

$$ I ( f ) = \ \int\limits _ \Omega p ( x) f ( x) dx. $$

The integration is performed over a set $ \Omega $ in the Euclidean space $ \mathbf R ^ {n} $, $ x = ( x _ {1} \dots x _ {n} ) $. A cubature formula is an approximate equality

$$ \tag{1 } I ( f ) \cong \ \sum _ {j = 1 } ^ { N } C _ {j} f ( x ^ {(} j) ). $$

The integrand is written as the product of two functions: the first, $ p ( x) $, is assumed to be fixed for each specific cubature formula and is known as a weight function; the second, $ f ( x) $, is assumed to belong to some fairly broad class of functions, e.g. continuous functions such that the integral $ I ( f ) $ exists. The sum on the right-hand side of (1) is called a cubature sum; the points $ x ^ {(} j) $ are known as the interpolation points (knots, nodes) of the formula, and the numbers $ C _ {j} $ as its coefficients. Usually $ x ^ {(} j) \in \Omega $, though this condition is not necessary. In order to compute the integral $ I ( f ) $ via formula (1), one need only calculate the cubature sum. If $ n = 1 $ formula (1) and the sum on its right-hand side are known as a quadrature formula and sum (see Quadrature formula).

Let $ \alpha = ( \alpha _ {1} \dots \alpha _ {n} ) $ be a multi-index, where the $ \alpha _ {i} $ are non-negative integers; let $ | \alpha | = \alpha _ {1} + \dots + \alpha _ {n} $; and let $ x ^ \alpha = x _ {1} ^ {\alpha _ {1} } \dots x _ {n} ^ {\alpha _ {n} } $ be a monomial of degree $ | \alpha | $ in $ n $ variables; let

$$ \mu = \ M ( n, m) = \ \frac{( n + m)! }{n!m! } $$

be the number of monomials of degree at most $ m $ in $ n $ variables; let $ \phi _ {j} ( x) $, $ j = 1, 2 \dots $ be an ordering of all monomials such that monomials of lower degree have lower subscript while the monomials of equal degree have been ordered arbitrarily, e.g. in lexicographical order. In this enumeration $ \phi _ {1} ( x) = 1 $, and the $ \phi _ {j} ( x) $, $ j = 1 \dots \mu $, include all monomials of degree at most $ m $. Let $ \phi ( x) $ be a polynomial of degree $ m $. The set of points in the complex space $ \mathbf C ^ {n} $ satisfying the equation $ \phi ( x) = 0 $ is known as an algebraic hypersurface of degree $ m $.

One way to construct cubature formulas is based on algebraic interpolation. The points $ x ^ {(} j) \in \Omega $, $ j = 1 \dots \mu $, are so chosen that they do not lie on any algebraic hypersurface of degree $ m $ or, equivalently, they are chosen such that the Vandermonde matrix

$$ V = \ [ \phi _ {1} ( x ^ {(} j) ) \dots \phi _ \mu ( x ^ {(} j) )] _ {j = 1 } ^ \mu $$

is non-singular. The Lagrange interpolation polynomial for a function $ f ( x) $ with knots $ x ^ {(} j) $ has the form

$$ {\mathcal P} ( x) = \ \sum _ {j = 1 } ^ \mu {\mathcal L} _ {j} ( x) f ( x ^ {(} j) ), $$

where $ {\mathcal L} _ {j} ( x) $ is the polynomial of the influence of the $ j $- th knot: $ {\mathcal L} _ {j} ( x ^ {(} i) ) = \delta _ {ij} $( $ \delta _ {ij} $ is the Kronecker symbol). Multiplying the approximate equality $ f ( x) \cong {\mathcal P} ( x) $ by $ p ( x) $ and integrating over $ \Omega $ leads to a cubature formula of type (1) with $ N = \mu $ and

$$ \tag{2 } C _ {j} = I ( {\mathcal L} _ {j} ),\ \ j = 1 \dots \mu . $$

The existence of the integrals (2) is equivalent to the existence of the moments of the weight function, $ p _ {i} = I ( \phi _ {i} ) $, $ i = 1 \dots \mu $. Here and below it is assumed that the required moments of $ p ( x) $ exist. A cubature formula (1) which has $ N = \mu $ knots not contained in any algebraic hypersurface of degree $ m $ and with coefficients defined by (2), is called an interpolatory cubature formula. Formula (1) has the $ m $- property if it is an exact equality whenever $ f ( x) $ is a polynomial of degree at most $ m $; an interpolatory cubature formula has the $ m $- property. A cubature formula (1) with $ N \leq \mu $ knots which has the $ m $- property is an interpolatory formula if and only if the matrix

$$ [ \phi _ {1} ( x ^ {(} j) ) \dots \phi _ \mu ( x ^ {(} j) ) ] _ {j = 1 } ^ {N} $$

has rank $ N $. This condition holds when $ n = 1 $, so that a quadrature formula with $ N \leq m + 1 $ knots that has the $ m $- property is an interpolatory formula. The actual construction of an interpolatory cubature formula reduces to a selection of the knots and a calculation of the coefficients. The coefficients $ C _ {j} $ are determined by the linear algebraic system of equations

$$ \sum _ {j = 1 } ^ \mu C _ {j} \phi _ {i} ( x ^ {(} j) ) = p _ {i} ,\ \ i = 1 \dots \mu , $$

which is simply the mathematical expression of the statement that (1) (with $ N = \mu $) is exact for all monomials of degree at most $ m $. The matrix of this system is precisely $ V ^ { \prime } $( $ = V ^ {T} $).

Now suppose it is necessary to construct a cubature formula (1) with the $ m $- property, but with less than $ \mu $ knots. Since this cannot be done by merely selecting the coefficients, not only the coefficients but also the knots are unknowns in (1), giving $ N ( n + 1) $ unknowns in total. Since the cubature formula must have the $ m $- property, one obtains $ \mu $ equations

$$ \tag{3 } \sum _ {j = 1 } ^ { N } C _ {j} \phi _ {i} ( x ^ {(} j) ) = p _ {i} ,\ \ i = 1 \dots \mu . $$

It is natural to require the number of unknowns to coincide with the number of equations: $ N ( n + 1) = \mu $. This equation gives a tentative estimate of the number of knots. If $ N = \mu /( n + 1) $ is not an integer, one puts $ N = [ \mu /( n + 1)] + 1 $, where $ [ \mu / ( n + 1)] $ denotes the integer part of $ \mu /( n + 1) $. A cubature formula with this number of knots need not always exist. If it does exist, its number of knots is $ 1/( n + 1) $ times the number of knots of an interpolatory cubature formula. In that case, however, the knots themselves and the coefficients are determined by the non-linear system of equations (3). In the method of undetermined parameters, one constructs a cubature formula by trying to give it a form that will simplify the system (3). This can be done when $ \Omega $ and $ p ( x) $ have symmetries. The positions of the knots are taken compatible with the symmetry of $ \Omega $ and $ p ( x) $, and in that case symmetric knots are assigned the same coefficients. The simplification of the system (3) involves a certain risk: While the original system (3) may be solvable, the simplified system need not be.

Example. Let $ \Omega = K _ {2} = \{ - 1 \leq x _ {1} , x _ {2} \leq 1 \} $, $ p ( x _ {1} , x _ {2} ) = 1 $. One is asked to construct a cubature formula with the $ 7 $- property; $ n = 2 $, $ \mu = M ( 2, 7) = 36 $, and 12 knots. The knots are located as follows. The first group of knots consists of the intersection points of the circle of radius $ a $, centred at the origin, with the coordinate axes. The second group consists of the intersection points of the circle of radius $ b $, also centred at the origin, with the straight lines $ x _ {1} = \pm x _ {2} $. The third group is constructed similarly, with radius denoted by $ c $. The coefficients assigned to knots of the same group are identical and are denoted by $ A, B, C $ for knots of the first, second and third group, respectively. This choice of knots and coefficients implies that the cubature formula will be exact for monomials $ x _ {1} ^ {i} x _ {2} ^ {j} $ in which at least one of $ i $ or $ j $ is odd. For the cubature formula to possess the $ 7 $- property, it will suffice to ensure that it is exact for $ 1 $, $ x _ {1} ^ {2} $, $ x _ {1} ^ {4} $, $ x _ {1} ^ {2} x _ {2} ^ {2} $, $ x _ {1} ^ {6} $, $ x _ {1} ^ {4} x _ {2} ^ {2} $. This yields a non-linear system of six equations in the six unknowns $ a, b, c $, $ A, B, C $. Solving this system, one obtains a cubature formula with positive coefficients and with knots lying in $ K _ {2} $.

Let $ G $ be a finite subgroup of the group of orthogonal transformations $ \textrm{ O } ( n) $ of the space $ \mathbf R ^ {n} $ which leave the origin fixed. A set $ \Omega $ and a function $ p ( x) $ are said to be invariant under $ G $ if $ g ( \Omega ) = \Omega $ and $ p ( g ( x)) = p ( x) $ for $ x \in \Omega $ and any $ g \in G $. The set of points of the form $ ga $, where $ a $ is a fixed point of $ \mathbf R ^ {n} $ and $ g $ runs through all elements of $ G $, is known as the orbit containing $ a $. A cubature formula (1) is said to be invariant under $ G $ if $ \Omega $ and $ p ( x) $ are invariant under $ G $ and if the set of knots is a union of orbits, with knots belonging to the same orbit being assigned identical coefficients. Examples of sets invariant under $ G $ are the entire space $ \mathbf R ^ {n} $, and any ball or sphere centred at the origin; if $ G $ is the group of transformations of a regular polyhedron $ U $ onto itself, then $ U $ is also invariant. Thus, one can speak of invariant cubature formulas when $ \Omega $ is $ \mathbf R ^ {n} $, a ball, a sphere, a cube or any regular polyhedron, and when $ p ( x) $ is any function invariant under $ G $, e.g. $ p ( r) $, where $ r = \sqrt {x _ {1} ^ {2} + \dots + x _ {n} ^ {2} } $.

Theorem 1) A cubature formula which is invariant under $ G $ possesses the $ m $- property if and only if it is exact for all polynomials of degree at most $ m $ which are invariant under $ G $( see [5]). The method of undetermined coefficients may be defined as the method of constructing invariant cubature formulas possessing the $ m $- property. In the above example, the role of the group $ G $ may be played by the symmetry group of the square. Theorem 1 is of essential importance in the construction of invariant cubature formulas.

For simple domains of integration, such as a cube, a simplex, a ball, or a sphere, and for the weight $ p ( x) = 1 $, one can construct cubature formulas by repeatedly using quadrature formulas. For example, when $ \Omega = K _ {n} = \{ {- 1 \leq x _ {i} \leq 1 } : {i = 1 \dots n } \} $ is the cube, one may use the Gauss quadrature formula with $ k $ knots $ t _ {i} $ and coefficients $ A _ {i} $ to obtain a cubature formula

$$ \int\limits _ {K _ {n} } f ( x) dx \cong \ \sum _ {i _ {1} \dots i _ {n} = 1 } ^ { k } A _ {i _ {1} } \dots A _ {i _ {n} } f ( t _ {i _ {1} } \dots t _ {i _ {n} } ) $$

with $ k ^ {n} $ knots; this is exact for all monomials $ x ^ \alpha $ such that $ 0 \leq \alpha _ {i} \leq 2k - 1 $, $ i = 1 \dots n $, and in particular for all polynomials of degree at most $ 2k - 1 $. The number of knots of such cubature formulas increases rapidly, a fact which limits their applicability.

Throughout the sequel it will be assumed that the weight function is of fixed sign, say

$$ \tag{4 } p ( x) \geq 0 \ \ \mathop{\rm in} \Omega \ \ \textrm{ and } \ \ p _ {1} > 0. $$

The fact that the coefficients of a cubature formula with such a weight function are positive, is a valuable property of the formula.

Theorem 2) If the domain of integration $ \Omega $ is closed and $ p ( x) $ satisfies (4), there exists an interpolatory cubature formula (1) possessing the $ m $- property, $ N \leq \mu $, with positive coefficients and with knots in $ \Omega $. The question of actually constructing such a formula is as yet open.

Theorem 3) If a cubature formula with a weight satisfying (4) has real knots and coefficients and possesses the $ m $- property, then at least $ \lambda = M ( n, l) $ of its coefficients are positive, where $ l = [ m/2] $ is the integer part of $ m/2 $. Under the assumptions of Theorem 3, the number $ \lambda $ is a lower bound for the number of knots:

$$ N \geq \lambda . $$

This inequality remains valid without the assumption that $ x ^ {(} j) $ and $ C _ {j} $ are real.

Regarding cubature formulas with the $ m $- property, one is particularly interested in those having a minimum number of knots. When $ m = 1, 2 $ it is easy to find such formulas for any $ n $, arbitrary $ \Omega $ and $ p ( x) \geq 0 $; the minimum number of knots is precisely the lower bound $ \lambda $: It is equal to 1 in the first case, and to $ n + 1 $ in the second. When $ m \geq 3 $, the minimum number of knots depends on the domain and the weight. For example, if $ m = 3 $, the domain is centrally symmetric, and if $ p ( x) = 1 $, the number of knots is $ 2n $; for a simplex and $ p ( x) = 1 $, it is $ n + 2 $.

By virtue of (4),

$$ \tag{5 } ( \phi , \psi ) = \ I ( \phi \overline \psi \; ) $$

is a scalar product in the space of polynomials. Let $ {\mathcal P} _ {k} $ be the vector space of polynomials of degree $ k $ which are orthogonal in the sense of (5) to all polynomials of degree at most $ k - 1 $. This space has dimension $ M ( n - 1, k) $— the number of monomials of degree $ k $. Polynomials in $ {\mathcal P} _ {k} $ are called orthogonal polynomials for $ \Omega $ and $ p ( x) $.

Theorem 4) There exists a cubature formula (1) possessing the $ ( 2k - 1) $- property and having $ N = M ( n, k - 1) $ knots (the lower bound) if and only if the knots are the common roots of all orthogonal polynomials for $ \Omega $ and $ p ( x) $ of degree $ k $.

Theorem 5) If $ n $ orthogonal polynomials of degree $ k $ have $ k ^ {n} $ finite and pairwise distinct common roots, these roots can be chosen as knots for a cubature formula (1) possessing the $ ( 2k - 1) $- property.

The error of a cubature formula (1) in which $ p ( x) = 1 $ and $ \Omega $ is bounded is defined by

$$ l ( f ) = \ \int\limits _ \Omega f ( x) dx - \sum _ {j = 1 } ^ { N } C _ {j} f ( x ^ {(} j) ). $$

Let $ B $ be a Banach space of functions such that $ l ( f ) $ is a linear functional on $ B $. The norm of the functional $ \| l \| = \sup _ {\| f \| = 1 } l ( f ) $ characterizes the quality of a given cubature formula for all functions of $ B $. Another approach to the construction of cubature formulas is based on minimizing $ \| l \| $ as a function of the knots and the coefficients of the (unknown) cubature formula (with only the number of knots fixed). Implementation of this approach, however, involves difficulties even for $ n = 1 $. Important results for any $ n \geq 2 $ have been obtained by S.L. Sobolev [4]. The question of minimizing $ \| l \| $ as a function of the coefficients for a given set of knots has been solved completely; the problem of choosing the knots is based on the assumption that they form a parallelepipedal grid and that the minimization depends exclusively on the parameters of this grid. The space $ B $, in particular, may be $ L _ {2} ^ {m} ( \mathbf R ^ {n} ) $, where $ m > n/2 $, and in that case the desired cubature formula is assumed to be exact for all polynomials of degree at most $ m - 1 $.

References

[1]	N.M. Krylov, "Approximate calculation of integrals" , Macmillan (1962) (Translated from Russian)
[2]	V.I. Krylov, L.T. Shul'gina, "Handbook on numerical integration" , Moscow (1966) (In Russian)
[3]	A.H. Stroud, "Approximate calculation of multiple integrals" , Prentice-Hall (1971)
[4]	S.L. Sobolev, "Introduction to the theory of cubature formulas" , Moscow (1974) (In Russian)
[5]	S.L. Sobolev, "Formulas for mechanical cubature on the surface of a sphere" Sibirsk. Mat. Zh. , 3 : 5 (1962) pp. 769–796 (In Russian)
[6]	I.P. Mysovskikh, "Interpolatory cubature formulas" , Moscow (1981) (In Russian)

Comments

The polynomial $ {\mathcal L} _ {j} ( x) $" of the influence of the j-th knot" (i.e. defined by $ {\mathcal L} _ {j} ( x ^ {(} i) ) = \delta _ {ij} $) is also called the basic Lagrangian (for $ x ^ {(} j) $).

The "m-property" is also known in Western literature as the degree of precision; i.e. a cubature formula has the $ m $- property if it has degree of precision $ m $.

Reference [a1] is both an excellent introduction as well as an advanced treatment of cubature formulas.

References

[a1]	H. Engels, "Numerical quadrature and cubature" , Acad. Press (1980)
[a2]	P.J. Davis, P. Rabinowitz, "Methods of numerical integration" , Acad. Press (1984)