Construction of Multivariate Interpolation Hermite Polynomials for Finite Element Method

A new algorithm for constructing multivariate interpolation Hermite polynomials in analytical form in a multidimensional hypercube is presented. These polynomials are determined from a specially constructed set of values of the polynomials themselves and their partial derivatives with continuous derivatives up to a given order on the boundaries of the finite elements. The efficiency of the finite element schemes, algorithms and programs is demonstrated by solving the Helmholtz problem for a cube.


Introduction
In recent decade we presented a new algorithm for the calculation of high-order Lagrange and Hermite interpolation polynomials (LIP and HIP) [1] of the simplex in analytical form, their classification and a typical example of the triangle element for high-accuracy finite element method (FEM) schemes [2][3][4]. We have illustrated the efficiency of the schemes using high-order accuracy LIP and HIP on benchmark calculations of exactly solvable boundary value problems(BVPs)for a triangle membrane, a hypercube and a helium atom [5,6].
In this paper we present a new symbolic algorithm implemented in Maple for constructing Hermitian finite elements or piece-wise multivariate Birkhoff interpolants in a standard d-dimensional cube that generalizes the construction and algorithm proposed for a three and four dimensional cube [7][8][9]. Our construction does not use explicit inverse matrices in the solution of the algebraic problem for a set of algebraic equations with respect to the unknown coefficients of polynomials of d variables [10]. The algorithm yields explicit expressions in analytical form for the Hermite interpolation polynomials (HIPs). The basis functions of the finite elements are high-order polynomials determined from a specially constructed set of values of the polynomials themselves and their partial derivatives up to a given order at the vertices of the hypercube. Such a choice of values allows us to construct a piecewise polynomial basis continuous at the boundaries of the finite elements together with the derivatives up to a given order. In the case of a d-dimensional cube, it is shown that the basis functions are determined by products of d HIPs depending on each of the d variables given by means of the recurrence relations in analytical form [2]. Using this fact we propose a new symbolic algorithm implemented in Maple for analytical calculation of the basis functions, i.e. HIPs of d variables with the derivatives continuous at the faces of the standard d hypercube. The method can be used to solve elliptic BVPs by means of a high-accuracy FEM. Its advantages are the reduced computational cost and the availability of accurate derivatives of the function interpolated. The present study is motivated by possible applications of the FEM on simplexes [4] or hypercubes [10] for solving a BVP in the collective nuclear model with tetrahedral symmetry [11], as well as by other applications in different areas, e.g., flow dynamics in unsteady fluid systems [12,13] and molecular dynamics [14].

Formulation of the problem
Consider a self-adjoint boundary-value problem (BVP) for the elliptic differential equation For the principal part coefficients of Eq. (1), the condition of uniform ellipticity holds in the bounded and V(x) are real-valued functions, continuous together with their generalized derivatives up to a given order in the domain x ∈Ω = Ω ∪ ∂Ω ∈ R d with the piecewise continuous boundary S = ∂Ω, which provides the existence of nontrivial solutions Φ(x) obeying the boundary conditions (BCs) of the first (I) or the second (II) kind [4]. For a discrete spectrum problem, the functions Φ m (x) from the Sobolev space H 2 2 (Ω), Φ m (x) ∈ H 2 2 (Ω), corresponding to the real eigenvalues E: E 1 ≤ E 2 ≤ · · · ≤ E m ≤ · · · obey the orthonormalization conditions [4].
The polyhedral domainΩ ⊂ R d is decomposedΩ =Ω h (x) = Q q=1 ∆ q in finite elements in the form of d dimensional simplexes or hypercubes ∆ q with HIPs or LIPs ϕ κ r (x), x ∈ R d , calculated using the algorithms of Refs. [4,10]. The piecewise polynomial functions P l (x) ∈ C κ of order p with continuous derivatives up to order κ ≤ κ max −1 are constructed by joining the polynomials, ϕ κ r (x), on the finite elements ∆ q ∈Ω h (x). The expansion of the sought solution Φ h m (x) from the Sobolev space H s≥1 2 (Ω) in the basis of piecewise polynomial functions P l (x), Φ h m (x) = L l=1 P l (x)Φ h lm leads to the generalized algebraic eigenvalue problem (AEP), (A − BE h m )Φ h m = 0, solved using the standard method [15]. The elements of the symmetric matrices of stiffness A and mass B comprise the integrals, which are calculated on the elements ∆ q in the domainΩ =Ω h (x) = Q q=1 ∆ q , recalculated into the local coordinates on the element ∆. The deviation of the approximate solution h is the maximal size of the finite element ∆ q , p is the order of the FEM scheme, c 1 and c 2 are coefficients independent of h.

Algorithm for constructing Hermitian finite elements in hypercube
The HIPs ϕ κ r (x) ≡ ϕ κ 1 ···κ i ···κ d r 1 [3], where the multiplicity of the degeneracy is shown in square brackets. Results of FEM are marked: for cubic elements with HIPs of the ninth order p = p 1 p 2 p 3 = 9 by M, and for tetrahedral elements with LIPs of the third order by T3, fourth order by T4, and fifth order by T5.
x r 1 ···r i ···r d = (x 1r 1 , . . . , x ir i , . . . , x dr d ), x ir i = ((p − r i )x i;min + r i x i;max )/p; r i = 0, . . . , p, i = 1, . . . , d are determined by relations [1,7] As an example of applying the above algorithms, we present the results of solving the Helmholtz problem (1) with BC (II) for a 3D cube with the edge length π. Figure 1 shows the discrepancy δE m = E h m − E m , m = 0, 1, . . . , 11 (except δE 0 ) of the numerical eigenvalues E h m of this problem from their exact values. The results were calculated using FEM with HIPs of the ninth order (M (n = 4, L = 1000)) for cubic elements and LIPs [3] of the third order (T3 (n = 3, L=1000)), fourth order (T4 (n = 2, L = 729)), and fifth order (T5 (n = 2, L = 1331)) for tetrahedral elements. The number of parts n into which we divide the edge of the cube and the length of the eigenvector L = (pn + 1) d κ d max of AEP are shown in parentheses. One can see that the accuracy of the FEM with HIPs of the ninth order p = p 1 p 2 p 3 = 9 with continuous first derivatives κ max − 1 = 1 at each face of the 3-dimensional cube is higher than with the LIPs of the third order, but lower than with the LIPs of the fourth and fifth order.

Conclusions
In this paper we present a new symbolic algorithm implemented in Maple for constructing in analytical form the Hermitian finite elements with d variables in a standard d-dimensional hypercube. In the construction of this algorithm we used the fact that the basis functions are determined by a product of d one-dimensional HIPs depending on each of the d variables with continuous derivatives up to a given order κ on the boundaries of the the finite elements using the recurrence relations (3). It allows us to construct a piecewise polynomial basis continuous on the boundaries of finite elements together with the derivatives up to the given order, which can be used to solve elliptic BVPs as well as other problems with partial derivatives of a high order by means of the high-accuracy finite element method.