Finite Difference Preconditioners for Legendre Based Spectral Element Methods on Elliptic Boundary Value Problems *

Finite difference type preconditioners for spectral element discretizations based on Legendre-Gauss-Lobatto points are analyzed. The latter is employed for the approximation of uniformly elliptic partial differential problems. In this work, it is shown that the condition number of the resulting preconditioned system is bounded independently of both of the polynomial degrees used in the spectral element method and the element sizes. Several numerical tests verify the h-p independence of the proposed preconditioning.


Introduction
Linear systems engendered by spectral element discretizations of a simple second-order elliptic boundary value problem have large condition numbers depending on the number of elements E and the polynomial degree N employed.Convergence of iterative solvers thus deteriorates as E and N increases.Regardless of these restrictions, the spectral element method is very popular, accurate and used in many engineering problems.However, it is widely known that an efficient preconditioner is necessary in order to improve the convergence of Krylov subspace methods traditionally used to solve the resulting linear system (see [1][2][3][4][5][6][7]). Since the work of [8], it is numerically known that finite difference preconditioning of the spectral element matrix leads to satisfactory results in terms of convergence rates.Multigrid methods are optimal in terms of convergence rate and have linear cost for finite difference problems.The Algebraic multigrid (=: AMG) method can be easily applied to finite difference discretizations of elliptic operators.If it is instead applied directly to high-order discretizations, such as spectral ele-ments, some outstanding issues still need to be addressed.The idea of employing a low-order discretization combined with multigrid as the preconditioner of a high-order problem was investigated in [9] where P 1 finite elements were employed instead of finite differences.Other efforts concerning the computational cost of P 1 finite element based two levels additive Schwarz preconditioners can be found in [10,11].In both approaches an intermediate problem for the Laplace equation was constructed using the high-order Legendre-Gauss-Lobatto (=: LGL) nodes.Furthermore, analytic work was performed in [12] for a second-order uniformly elliptic boundary value problem using LGL nodes and also the analytic research based on Chebyshev-Gauss-Lobatto (=: CGL) nodes was done in [13], in which the various node configurations (LGL and CGL nodes) were employed for the construction of the P 1 finite element preconditioner.In the 1-and 2-dimensional case, the approach was proven optimal and scalable, respectively.
Thus an efficient and optimal algorithm, with linear cost, for solving problems based on spectral element discretizations, which guarantees the convergence of the overall preconditioning strategy, is readily available.The P 1 finite element matrix lowers the complexity of the system to invert since the matrices, representing the Laplacian, are tri, penta or hepta diagonal and are easily solvable using the multigrid method.
Recently Canuto et al. [14,15] have discussed preconditioning with low order finite-elements methods.They show the equivalence between spectral elements and low order finite-elements.The proof is shown for finite-elements and so called numerical integration (NI).In our approach, following motivations found in [4], we propose another preconditioning approach using a simple finite difference scheme based on the LGL nodes in a pointwise approximation perspective.The latter is represented as tensor product matrices, and finite element related properties are employed to prove its optimality.Additionally, it is shown that the proposed preconditioners herein are optimal from the theoretical analysis standpoint and that optimality is confirmed with numerical experiments.
The paper is organized as follows.First, the problem is described in Section 2 and basic definitions are recalled in the same Section 2 also.In the same section, the finite difference preconditioner is introduced.In Section 3, a uniform bound for , the ratio of LGL weights to the distance between the adjacent LGL nodes, is derived which will be used for the 2D case.The optimality results are presented in Sections 4 and 5 with the numerical experiments which support the developed theories.Finally, conclusions are drawn.q 

Description of the Problem
Consider the classical boundary value problem of a second-order elliptic operator where (2. 2) The resulting matrix for the operator A, with variable coefficients, discretized using spectral element methods based on LGL nodes is represented as ˆhN A and 2 2 h N in one-and two-dimensional case, respectively.Twodimensional matrix is the result of tensor products using the matrices obtained in the one-dimensional case.We shall consider the preconditioner h h N , respectively.Moreover, the goal is to show that the preconditioners are optimal in the sense that the condition numbers of the preconditioned systems hN and for each case are independent of mesh sizes

elements and degrees
k of the piecewise polynomials employed in spectral element discretization.

