# 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)

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$.