Finite Element Preconditioning on Spectral Element Discretizations for Coupled Elliptic Equations

The uniform bounds on eigenvalues of (cid:2) B − 1 h 2 (cid:2) A N 2 are shown both analytically and numerically by the P 1 ﬁnite element preconditioner (cid:2) B − 1 h 2 for the Legendre spectral element system (cid:2) A N 2 u (cid:2) f which is arisen from a coupled elliptic system occurred by an optimal control problem. The ﬁnite element preconditioner is corresponding to a leading part of the coupled elliptic system.


Introduction
Optimal control problems constrained by partial differential equations can be reduced to a system of coupled partial differential equations by Lagrange multiplier method 1 .In particular, the needs for accurate and efficient numerical methods for these problems have been important subjects.Many works are reported for solving coupled partial differential equations by finite element/difference methods; or finite element least-squares methods 2-5 , etc. .But, there are a few literature for examples, 6, 7 on coupled partial differential equations using the spectral element methods SEM despite of its popularity and accuracy see, e.g., 8 .
One of the goals in this paper is to investigate a finite element preconditioner for the SEM discretizations.The induced nonsymmetric linear systems by the SEM discretizations from such coupled elliptic partial differential equations have the condition numbers which are getting larger incredibly not only as the number of elements and degrees of polynomials increases but also as the penalty parameter δ decreases see 5 and Section 4 .Hence, an efficient preconditioner is necessary to improve the convergence of a numerical method whose number of iterations depends on the distributions of eigenvalues see 9-12 .Particularly, the lower-order finite element/difference preconditioning methods for spectral collocation/element methods have been reported 9, 10, 13-17 , etc. .The target-coupled elliptic type equations are as follows: which is the result of Lagrange multiplier rule applied to a L 2 optimal control problem subject to an elliptic equation see 1 .Applying the P 1 finite element preconditioner to our coupled elliptic system discretized by SEM using LGL Legendre-Gauss-Lobatto nodes, we show that the preconditioned linear systems have uniformly bounded eigenvalues with respect to elements and degrees.
The field of values arguments will be used instead of analyzing eigenvalues directly because the matrix representation of the target operator, even with zero convection term, is not symmetric.We will show that the real parts of eigenvalues are positive, uniformly bounded away from zero, and the absolute values of eigenvalues are uniformly bounded whose bounds are only dependent on the penalty parameter δ in 1.1 and the constant vector b in 1.1 .Because of this result, one may apply a lower-order finite element preconditioner to a real optimal control problem subject to Stokes equations which requires an elliptic type solver.
This paper is organized as follows: in Section 2, we introduce some preliminaries and notations.The norm equivalences of interpolation operators are reviewed to show the norm equivalence of an interpolation operator using vector basis.The preconditioning results are presented theoretically and numerically in Section 3 and Section 4, respectively.Finally, we add the concluding remarks in Section 5.

Coupled Elliptic System
Because we are going to deal with a coupled elliptic system, the vector Laplacian, gradient and divergence operators for a vector function u u, v T , where T denotes the transpose, are defined by With the usual L 2 inner product •, • and its norm • , for vector functions u : u, v T and w : w, z T , define and, for matrix functions U and V , define We use the standard Sobolev spaces like H 1 Ω and H 1 0 Ω on a given domain Ω with the usual Sobolev seminorm | • | 1 and norm • 1 .The main content of this paper is to provide an efficient low-order preconditioner for the system 1.1 .
Multiplying the second equation by 1/δ, the system 1.1 can be expressed as where with the zero boundary condition u 0 on ∂Ω.Let B be another decoupled uniformly elliptic operator such that with the zero boundary condition.

LGL Nodes, Weights, and Function Spaces
Let {η k } N k 0 and {ω k } N k 0 be the reference LGL nodes and its corresponding LGL weights in I −1, 1 , respectively, arranged by −1 : We use {t j } E j 0 as the set of knots in the interval I such that −1 : Here E denotes the number of subintervals of I. Denote N by the degree of a polynomial on each subinterval I j : t j−1 , t j and G : {ξ j,k } E,N j 1,k 0 by the set of kth-LGL nodes ξ j,k in each subinterval I j j 1, . . ., E arranged by where and the corresponding LGL weights {ρ j,k } E,N j 1,k 0 are given by ρ j,k : 2.9 Let P k be the space of all polynomials defined on I whose degrees are less than or equal to k.
The Lagrange basis for P N on I is given by { φ i t } N i 0 satisfying where δ ij denotes the Kronecker delta function.
We define P h N as the subspace of continuous functions whose basis {φ μ } N 1 μ 0 is piecewise continuous Lagrange polynomials of degree N on I j with respect to G and define V h N as the space of all piecewise Lagrange linear functions {ψ μ } N 1 μ 0 with respect to G. Note that these basis functions have a proper support.See 9, 18 for detail.Let P 0 N,h : h be tensor product function spaces of one-dimensional function spaces and let P 0 N,h , respectively.Now let us order the interior LGL points in Ω by horizontal lines as , where μ μ N ν−1 for μ, ν 1, 2, . . ., N. Accordingly, the basis functions for P 0 N,h and V 0 N,h are also arranged as the same way.Then, with the notations φ μ x, y : φ μ x φ ν y and ψ μ x, y : ψ μ x ψ ν y , the basis functions of P 0 N,h and V 0 N,h are given, respectively, by

