Equivalent Finite Element Formulations for the Calculation of Eigenvalues Using Higher-Order Polynomials

This paper investigates higher-order approximations in order to extract Sturm-Liouville eigenvalues in one-dimensional vibration problems in continuum mechanics. Several alternative global approximations of polynomial form such as Lagrange, Bernstein, Legendre as well as Chebyshev of first and second kind are discussed. In an instructive way, closed form analytical formulas are derived for the stiffness and mass matrices up to the quartic degree. A rigorous proof for the transformation of the matrices, when the basis changes, is given. Also, a theoretical explanation is provided for the fact that all the aforementioned alternative pairs of matrices lead to identical eigenvalues. The theory is sustained by one numerical example under three types of boundary conditions.


Introduction
In the framework of the standard Galerkin/Ritz method for the solution of one-dimensional boundary value problems governed by a differential equation within the domain [0,L], the usual procedure consists of subdividing [0,L] into a certain number of finite elements for which piecewise-linear (i.e., local) interpolation is assumed [1]. In general, the shortest the elements are the more accurate the numerical solution is (h-version). Alternatively, higher order p-methods [2] suggest the introduction of nodeless basis functions based on differences of Legendre polynomials (up to the seventh degree) that cooperate with the two linear shape functions, i.e., , the latter associated to the ends x = 0 and x = L. A literature survey suggests that for a certain discretization of the domain, the corresponding p-version is generally more accurate than the h-version [3][4][5][6].
The matter of using higher order approximations through computer-aided-design (CAD) based Coons-Gordon macroelements has been recently discussed for two-and three-dimensional problems [7][8][9][10][11]. In those works some similarities and differences of the so-called 'Coons macroelements' with respect to the 'higher order p-method' have been reported in detail. Moreover, alternative CAD based NURBS or/and Bézier techniques have been proposed within the last eighteen years [12][13][14][15][16].
In this context, this paper continues the investigation on the eigenvalue problem by moving from 2-D and 3-D to 1-D eigenvalue problems, and seeks for any similarities or essential differences between five alternative methods. The study includes classical Lagrange polynomials and extends to the Bernstein polynomials that are inherent in the definition of Bézier CAD curves [17], as well as to Chebyshev polynomials that have been previously used in spectral and collocation methods (e.g. [18]). In this paper it was found that all the aforementioned polynomials are equivalent in the sense that (after the proper transformation) they symbolically coincide with the classical higher order p-method (or p-version) [2] as well as with the class { } n x (Taylor series).

