Cubature formulas of multivariate polynomials arising from symmetric orbit functions

The paper develops applications of symmetric orbit functions, known from irreducible representations of simple Lie groups, in numerical analysis. It is shown that these functions have remarkable properties which yield to cubature formulas, approximating a weighted integral of any function by a weighted finite sum of function values, in connection with any simple Lie group. The cubature formulas are specialized for simple Lie groups of rank two. An optimal approximation of any function by multivariate polynomials arising from symmetric orbit functions is discussed.


Introduction
The purpose of this paper is to extend the results of [24,25], where cubature formulas for numerical integration connected with three types of multivariate Chebyshev-like polynomials arising from Weyl group orbit functions are developed. The specific goal of this article is to derive the cubature rules and corresponding approximation methods for the family of the polynomials arising from symmetric exponential Weyl group orbit sums [2,15] and detail specializations of the general results for the twovariable polynomials.
The family of polynomials induced by the symmetric Weyl group orbit functions (C−functions) forms one of the most natural generalizations of the classical Chebyshev polynomials of one variable -indeed, the lowest symmetric orbit function arising from the Weyl group of A 1 coincides with the common cosine function of one variable and thus induces the family of Chebyshev polynomials of the first kind [27]. The continuous and discrete orthogonality of the sets of cosine functions cos(nx) generalize to the families of multivariate C−functions [13,15,26]. This generalization serves as an essential starting point for deriving the cubature formulas and approximation methods.
Cubature formulas for numerical integration constitute multivariate generalizations of classical quadrature formulas for functions of one variable. A weighted integral over some domain inside R n of any given function is estimated by a finite weighted sum of values of the same function on a specific set of points (nodes). A standard requirement is imposed: the cubature formula has to hold as an exact equality for polynomials up to a certain degree. Numerous types of cubature formulas with diverse shapes of the integration domains and various efficiencies exist [8]. The efficiency of a given cubature formula reflects how the achieved maximal degree of the polynomials relates to the number of the necessary nodes. Optimal cubature formulas of the highest possible efficiency (Gaussian formulas) are for multivariate functions obtained for instance in [23,25,28].
The sequence of Gaussian cubature formulas derived in [23] arises from the antisymmetric orbit functions (S−functions) of the Weyl groups of type A n , n ∈ N. Generalization of these cubature formulas from [23] to polynomials of the S−functions of Weyl groups of any type and rank is achieved in [25]. A crucial concept, which allows the generalization of the A n formulas, is a novel definition of a degree of the underlying polynomials. This generalized degree (m−degree) is based on invariants of the Weyl groups and their corresponding root systems. Besides the polynomials corresponding to the C− and S−functions, two additional families of multivariate polynomials arise from Weyl group orbit functions of mixed symmetries [24]. These hybrid orbit functions (S s − and S l −functions) exist only for root systems of Weyl groups with two different lengths of roots -B n , C n , F 4 and G 2 . The cubature formulas related to the polynomials of these S s − and S l −functions are developed in [24]. Deduction of the remaining cubature formulas, which correspond to the polynomials of the C−functions, completes in this paper the results of [23][24][25]. The integration domains and nodes of these cubature formulas are constructed in a similar way as one-dimensional Gauss-Chebyshev formulas and their Chebyshev nodes.
Instead of a classical one-dimensional interval, the multivariate C−functions are considered in the fundamental domain of the affine Weyl group -a simplex F ⊂ R n . The discrete orthogonality relations of C−functions are performed over a finite fragment of a grid F M ⊂ F , with the parameter M ∈ N controlling the density of F M inside F . The simplex F together with the set of points F M have to be transformed via a transform which induces the corresponding family of polynomials (X−transform). This process results in the integration domain Ω of non-standard shape and the set of nodes Ω M , with specifically distributed points inside Ω. The last ingredient, needed for successful practical implementation, is the explicit form of the weight polynomial K. For practical purposes, the explicit construction of all two-variable cases is presented.
Except for direct numerical integration, one of the most immediate applications of the developed cubature formulas is related multivariate polynomial approximation [4]. The Hilbert basis of the orthogonal multivariate polynomials induced by the C− functions guarantees that any function from the corresponding Hilbert space is expressed as a series involving these polynomials. A specific truncated sum of this expansion provides the best approximation of the function by the polynomials. Among other potential applications of the developed cubature formulas are calculations in fluid flows [9], laser optics [5], stochastic dynamics [33], magnetostatic modeling [34], micromagnetic simulations [6], electromagnetic wave propagation [29], liquid crystal colloids [31] and quantum dynamics [20].
The paper is organized as follows. In Section 2, notation and pertinent properties of Weyl groups, affine Weyl groups and C−functions are reviewed. In Section 3, the cubature formulas related to C−functions are deduced. In Section 4, the explicit cubature formulas of the rank two cases A 2 , C 2 , G 2 are constructed. In Section 5, polynomial approximation methods are developed.