N
For the above goal, we will introduce some notations and definitions from now on.Let be the reference LGL points in   For convenience, we denote j N as the degree of Legendre polynomial on each subinterval 1 : , will be the sum of the degrees of each element up to k: : and, respectively,   , 0 : For simplicity, all LGL points are arranged as 0 1 and the corresponding LGL weights are ordered as Let be the space of all polynomials whose degrees are less than or equal to k and let N  be the subspace of , which consists of piecewise polynomials  whose degrees are less than or equal to j N .Let N  be the space of all piecewise Lagrange linear functions H I with the basis functions   and   . To communicate between the space of piecewise linear functions and the space of piecewise polynomial, we use the interpolation operator such that for .

 
v C I  For the two-dimensional case, the LGL points are re- : , , where  .For the two-dimensional high-order and piecewise linear space, let 2 0 0 0 , , , : , : define the one-dimensional finite difference operator as following: , Finally, the notation for any two real quantities a and b is a shorthand notation corresponding to the existence of two positive constants c and C, which do not depend on the mesh sizes and the degrees of polynomials

 for any two vectors and d
where the superscript T denotes the transpose of a vector or matrix.The standard Sobolev space H 1 and L 2 will be used.And we will use  and 0  as the standard Sobolev H 1 -norm, H 1 -seminorm and the usual L 2 -norm, respectively.

Basic Estimates
We begin by recalling the relations between the distances of LGL nodes and the LGL weights in the reference , : , 1 , , Then the k are uniformly bounded for all q 0,1, , k  .In particular,    are constants independent of the polynomial degree N.
Proof.See Lemma 3.1 in [6].□ Since the goal is to develop preconditioners on spectral elements considering different polynomial degrees on each elements of which sizes are not identical, we need an advanced version or -version of Lemma 3.1.Hence the modifications of (3.1) : and  N independently of both degrees k of polynomials and the mesh sizes .
Since the first and second cases of (3.6) are   it is clear that , are uniformly bounded by Lemma 3.1.
For the case , , it can be easily shown that 2 where are the absolute positive constants that appear in Lemma 3.1.Therefore, which completes the proof.□ Define the two matrices W and H, which are made up of, respectively, the quadrature weights and the distances between the LGL points: where