General
A general Sturm-Liouville problem can be written in the following differential equation It can be reduced to a study of the canonical Liouville normal form Without loss of generality, in this paper we deal with the particular case that ( ) 0 q x = , for which Eq(1) degenerates to the well-known 'Helmholtz equation': Equation (2) covers many practical problems in physics and engineering such as the axial elastic vibrations of a beam or the sound propagation along a straight acoustic pipe.
The solution of Eq(2) is usually expressed as a series expansion in two alternative ways: where f i (x) are the basis functions and a j the generalized co-ordinates that refer to n+1 nodeless parameters. Also, N j (x) are the shape functions and U j the nodal displacements at the positions x=x j , j=0,…,n. It is reminded that N j (x) are cardinal functions, i.e. N j (X i )=δ ij (= Kronecker's delta) and also partition the unity, i.e.

( )
. There are also some formulations in which both shape and basis functions participate such as higher order p-methods [2]. In any case, the Galerkin/Ritz procedure [1] consists first of the Galerkin method, i.e.
denoting the weighting function, and second by the partial integration of the second order term , which finally leads to the alternative matrix formulations: where the elements of the mass ( , M M ) and stiffness ( , K K ) matrices are given by: According to the standard literature [1], for a finite element of length l k , the two linear shape functions, that is

Basis Functions and Shape Functions
, are associated to the ends x = 0 and x = l k . It is well known [1] that the basis functions related to the aforementioned shape functions are: giving unity at k x and passing through n points. Typical sets of Lagrange polynomials are shown in Table 1. Therefore, Equation (3a) holds with denoting the jth Lagrange polynomial ( ). As for each k the denominator in Eq(6) is a constant, whereas the nominator is a polynomial of degree n, the latter can be written in terms of its roots ( 1   0  2  1  1  2  1  1 , , , , , , follows:       According to standard computer-aided-geometry knowledge, for example [17], a Bézier curve of n-th degree is defined in terms of the normalized coordinate u x L = as In standard mathematical texts the variable (7) is called "Bernstein polynomial of degree n" and is denoted by [19, p. 36 The geometrical coefficients, { } i P , are called control points. Equation (8) denotes that for a curve that is defined by (n + 1) control points, the highest degree is , which means that the degree of a Bézier curve is defined by the number of control points. It has been shown that Bernstein polynomials have the property of partition of unity and also the first ( ) and the last of them ( ) give unity at the first ( of the interval [0,1], respectively. It is remarkable that at any internal control point is less than the unity [17], and also vanishes at the endpoints u = 0, 1. It is trivial to prove that, at least for the case of a straight segment, when the interval [0,L] is uniformly divided into a number of so-called breakpoints, the control points coincide with them. Following the ideas of CAD/CAE integration previously proposed by the author in the context of closely related isoparametric Coons interpolation [7][8][9][10][11] as well as the "isogeometric" idea of Hughes et al. [14] referring to NURBS representation, in this work we propose to substitute the Cartesian coordinate ( )

Higher Order P-Method
Following Szabó and Babuška [2], the variable is expanded to a series like that in Eq(3b). According to this technique [2, p.38], the two first basis functions, which correspond to x=0 and x=L, are identical with the linear shape functions associated to the end points of the entire interval [0,L], i.e., are called nodal shape functions. The rest (n-1) basis functions are called internal shape functions and are taken as suitable integrals of Legendre polynomials properly multiplied by specific coefficients dependent on the ascending order. It must become clear that the latter refer to nodeless quantities (general coefficients) and, therefore, they are also called "bubble modes".

Chebyshev Polynomials
Chebyshev polynomials are categorized as first and second kind.

General Remarks
From the above analysis, it becomes obvious that Lagrange, Bernstein, Legendre, and Chebyshev polynomials are different options with the following characteristics: I. Lagrange polynomials are cardinal shape functions that operate directly on the nodal values, IV. Chebyshev polynomials (both of first and second kind) differ from the above cases in the sense that they neither become unit nor vanish at the ends of the domain.
A careful inspection of all above five polynomial interpolations [Eqs. (6), (8), (9) and (10) is easily concluded that all these polynomials can be expanded in a Taylor (power) series of the form: For each of the five abovementioned polynomials, different coefficients j a [Eq (11)] are derived.

Special Cases
In order to bring insight in the topic of subsection 2.2.6, some special cases of small size will be discussed in full detail.
I. First, in the first four (Lagrange, Bernstein, P-method, Chebyshev of 1 st kind) alternative types of polynomials, the case of using only one segment for the discretization (n = 1) is strictly related to the linear interpolation between the two endpoints. Therefore, the same classical finite element matrices, which are given below by Eq.(12a) and (12b), are obtained. For Chebyshev of 2 nd kind ( ( ) ( ) the factor "2" in induces a slight modification.
II. Second, in the case of two subdivisions, equivalent polynomials which include terms up to the second degree, is represented in Table 2. From Table 2, it becomes obvious that all functional sets are complete, as they include all the terms of the set {1, u, u 2 }. Moreover, while Lagrange and Bernstein (basis) polynomials ensure the partition of unity (rigid body) property, the p-method and Chebyshev polynomials do not.

Analytically Calculated Matrices
Based on Eqs(5), the stiffness and mass matrices are analytically calculated through manipulation and are shown below.

Linear Interpolation (n = 1)
In this case, all three formulations of this paper degenerate to the simple case of the conventional linear finite element of the literature [1]. The matrices are given by: where denotes the length of the finite element. It is well known that accuracy increases when the total length L of the domain is divided into an increasing number of such elements (h-version). In this case, the total mass and stiffness matrices are obtained from the assemblage of all finite elements involved [1].

Cubic Approximation (n = 3)
Based on Eqs (5), the stiffness and mass matrices are given by: 1. Lagrange polynomials: 3. Higher order P-method: One can notice that in the case of the p-method, no symmetry exists along the diagonal terms of the mass matrix which correspond to internal modes (

Quartic Approximation (n = 4)
Based on Eqs (5) 3. Higher order P-method: One can notice that in the case of the p-method, again no symmetry exists along the diagonal terms of the mass matrix which correspond to internal modes ( 2 5 2 45 ≠ ). Table 2. Alternative polynomials involved in case of n = 2 (series expansion is included in the lower row). 2 u u + Table 3. Alternative polynomials involved in case of n = 3.

Implementation of the Boundary Conditions
In this two-point boundary value problem, three alternative types of boundary conditions may exist, as follows: i). Free-free (F-F) boundary conditions: ii). Both ends under Dirichlet (D-D) boundary conditions: iii). One end under Dirichlet, the other free (D-F): In general, the mathematical problem becomes: det 0 It is well known that for small problems Eq(17) can be solved even through the characteristic polynomial, whereas for large-scale problems it is usually elaborated using any known algorithm such as subspace iteration, Lanczos, QR, etc.
Details about the treatment of boundary conditions will be discussed in Section 4.

Numerical Results
The domain [0, L] is subjected to either of boundary conditions (b.c.) given by Equations (16)  All numerical and symbolic (Symbolic Toolbox) calculations were performed on a PC (DELL Latitude E6500) using MATLAB 7.12.0 (R2011a). Given the stiffness K and mass matrices M, the symbolic eigenvalues were found using the command eig(inv(M)*K), whereas the alternative command eig(K,M) works for the numerical operations only.
For a given polynomial of degree n, the error norm of the i-th calculated eigenvalue is determined as:

Free-Free (F-F) Boundary Conditions
This case requires no special care as none row or column has to be eliminated.
In all six cases, i.e. (i) Lagrange, (ii) Bernstein, (iii) p-method, (iv) Chebyshev of 1st kind, (v) Chebyshev of 2nd kind, and finally (vi) Taylor series, the eigenvalues were found to be identical. Moreover, the convergence quality is excellent (the same quality with the results pre-sented in Table 4) but it is not presented for the sake of briefness.

One Dirichlet (D-F) or Two Dirichlet (D-D) Boundary Conditions
For the last two types of boundary conditions [D-D and D-F: Eq(16)b,c], and for the first three types of polynomials [i.e. (i) Lagrange, (ii) Bernstein, and (iii) p-method], the row(s) and column(s) that correspond to the restricted end(s) is (are) simply eliminated. As previously happened, in all these three types of polynomials the eigenvalues were found to be identical. Moreover, convergence quality is shown in Table 4, labeled as D-D and D-F.
In contrast, in case of either Chebyshev polynomials or Taylor series the way of elimination of row(s) and column(s) is not apparent yet.

A Theoretical Explanation
Since all types of polynomials dealt in this study [i.e. (i) Lagrange, (ii) Bernstein, (iii) p-method and relevant Legendre polynomials, as well as (iv) Chebyshev of first and (v) Chebyshev of second kind] span the basis functions included in the class { } n x (Taylor power series): it becomes necessary to obtain the corresponding stiffness and mass t matrices. It is clarified that henceforth the upper left superscript 't' in Eq(20) will stand for the word "Taylor".

ERRORS (in %)
Degree of polynomial n = 2 n = 3 n = 4 n =5 n = 6 n = 7 n = 8 n = 9 n = 10 n = 11 n = 12 In this case the variable U is expanded in a Taylor series: According to Eq.(5b), the elements of the mass and stiffness matrices will be given by:

A Proof for the Identical Eigenvalues
For a single square matrix of order (n+1), Ralston and Rabinowitz [19, p. 484] (Theorem 10.3) have demonstrated that any similarity transformation A 1 − PAP applied to leaves the eigenvalues of the matrix unchanged. The same is found in the textbook of Bathe [20, p. 45] where, provided is orthogonal, a poof is given through the characteristic polynomial.
Based on Eq(23), the corresponding mass matrices [cf. Eq(5b)] are given in compact form as follows: Second, we consider the equivalent expressions of U(x) using either Taylor series or the μ -type polynomial basis, for which it holds: Since U(x) is a scalar, it equals to its transpose, therefore: Left-multiplication of Eq(25) by Eq(24) in parts, leads to: Substituting Eq(24) into Eq(28), one obtains the desired matrix identity: Equation (29) constitutes the proof of the theorem. Remark: It is reminded that Eq(29) is the well known matrix transformation (from local to global system) in the finite element praxis (e.g., [1,20]). Nevertheless, in this study only mathematical considerations have been taken into account, whereas in engineering books it is usually derived on the basis of virtual work (the latter is considered to be invariable in both local and global orthogonal co-ordinate systems). In contrast, in this work no orthogonality relationship between the different bases was considered.
THEOREM-2. The contents of Theorem-1 is extended to the stiffness matrices and t as well. μ K K PROOF. Taking the derivatives in x of the vectors in Eq(23), we obtain the two new vectors: Therefore, the corresponding stiffness matrix [Eq(5b)] is written as: Henceforth the proof is the same as in Theorem-1, leading to the identity: Therefore, we derive a relationship in the form: which in vector form is written as: or equivalently: Substituting Eq(41) into Eq(35) we derive: ( )

Discussion
The motivation of this work was the interesting paper by Çelik [18] that deals with the collocation method using Chebyshev polynomials, which constitute a complete functional set. Previously, the state-of-the-art was the use of basis functions in the form of algebraic polynomials in the form  [21]). At that period the author had accumulated numerical experience on the excellent behavior of global approximation, either using B-splines or Lagrange polynomials ( [7][8][9][10][11], among others). Later, when he tried to compare the eigenvalues obtained using (i) Lagrange polynomials, (ii) Bernstein (basis) ones, and (iii) Taylor series, a numerical coincidence was remarked when the same collocation points were used [22].
The abovementioned numerical coincidence pushed the author to extend his research from collocation to the most popular finite element method in which closed form analytical formulas of the matrices can be derived. In this formulation, not only the previously found coincidence[22] was repeated, but also it was further found that the famous P-method [2] is also identical with all others.
Although one could simply state that all these polynomials (Lagrange, Bernstein, Legendre, and Chebyshev) have the same basic basis, which is the class { } n x (Taylor series), the possible statement that the Rayleigh quotient is the same has to be mathematically proven.
Another interesting point is the way that the boundary conditions are imposed. As previously mentioned, Lagrange polynomials are associated directly to the nodal values of the primary variable U all-over the domain, whereas the Bernstein (basis) polynomials refer to control points. However, since the extreme control points, P 0 and P n , coincide with the two ends, there is no problem to impose free-free or Dirichlet boundary conditions. Similarly, the p-method is based on the well known linear shape functions associated to the ends of the interval [a,b], and simply its functional set is enriched by the nodeless bubble functions.
In contrast to the two aforementioned polynomials, similar to the Taylor series, Chebyshev polynomials have a more spectral character, as they refer to arbitrary coefficients instead of the pure values of the variable U, i.e. U 0 and U n . It is worthy to mention a common characteristic similarity between Chebyshev polynomials and P-method. In more details, as the degree n of the polynomial increases, the stiffness matrix K n and mass matix M n (both of order (n+1)) can be immediately derived from the previously calculated submatrices K n-1 and M n-1 (both of order n), by simply completing the (n+1)-th row and the (n+1)-th column (it is reminded that in Section 2 the indices vary between 0,…,n).
Despite the coincidence found in one-dimensional problems, preliminary comparisons in two-and three-dimensional problems suggest that there is a significant difference between Lagrangian type finite elements and the p-method, as the first occupy a broader space than the second ones [9,10].

Conclusions
In this work, the entire one-dimensional domain [0,L] was considered as a single macroelement (global approximation), at the ends of which three (all possible) different types of boundary conditions were imposed. Also, for the spatial approximation of the variable U five different types of global higher order polynomials were considered. It was found that all these five basis functions, i.e. (i) classical Lagrange polynomials, (ii) Bernstein polynomials (useful in CAD curve representation), (iii) bubble functions based on Legendre polynomials (p-method), as well as Chebyshev polynomials of (iv) first kind and (v) second kind, share the same func- . Although the coefficients of series expansions and relevant matrices highly depend on the particular polynomial chosen, the calculated eigenvalues were found always the same, a fact that was also rigorously proven initially in case of free-free and then for arbitrary boundary conditions.