Root systems and polynomials
2.1. Pertinent properties of root systems and weight lattices.
The notation, established in [13], is used. Recall that, to the Lie algebra of the compact, connected, simply connected simple Lie group G of rank n, corresponds the set of simple roots ∆ = (α 1 , . . . , α n ) [1,2,14,32]. The set ∆ spans the Euclidean space R n , with the scalar product denoted by , . The following standard objects related to the set of simple roots ∆ are used: • The marks m 1 , . . . , m n of the highest root ξ ≡ −α 0 = m 1 α 1 + · · · + m n α n .
• n reflections r α , α ∈ ∆ in (n − 1)-dimensional 'mirrors' orthogonal to simple roots intersecting at the origin denoted by r 1 ≡ r α 1 , . . . , r n ≡ r αn . Following [25], we define so called m−degree of λ ∈ P + as the scalar product of λ with the highest dual root η, i.e. by the relation |λ| m = λ, η = λ 1 m ∨ 1 + · · · + λ n m ∨ n . Let us denote a finite subset of the cone of the positive weights P + consisting of the weights of the m−degree not exceeding M by P + M , i.e. P + M = λ ∈ P + | |λ| m ≤ M . Recall also the separation lemma which asserts for λ ∈ P + , λ = 0 and any M ∈ N that Note that this lemma is proved in [25] for M > m only -the proof, however, can be repeated verbatim with any M ∈ N. For two dominant weights λ, ν ∈ P + for which ν ≤ λ we have for their m−degrees Taking into account equation (2), we have the following proposition.

Affine Weyl groups.
The Weyl group W is generated by n reflections r 1 , . . . , r n and its order |W | can be calculated using the formula |W | = n! m 1 . . . m n c. (6) The affine Weyl group W aff is the semidirect product of the Abelian group of translations Q ∨ and of the Weyl group W , The fundamental domain F of W aff , which consists of precisely one point of each W aff -orbit, is the convex hull of the points 0, ω ∨ 1 m 1 , . . . , ω ∨ n mn . Considering n + 1 real parameters y 0 , . . . , y n ≥ 0, we have F = y 1 ω ∨ 1 + · · · + y n ω ∨ n | y 0 + y 1 m 1 + · · · + y n m n = 1 .
The volumes vol(F ) ≡ |F | of the simplices F are calculated in [13].
Considering the standard action of W on R n , we denote for λ ∈ R n the isotropy group and its order by Stab(λ) = {w ∈ W | wλ = λ} , h λ ≡ |Stab(λ)|, and denote the orbit by W λ = {wλ ∈ R n | w ∈ W } . Then the orbit-stabilizer theorem gives for the orders Considering the standard action of W on the torus R n /Q ∨ , we denote for x ∈ R n /Q ∨ the order of its orbit by ε(x), i.e.
For an arbitrary M ∈ N, the grid F M is given as cosets from the W −invariant group 1 M P ∨ /Q ∨ with a representative element in the fundamental domain F The representative points of F M can be explicitly written as The numbers of elements of F M , denoted by |F M |, are also calculated in [13] for all simple Lie algebras.