HU U WU U 
The next result is a clear consequence but we write down for convenience.
Lemma 3.4.For any diagonal matrix S with nonnegative entries whose size is the same as W (see definition in where the equivalence depends on the minimum and maximum entries of S.
The matrix W will be replaced by for twodimensional problems.

W W 
Proposition 3.5.Let be symmetric and positive definite matrices such that for any nonzero vector .Then for any symmetric and positive definite matrix , we have Proof.Since all the matrices are symmetric and positive definite, it is enough to discuss (3.9) in terms of eigenvalues.
Consider the eigenvalue problem .

EU FU
where 1 is a vector consisting of element 1 with length , the vectors and eigenvalues in (3.10) are complete sets of eigenvectors and eigenvalues of the eigenvalue problem From (3.8), since   are positive, bounded and inde- By the same reasoning, it follows Finally, the known results stated in Theorem 3.3 and 3.5 from [12] are recalled here for completeness.
Theorem 3.6.It follows that for all and .

One-Dimensional Case
Expanding and in (4.1) on the space as , the matrix representation of the operator A in (4.1) by the spectral element discretization based on LGL points is now given by where , , D is the differentiation matrix defined as and W is the matrix defined in (3.7).
Since P and Q are the diagonal matrices with nonnegative elements, by Lemma 3.4 we have for any vector U, More precisely, it follows that Besides, we can see that for where .
be the matrix representation of , which is defined in (2.5).For , the easy calculation leads to where .
, , , d U u u u   Note that the constants such as 0 , 1 appeared in this and next sections are generic positive constants, independent of the number of elements i and the degree C E C j N of polynomials applied to spectral discretization.Now to provide a finite difference type preconditioner, consider the bilinear form with a positive constant . This induces the matrix : , , , , , Lemma 4.1.For any vector , it follows that Proof.Note that can be expressed as , so that its piecewise polynomial The last inequality is due to Poincare's inequality.
Therefore, (4.7), (4.6) and Theorem 3.6 complete the proof.□ Remark 1. From the above lemma, one may notice that the upper bound for the condition numbers of the matrix does not depend on the mesh size h j of an element and the degree N i of a polynomial.Further, it does not depend on 1 ,0 ĥ h N L A   .Hence, even if (4.1) has coeffi- cient functions and , it is enough to take the preconditioner with Next, we consider another bilinear form , , , , , , Now the goal is to analyze the validity of the matrix operator for the preconditioner to ˆh L ˆhN A .Theorem 4.2.For any vector , it follows that where the equivalence depends only on , , , , ˆu p ˆl p ˆu q ˆl q  and  .
Proof.For such that , its piecewise polynomial interpolation is On the other hand, applying Corollary 3.3, (4.5) and Theorem 3.6 we have Hence, using Theorem 3.6 again, we have To guarantee the positivity of the lower bound, if , we will take ˆ0 l q  0   .Applying Poincare's inequality with (4.9) and (4.6), we have for some positive constant C. Therefore, using (4.3)-(4.5)and Theorem 3.6, we get which leads to the conclusion.The latter, combined with the min-max theorem, yields the next result.  of are all positive real and bounded above and below.The bounds are independent of the mesh sizes h j and the degrees of the polynomials N k .That is, there are absolute constants and such that 2 reveals that the condition numbers of are bounded above by . Hence one    may notice that the condition numbers do not depend on the degrees of polynomials and the mesh sizes k N j h .According to Remark 2, we may summarize the behavior of the condition numbers of in the following corollary.We will investigate the efficiency of the preconditioners for several problems with constant and variable coefficients.Moreover, for a variable coefficients problem, we will compare the developed preconditioners with the 1 finite element preconditioner (see [12]) in terms of iteration numbers using the preconditioned conjugate gradient method (=: PCG).P Example 1.Consider the operator with constant coefficients: with homogeneous Dirichlet boundary conditions, where and are constants.Note that the spectral element discretized matrix corresponding to (4.2) becomes Since the condition numbers are dependent on the ratios q p and  by the definition, we take 1   and focus mainly on the role of  , where ˆH ĥ q p -value equal or less than 1 (see Table 1).
Second, to show that the preconditioning work is not independent of the polynomial degrees N and the numbers of element E, applied in spectral element method, it is tested for the cases with  with several  -values for each q p lue (see Figure 3).When -va  is t n to be equal to ake q p -value or   , where ˆ0 p  the spectra ment discretization yield the matrix , where P and are the Q ma es of 2   q/p=10 3 q/p=10 2 q/p=10 q/p=1 q/p=10 −1 q/p=10 −2 q/p=10 −3 element discretization, respectively.Figure 5

Figure 5. Condition numbers of where
we compute the iteration number using PCG preconditioned by ˆĥ L B  .Also we compare the developed preconditioners with the 1 finite element preconditioner ˆh L h P F (see [10,12] for example) in terms of iteration numbers using PCG.
The computational results are shown in Table 2 with various numbers of elements E and polynomial degrees N.While the number of iterations increases as N or E becomes large in case of CG iterations, the PCG preconditioned by gives relatively small.In particular, it can be predicted that the preconditioning effects become stronger as N or E are made larger.Furthermore we note that the results are comparable to the ones for finite element preconditioner.ˆh L

Two-Dimensional Case
In this section the tensor notation is employed.It has the form and the superscripts denote the spatial dimension on which it acts.Their order will always be the same and thus we omit them.For the tensor product representation, we refer to [16].
Consider the elliptic operator corresponding to 2D case with zero boundary conditions: where q q x y q     , .As in the 1D case, by expanding the variable coefficients   , p x y and   , q x y in terms of the basis , two coefficients matrices are obtained: Then the matrix from the spectral element discretization, using one-dimensional matrices and W, becomes where M is the 1D mass matrix, which becomes the identity for the Lagrange basis From the similar argument in 1 (4.4), we have

 for any nonzero vector
In order to describe the preconditioner for the defined in (2.6), we c si V .
2D on-case, with the operator h L 2 der a bilinear form where As in the 1D case, we compare the number of PCG iterations for the preconditioners 2 ˆh L and 2 ˆh F which provided [12].In Table 3, we can see the suggested finite difference preconditioner is more is in that ← β = q 1 /4p

Conclusion
We have proposed finite differences prec rs and for spectral element discretizations, and c on LGL nodes for uniformly e ptic pr n-d ension, , resp ively. he tw o tioners optim he c ponding spectral ele ments problems was emonstrated th ugh the theor cal proofs and the c putational res s.Th rden the efficiency is now on t ultigrid algorithm for solving finite-differe es pr ms and not G igh-order elements.

REFERENCES
and 2D case, respectively.Then the problem is to demonstrate the validity of the preconditioners for translated LGL points and corresponding weights from and   are N j h of each element.Proof.Since the LGL nodes in each subinterval are defined by an affine map transformation from the ones in reference interval   1,1  , they depend on the mesh sizes of each element thus they are distinguished as follows Copyright © 2013 SciRes.AM S. KIM ET AL.

Corollary 4 . 3 .
The eigenvalues   1 d   investigated the condition numbers of the matrices by preconditioned by ˆĥ

Figure 1 .
Figure 1.Condition numbers of where
Figure 4. Condition numbers of , where

Table 3 . Iteration numbers
1 Figure 6.Condition numbers of