Interpolation Operators
We denote C Ω as the set of continuous functions in Ω : I × I. Let I N 2 : C Ω → P N : P N ⊗ P N be the usual reference interpolation operator such that I N 2 u η i , η j u η i , η j for u ∈ C Ω see, e.g., 17 .The global interpolation operator μ 1 u μ φ μ x, y , where u μ u ξ μ , ξ ν .Hence, it follows that Let Φ i x, y φ ij x, y φ i x φ j y and Ψ i x, y ψ ij x, y ψ i x ψ j y , where i i N 1 j and i, j 0, . . ., N, be the basis of P N and V N , respectively.Let us denote M N 2 and M h 2 by the mass matrices such that and denote S N 2 and S h 2 by the stiffness matrices such that where i, j 1, . . ., N 1 2 .According to Theorems 5.4 and 5.5 in 17 , there are two absolute positive constants c 0 and c 1 such that for any U u 1 , . . ., u N 1 2 T , and for all u ∈ V N , The extension of 2.17 to the interpolation operator I h N 2 leads to for all u ∈ V 0 N,h , where the constants c 0 and c 1 are positive constants independent of E and N see 18 .
Theorem 2.1.For all u ∈ V 0 N,h , there are positive constants c 0 and c 1 independent of E and N such that

2.19
Proof.By the definitions of the interpolation operator I h N 2 and the norms, we have which completes the proof because of 2.18 .

Analysis on P 1 Finite Element Preconditioner
The bilinear forms corresponding to 2.4 and 2.6 are given by where u : u x, y , v x, y T , w : w x, y , z x, y T and A, B, C are the same matrices in 2.4 and f 0, u x, y T .Note that the bilinear form β •, • in 3.2 is symmetric but the bilinear form α •, • in 3.1 is not symmetric.The following norm equivalence guarantees the existence and uniqueness of the solution in H 1 0 Ω × H 1 0 Ω for the variational problem 3.1 .
Proposition 3.1.For a real valued vector function u x, y u, v T ∈ H 1 0 Ω × H 1 0 Ω , we have , and u is a real valued vector function, we have
Lemma 3.2.If u u, v T is a complex valued function in H 1 0 Ω × H 1 0 Ω , then we have the following estimates: Proof.Let us decompose u and v in u as u x, y p x, y iq x, y and v x, y r x, y is x, y , where p, q, r, and s are real functions and i 2 −1.Then we have Hence, one may see that the real part of α u, u is β u, u and the pure imaginary part is the sum of 3.6 , 3.7 , and 3.8 .By Cauchy-Schwarz inequality, -inequality, and the range of 3.9 Hence 3.5 is proved.
Let σ A and F A be the spectrum or set of eigenvalues and field of values of the square matrix A, respectively.Let A N 2 and B h 2 be the two dimensional stiffness matrices on the spaces P 0 N,h and V 0 N,h induced by 3.1 and 3.2 , respectively.Then, we have the following.

Lemma 3.3. For
where Proof.Since B h 2 is symmetric positive definite, there exists a unique positive definite square root B 1/2 h 2 of B h 2 .So, we have , where V B which completes the proof.
The following theorem shows the uniform bounds of eigenvalues which is independent of both N and E for our preconditioned system

3.14
Theorem 3.4.Let {λ p } 2N 2 p 1 be the set of eigenvalues of Journal of Applied Mathematics 9 then, there are constants c 0 , C 0 , and Λ δ independent of E and N, such that where be represented as u x, y 2N 2 p 1 u p Ψ p x, y .Then, its piecewise polynomial interpolation can be written as 3.17 Let U u 1 , . . ., u 2N 2 T .From the definitions of the bilinear forms, we have This implies that w λ : By Theorem 2.1 and Lemma 3.2, the real part of w λ satisfies Re w λ where the notation a ∼ b means the equivalence of two quantities a and b which does not depend on E and N. Again, from Theorem 2.1 and Lemma 3.2, the absolute value of w λ satisfies