Orbit functions.
Symmetric orbit functions [16] are defined as complex functions C λ : R n → C with the labels λ ∈ P + , Note that in [13] the results for C−functions are formulated for the normalized C−functions Φ λ which are related to the orbit sums (12) as Φ λ = h λ C λ . Due to the symmetries with respect to the Weyl group W as well as with respect to the shifts from Q ∨ it is sufficient to consider C−functions restricted to the fundamental domain of the affine Weyl group F . Moreover, the C−functions are continuously orthogonal on F , and form a Hilbert basis of the space L 2 (F ) [16], i.e. any function f ∈ L 2 (F ) can be expanded into the series of C−functions Special case of the orthogonality relations (14) is when one of the weights is equal to zero, For any M ∈ N , the C−functions from a certain subset of P + are also discretely orthogonal on F M and form a basis of the space of discretized functions C F M of dimension |F M | [13]; special case of these orthogonality relations is when one of the weights is equal to zero modulo the lattice M Q, The key point in developing the cubature formulas is comparison of formulas (16) and (17) in the following proposition.
Proposition 2.2. For any M ∈ N and λ ∈ P + 2M −1 it holds that Proof. Suppose first that λ = 0. Then from (16) and (17) we obtain Secondly let λ = 0 and |λ| m < 2M . The from the separation lemma (4) we have that λ / ∈ M Q and thus Let us denote for convenience the C−functions corresponding to the basic dominant weights ω j by Z j , i.e.
Recall from [2], Ch. VI, §4 that any W −invariant sum of the exponential functions e 2πi ν, a can be expresssed as a linear combination of some functions C λ with λ ∈ P + . Also for any λ ∈ P + a function of the monomial type Z λ 1 1 Z λ 2 2 . . . Z λn n can be expressed as the sum of C−functions by less or equal dominant weights than λ, i.e.
Conversely, any function C λ , λ ∈ P + can be expressed as a polynomial in variables Z 1 , . . . , Z n , i.e. there exist a multivariate polynomials p λ ∈ C[y 1 , . . . , y n ] such that Antisymmetric orbit functions [16] are defined as complex functions S λ : R n → C with the labels λ ∈ P ++ , The antisymmetry with respect to the Weyl group W and the symmetry with respect to the shifts from Recall that Proposition 9 in [16] states that for the lowest S−function S it holds that Since the square of the absolute value |S | 2 = S S is a W −invariant sum of exponentials it can be expressed as a linear combinations of C−functions. Each C−function in this combination is moreover a polynomial of the form (20). Thus there exist a unique polynomial K ∈ C[y 1 , . . . , y n ] such that 3. Cubature formulas The key component in the development of the cubature formulas is the integration by substitution. The X−transform transforms the fundamental F ⊂ R n domain to the domain Ω ⊂ R n on which are the cubature rules defined. In order to obtain a real valued transform we first need to examine the values of the C−functions.
The C−functions of the algebras are real-valued [15]. Using the notation (3), for the remaining cases it holds that Specializing the relations (26) for the C−functions corresponding to the basic dominant weights Z j , we obtain that the functions Z j are real valued, except for the following cases for which it holds that Taking into account (27), we introduce the real-valued functions X j , j ∈ {1, . . . , n} as follows. For the cases (25) we set and for the remaining cases (27) we define Thus, we obtain a crucial mapping X : R n → R n given by The image Ω ⊂ R n of the fundamental domain F under the mapping X forms the integration domain on which the cubature rules will be formulated, i.e.
In order to use the mapping X for an integration by substitution we need to know that it is one-to-one except for possibly some set of zero measure. Since the image Ω M ⊂ R n of the set of points F M under the mapping X forms the set of nodes for the cubature rules, i.e.
a discretized version of the one-to-one correspondence of the restriced mapping X M of X to F M , i.e.
is also essential. Note that due to the periodicity of C−functions (13), the restriction (33) is well-defined for the cosets from F M .
Proposition 3.1. The mapping X : F → Ω, given by (30), is one-to-one correspondence except for some set of zero measure. For any M ∈ N is the restriction mapping X M : F M → Ω M , given by (33), one-to-one correspondence and thus it holds that Proof. Let us assume that there exists a set F ⊂ F of non-zero measure such that X(x) = X(y) with x, y ∈ F . Since the transforms (28), (29) are as regular linear mappings one-to-one correspondences, this fact implies that Z 1 (x) = Z 1 (y), . . . , Z n (x) = Z n (y) with x, y ∈ F . Then from the polynomial expression (20) we obtain for all λ ∈ P + that it holds that C λ (x) = C λ (y). Since the C−functions C λ , λ ∈ P + form a Hilbert basis of the space L 2 (F ) we conclude that for any f ∈ L 2 (F ) is valid that Retracing the steps of the continuous case above, let us assume that there exist two distinct points x, y ∈ F M , x = y such that X(x) = X(y). Since the transforms (28), (29) are as regular linear mappings one-to-one correspondences, this fact again implies that Z 1 (x) = Z 1 (y), . . . , Z n (x) = Z n (y). Then from the polynomial expression (20) we obtain for all λ ∈ P + that it holds that C λ (x) = C λ (y). The same equality has hold for those C−functions C λ which form a basis of the space C F M . We conclude that for The absolute value of the determinant of the Jacobian matrix of the X transform (30) is essential for construction of the cubature formulas -its value is determined in the following proposition.
where κ is defined as otherwise. (36) Proof. Note that the X transform can be composed of the the following two transforms: the transform ζ : x → (Z 1 (x), . . . , Z n (x)) and the transform R : (Z 1 , . . . , Z n ) → (X 1 , . . . , X n ) via relations (28), (29). To calculate the Jacobian of the transform ζ, let us denote by α ∨ the matrix of the coordinates (in columns) of the vectors α ∨ 1 , . . . , α ∨ n in the standard orthonormal basis of R n and by a 1 , . . . , a n the coordinates of a point x ∈ R n in α ∨ −basis, i.e. x = a 1 α ∨ 1 + · · · + a n α ∨ n . If a denotes the coordinates a 1 , . . . , a n arranged in a column vector then it holds that x = α ∨ a. The absolute value of the Jacobian of the mapping a → (Z 1 (α ∨ a), . . . , Z n (α ∨ a)) is according to equation (32) in [25] given by (2π) n |S (α ∨ a)|. Using the chain rule, this implies for the absolute value of the Jacobian |J x (ζ)| of the map ζ that It can be seen directly from formula (6) and Proposition 2.1 in [13] that The calculation of the absolute value of the Jacobian determinant κ = |J x (R)| is straightforward from definitions (28), (29).

The cubature formula.
We attach to any λ ∈ P + a monomial y λ ≡ y λ In order to investigate how the m−degree of a polynomial changes under the substitution of the type (29) we formulate the following proposition.
Proof. Since any polynomial p ∈ Π M is a linear combination of monomials y λ with |λ| m ≤ M , it is sufficient to prove p ∈ Π M for all monomials p of m−degree at most M . If p is a monomial y λ with |λ| m ≤ M , then Using the binomial expansion, we obtain Therefore, the m−degree of the polynomial p is given by Since we assume that m ∨ j = m ∨ k , we conclude that deg m p = n l=1 λ l m ∨ l = deg m p.
Having the X M transform (33), it is possible to transfer uniquely the values (10) of ε(x), x ∈ F M to the points of Ω M , i.e. by the relation Taking the inverse transforms of (28), (29) and substituting them into the polynomials (24), (20) we obtain the polynomials K, p λ ∈ C[y 1 , . . . , y n ] such that Theorem 3.4 (Cubature formula). For any M ∈ N and any p ∈ Π 2M −1 it holds that Proof. Proposition 3.1 guarantees that the X−transform is one-to-one except for some set of measure zero and Proposition 3.2 together with (23) gives that the Jacobian determinant is non-zero except for the boundary of F . Thus using the integration by substitution y = X(x) we obtain Ω p(y)K − 1 2 (y) dy = κ(2π) n |F ||W | F p(X(x)) dx.
Note that for practical purposes it may be more convenient to use the cubature formula (41) in its less developed form resulting from (42), This form may be more practical since the explicit inverse transform to X M is usually not available and, on the contrary, the calculation of the coefficients ε(x) and the points X(x), x ∈ F M is straightforward.
The values of ε j can be found in Table 2. The cubature rule for any p ∈ Π 2M −1 is of the form Ω p(y 1 , y 2 )K − 1 2 (y 1 , y 2 ) dy 1 dy 2 = π 2 9M 2 j∈I M ε j p(y The cubature rule (50) is an analogue of the formula deduced in [22] using the generalized cosine functions T C k and generalized sine functions T S k defined by It can be shown by performing the following change of variables and parameters that generalized cosine and sine functions actually coincide (up to scalar multiplication) with the symmetric C-functions and antisymmetric S-functions of A 2 . More precisely, we obtain  Table 3 shows the values of the finite weighted sums for M = 10, 20, 30, 50, 100.
The integral domain Ω can be described explicitly as The weight K-polynomial (40) is given explicitly as K(y 1 , y 2 ) = (y 2 1 − 4y 2 )((y 2 + 4) 2 − 4y 2 1 ) . The index set which will label the sets of points F M and Ω M is introduced via Thus the grid F M consists of the points If for j = [s 0 , s 1 , s 2 ] ∈ I M we denote (y then the set of nodes Ω M consists of the points The integration domain Ω is together with the set of nodes Ω 15 depicted in Figure 2. Similarly to the case A 2 , each point of F M as well as of Ω M is labeled by the index set I M and it is convenient for the point x ∈ F M and its image in Ω M labeled by j ∈ I M to denote The values of ε j can be found in Table 2. The cubature rule for any p ∈ Π 2M −1 takes the form Ω p(y 1 , y 2 )K − 1 2 (y 1 , y 2 )dy 1 dy 2 = π 2 4M 2 j∈I M ε j p(y 2 )) .  Table 2. The values of ε j are shown for A 2 , C 2 and G 2 with denoting the corresponding coordinate different from 0.  Table 3 shows the values of the finite weighted sums for M = 10, 20, 30, 50, 100.
Performing the transform (28), the resulting real-valued functions X 1 , X 2 are given by X 1 = 2(cos 2πa 1 + cos 2π(a 1 − 3a 2 ) + cos 2π(2a 1 − 3a 2 )) , The integral domain Ω can be described explicitly as Ω = (y 1 , y 2 ) ∈ R 2 | − 2((y 2 + 3) The weight K-polynomial (40) is given explicitly as K(y 1 , y 2 ) = (y 2 2 − 4y 1 − 12)(y 2 1 − 4y 3 2 + 12y 1 y 2 + 24y 1 + 36y 2 + 36) . The index set which will label the sets of points F M and Ω M is introduced via Thus the grid F M consists of the points The values of ε j can be found in Table 2. The cubature rule for any p ∈ Π 2M −1 is of the form As in the case A 2 , the cubature rule (52) is similar to the Gauss-Lobatto cubature formula derived in [21], where Xu et al study four types of functions CC k (t), SC k (t), CS k (t) and SS k (t) closely related to the orbit functions over Weyl groups. Concretely, the functions CC k (t) and SS k (t) are given by where the variable t is given by homogeneous coordinates, i.e.,  Table 3. It shows the estimations of the integrals of 1 over the regions Ω of A 2 , C 2 and G 2 resp. by finite weighted sums of the right-hand side of (50),(51) and (52) resp. for M = 10, 20, 30, 50, 100.
H . It can be verified that the linear transformations of variable and parameter give the connection with C−functions and S−functions of G 2 :  Table 3 shows the values of the finite weighted sums for M = 10, 20, 30, 50, 100.
Since the polynomial function (40) is continuous and strictly positive in Ω • , its square root K − 1 2 can serve as a weight function for the weighted Hilbert space L 2 K (Ω), i.e. a space of complex-valued cosets of measurable functions f such that Ω |f | 2 K − 1 2 < ∞ with an inner product defined by Our aim is to consctruct a suitable Hilbert basis of L 2 K (Ω). Taking the inverse transforms of (28), (29) and substituting them into the polynomials (20) we obtain the polynomials p λ ∈ C[y 1 , . . . , y n ] such that Moreover, successively applying Proposition 3.3 to perform the substitutions (29) in p λ and taking into account (37) we obtain that deg m p λ = |λ| m .
(55) Calculating the scalar product (53) for the p λ polynomials (54), we obtain that the continuous orthogonality of the C−functions (14) is inherited, i.e.
Assigning to any function f ∈ L 2 K (Ω) a function f ∈ L 2 (F ) by the relation f (x) = f (X(x)) and taking into account the expansion (15) we obtain for its expansion coefficients c λ that

Concluding remarks
• Due to the generality of the present construction of the cubature formulas, some of the cases presented in this paper appeared already in the literature. The case A 2 is closely related to twovariable analogues of Jacobi polynomials on Steiner's hypocycloid [19]; the case C 2 is related to two-variable analogues of Jacobi polynomials on a domain bounded by two lines and a parabola [17][18][19] and the corresponding Gaussian cubature formulas induced by these polynomials are studied for example in [28]. The case G 2 and its cubature formulas are detailed in [21]. • The Chebyshev polynomials of the first kind induce the cubature formula of the maximal efficiency -Gauss-Chebyshev quadrature [12,27]. The nodes of this formula are M roots of the Chebyshev polynomials of the first kind of degree M and it exactly evaluates a weighted integral for any polynomial of degree at most 2M − 1. The set of nodes Ω M of A 1 does not correspond to the set of roots of the Gauss-Chebyshev quadrature -Ω M of A 1 consists of M + 1 points and includes two boundary points of the interval. Thus, the number of points Ω M exceeds by one the minimum number of nodes and the resulting cubature formula is not optimal. This phenomenon, already observed for the A n sequence in [23], generalizes to all cases of Weyl groups. Even though the developed formulas are slightly less efficient than the optimal ones, they contain points on the entire boundary of Ω M -this might be useful for some applications. • The technique for construction of the cubature formulas in this article is based on discrete orthogonality of C−functions over the set F M . For the case A 1 the cosine transform corresponding to the discrete orthogonality of cosines is standardly known as discrete cosine transform of the type DCT-I [3]. The M roots of the Chebyshev polynomials of the first kind, which lead to the optimal cubature formula, enter the discrete orthogonality relations of the type DCT-II. This indicates that the discrete orthogonality relations of C-functions, which generalize the transforms of the type DCT-II, DCT-III and DCT-IV for admissible cases of Weyl groups [10], might lead to cubature formulas of higher efficiency. poses an open problem.
• The Clenshaw-Curtis method for deriving cubature formulas is based on expressing a given function into a series of the corresponding set of orthogonal polynomials and then integrating the series term by term [7,30]. The coefficients in the series are calculated using discrete orthogonality properties. The comparison of the resulting cubature formulas by the Clenshaw-Curtis technique and the method used in the present paper and in [24,25] deserves further study. • The polynomials of the orbit functions used in the present paper and in [24,25] to derive the cubature formulas are directly related to the Jacobi polynomials associated to root systems as well as to the Macdonald polynomials. The discrete orthogonality of the Macdonald polynomials with unitary parameters, achieved in [11], opens a possibility of an extension of the current cubature formulas to the corresponding subset of the Macdonald polynomials.