3.21
Therefore, from Lemma 3.3, the real parts and the absolute values of eigenvalues of B −1 h 2 A N 2 satisfy 3.16 .

Remark 3.5. Let the one dimensional 1D bilinear forms for u, w ∈ H
and denote A N and B h by the 1D stiffness matrices on the spaces P 0 N,h and V 0 N,h corresponding to the bilinear forms 3.22 and 3.23 , respectively.Then one can easily get the same results as Theorem 3.4 for 1D case.Since the proof is similar to 2D case, we omit the statements.Now, for actual numerical computations, we need 2D stiffness and mass matrices expressed as the tensor products of 1D matrices see 20 for details .For this, let us denote 1D spectral element matrices as and 1D finite element matrices where {φ μ } N μ 1 and {ψ ν } N ν 1 are the Lagrange basis of the spaces P 0 N,h and V 0 N,h , respectively.For actual computations for S N , R N , and M N , the inner product •, • on the space P 0 N,h will be computed using LGL quadrature rule at LGL nodes.Without any confusion, such approximate matrices can be denoted by same notations in the next section.We note that the approximate matrices S N , R N , and M N and the exact matrices S N , R N , and M N are equivalent, respectively, because of the equivalence of LGL quadrature on the polynomial space we used.For example, the mass matrix M N can be computed using the LGL weights {ρ j,k } only.

Matrix Representation
In this section we discuss effects of the proposed finite element preconditioning for the spectral element discretizations to the coupled elliptic system 1.1 .For this purpose, first we set up one dimensional matrices A N and B h corresponding to 3.22 and 3.23 using the matrices in 3.24 .One may have T be a constant vector in 1.1 , then the 2D matrices A N 2 and B h 2 can be expressed as where

Numerical Analysis on Eigenvalues
The linear system A N 2 u N f N and the preconditioned linear system B −1 N will be compared in the sense of the distribution of eigenvalues.As proved in Section 3, it is shown that the behaviors of spectra of B −1 h 2 A N 2 are independent of the number of elements and degrees of polynomials.
One may also see the condition numbers of these discretized systems by varying the penalty parameter δ.The condition numbers of A N 2 are presented in Figure 1 for fixed δ 1 left and fixed E 3 right as increasing the degrees N of polynomials.It shows that such condition numbers depend on N, E, and δ.In particular, the smaller δ is, the larger condition number it yields.Figures 2 and 4 show the spectra of the resulting preconditioned operator B −1 h 2 A N 2 for the polynomials of degrees N 4, 6, 8, 10 and E 4 when δ 1 and δ 10 −4 , respectively.Also, Figures 3 and 5 show the spectra of B −1 h 2 A N 2 for E 4, 6, 8, 10 and N 4 for δ 1 and δ 10 −4 , respectively.The same axis scales are presented for the same δ when b 1, 1 T .As proven in Theorem 3.4, the eigenvalues of B −1 h 2 A N 2 are independent of N and E, but they depend still on the penalty parameter δ.
By choosing the convection coefficient b 10, 10 T , in Figure 6, the distributions of eigenvalues of A N 2 left and B −1 h 2 A N 2 right are presented for penalty parameters δ 1, 10 −3 to examine their dependence.In this figure, we see that the distributions of eigenvalues both real and imaginary part of A N 2 are strongly dependent on δ.The real parts of such eigenvalues are increased in proportion to 1/δ.On the other hand, as predicted by Theorem 3.4, the real parts of the eigenvalues of B −1 h 2 A N 2 are positive and uniformly bounded away from 0. Moreover, the real parts are not dependent on δ and b see Figures 2-6 , and their moduli are uniformly bounded.The numerical results show that the imaginary parts of the eigenvalues are bounded by some constants which are dependent on δ and b.These phenomena support the theory proved in Theorem 3.4.

Concluding Remarks
An optimal control problem subject to an elliptic partial differential equation yields coupled elliptic differential equations 1.1 .Any kind of discretizations leads to a nonsymmetric linear system which may require Krylov subspace methods to solve the system.In this paper, the spectral element discretization is chosen because it is very accurate and popular, but the resulting linear systems have large condition numbers.This situation now becomes one of disadvantages if one aims at a fast and efficient numerical simulation for an optimal control problem subject to even a simple elliptic differential equation.To overcome such a disadvantage, the lower-order finite element preconditioner is proposed so that the preconditioned linear system has uniformly bounded condition numbers independent of the degrees of polynomials and the mesh sizes.One may also take various degrees of polynomials on subintervals with different mesh sizes.In this case, similar results can be obtained without any difficulties.This kind of finite element preconditioner may be used for an optimal control problem subject to Stokes flow see, e.g., 13 .