Journal of Applied MathematicsVolume 2012 (2012), Article ID 245051, 16 pagesdoi:10.1155/2012/245051

Research Article

Finite Element Preconditioning on Spectral Element Discretizations for Coupled Elliptic Equations

JongKyum Kwon, Soorok Ryu, Philsu Kim, and Sang Dong Kim

Department of Mathematics, Kyungpook National University, Daegu 702-701, Republic of Korea

Received 8 August 2011; Accepted 2 November 2011

Academic Editor: Massimiliano Ferronato

Copyright © 2012 JongKyum Kwon et al. JongKyum Kwon et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The uniform bounds on eigenvalues of are shown both analytically and numerically by the ₁ finite element preconditioner for the Legendre spectral element system which is arisen from a coupled elliptic system occurred by an optimal control problem. The finite element preconditioner is corresponding to a leading part of the coupled elliptic system.

1. 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 optimal control problem subject to an elliptic equation (see [1]). Applying the 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 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.

2. Preliminaries

2.1. 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 , where denotes the transpose, are defined by With the usual inner product and its norm , for vector functions and , define and, for matrix functions and , define

We use the standard Sobolev spaces like and on a given domain with the usual Sobolev seminorm and norm . The main content of this paper is to provide an efficient low-order preconditioner for the system (1.1).

Multiplying the second equation by , the system (1.1) can be expressed as where with the zero boundary condition . Let be another decoupled uniformly elliptic operator such that with the zero boundary condition.

2.2. LGL Nodes, Weights, and Function Spaces

Let and be the reference LGL nodes and its corresponding LGL weights in , respectively, arranged by . We use as the set of knots in the interval such that . Here denotes the number of subintervals of . Denote by the degree of a polynomial on each subinterval and by the set of -LGL nodes in each subinterval () arranged by where and the corresponding LGL weights are given by Let be the space of all polynomials defined on whose degrees are less than or equal to . The Lagrange basis for on is given by satisfying where denotes the Kronecker delta function.

We define as the subspace of continuous functions whose basis is piecewise continuous Lagrange polynomials of degree on with respect to and define as the space of all piecewise Lagrange linear functions with respect to . Note that these basis functions have a proper support. See [9, 18] for detail. Let and . With the notation , define and as subspaces of and , respectively. The basis functions for and are given respectively by where . For two dimensional (2D) case, let and be tensor product function spaces of one-dimensional function spaces and let and be subspaces of and , respectively. Now let us order the interior LGL points in by horizontal lines as , where for . Accordingly, the basis functions for and are also arranged as the same way. Then, with the notations and , the basis functions of and are given, respectively, by where .

2.3. Interpolation Operators

We denote as the set of continuous functions in . Let be the usual reference interpolation operator such that for (see, e.g., [17]). The global interpolation operator is given by , where . Hence, it follows that for . With this interpolation operator , let us define the vector interpolation operator such that, for , Let and , where and , be the basis of and , respectively. Let us denote and by the mass matrices such that and denote and by the stiffness matrices such that where .

According to Theorems 5.4 and 5.5 in [17], there are two absolute positive constants and such that for any , and for all , The extension of (2.17) to the interpolation operator leads to for all , where the constants and are positive constants independent of and (see [18]).

Theorem 2.1. For all , there are positive constants and independent of and such that

Proof. By the definitions of the interpolation operator and the norms, we have which completes the proof because of (2.18).

3. Analysis on Finite Element Preconditioner

The bilinear forms corresponding to (2.4) and (2.6) are given by where , and are the same matrices in (2.4) and . 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 for the variational problem (3.1).

Proposition 3.1. For a real valued vector function , we have

Proof. Since is a constant vector, , and is a real valued vector function, we have Hence, (3.3) is proved because .

Lemma 3.2. If is a complex valued function in , then we have the following estimates: where .

Proof. Let us decompose and in as and , where , and are real functions and . Then we have Hence, one may see that the real part of is 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 , we have Hence (3.5) is proved.

Let and be the spectrum (or set of eigenvalues) and field of values of the square matrix , respectively. Let and be the two dimensional stiffness matrices on the spaces and induced by (3.1) and (3.2), respectively. Then, we have the following.

Lemma 3.3. For , one has where

Proof. Since is symmetric positive definite, there exists a unique positive definite square root of . So, we have for . Now let be a short-hand notation for . Therefore, from the relation of spectrum and field of values (see [19] or [11]), it follows that which completes the proof.

The following theorem shows the uniform bounds of eigenvalues which is independent of both and for our preconditioned system

Theorem 3.4. Let be the set of eigenvalues of then, there are constants , and independent of and , such that where .

Proof. Let be represented as . Then, its piecewise polynomial interpolation can be written as Let . From the definitions of the bilinear forms, we have This implies that By Theorem 2.1 and Lemma 3.2, the real part of satisfies where the notation means the equivalence of two quantities and which does not depend on and . Again, from Theorem 2.1 and Lemma 3.2, the absolute value of satisfies Therefore, from Lemma 3.3, the real parts and the absolute values of eigenvalues of satisfy (3.16).

Remark 3.5. Let the one dimensional (1D) bilinear forms for be and denote and by the 1D stiffness matrices on the spaces and 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 and are the Lagrange basis of the spaces and , respectively. For actual computations for , , and , the inner product on the space 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 , , and and the exact matrices , , and are equivalent, respectively, because of the equivalence of LGL quadrature on the polynomial space we used. For example, the mass matrix can be computed using the LGL weights only.

4. Numerical Tests of Preconditioning

4.1. 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 and corresponding to (3.22) and (3.23) using the matrices in (3.24). One may have Let be a constant vector in (1.1), then the 2D matrices and can be expressed as where

4.2. Numerical Analysis on Eigenvalues

The linear system and the preconditioned linear system 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 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 are presented in Figure 1 for fixed (left) and fixed (right) as increasing the degrees of polynomials. It shows that such condition numbers depend on , , and . In particular, the smaller is, the larger condition number it yields.

(a)

(b)

Figure 1: The condition numbers of the matrix when .

Figures 2 and 4 show the spectra of the resulting preconditioned operator for the polynomials of degrees and when and , respectively. Also, Figures 3 and 5 show the spectra of for and for and , respectively. The same axis scales are presented for the same when . As proven in Theorem 3.4, the eigenvalues of are independent of and , but they depend still on the penalty parameter .

(a)

(b)

(c)

(d)

Figure 2: The eigenvalues of for and when .

(a)

(b)

(c)

(d)

Figure 3: The eigenvalues of for and when .

(a)

(b)

(c)

(d)

Figure 4: The eigenvalues of for and when .

(a)

(b)

(c)

(d)

Figure 5: The eigenvalues of for and when .

By choosing the convection coefficient , in Figure 6, the distributions of eigenvalues of (left) and (right) are presented for penalty parameters to examine their dependence. In this figure, we see that the distributions of eigenvalues (both real and imaginary part) of are strongly dependent on . The real parts of such eigenvalues are increased in proportion to . On the other hand, as predicted by Theorem 3.4, the real parts of the eigenvalues of are positive and uniformly bounded away from 0. Moreover, the real parts are not dependent on and (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 . These phenomena support the theory proved in Theorem 3.4.

(a)

(b)

(c)

(d)

Figure 6: The eigenvalues of (left) and (right) for (top) and (bottom) when and .

5. 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]).

Acknowledgments

The authors would like to thank Professor. Gunzburg for his kind advice to improve this paper. This work was supported by KRF under contract no. C00094.

References

M. D. Gunzburger, Perspectives in Flow Control and Optimization, vol. 5 of Advances in Design and Control, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2003.

R. Li, W. Liu, H. Ma, and T. Tang, “Adaptive finite element approximation for distributed elliptic optimal control problems,” SIAM Journal on Control and Optimization, vol. 41, no. 5, pp. 1321–1349, 2002.

P. Bochev and M. D. Gunzburger, “Least-squares finite element methods for optimality systems arising in optimization and control problems,” SIAM Journal on Numerical Analysis, vol. 43, no. 6, pp. 2517–2543, 2006.

A. Borzì, K. Kunisch, and D. Y. Kwak, “Accuracy and convergence properties of the finite difference multigrid solution of an optimal control optimality system,” SIAM Journal on Control and Optimization, vol. 41, no. 5, pp. 1477–1497, 2003.

S. Ryu, H.-C. Lee, and S. D. Kim, “First-order system least-squares methods for an optimal control problem by the Stokes flow,” SIAM Journal on Numerical Analysis, vol. 47, no. 2, pp. 1524–1545, 2009.

Y. Chen, N. Yi, and W. Liu, “A Legendre-Galerkin spectral method for optimal control problems governed by elliptic equations,” SIAM Journal on Numerical Analysis, vol. 46, no. 5, pp. 2254–2275, 2008.

R. Ghanem and H. Sissaoui, “A posteriori error estimate by a spectral method of an elliptic optimal control problem,” Journal of Computational Mathematics and Optimization, vol. 2, no. 2, pp. 111–135, 2006.

A. T. Patera, “A spectral element method for fluid dynamics: Laminar flow in a channel expansion,” Journal of Computational Physics, vol. 54, pp. 468–488, 1984.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods : Fundamentals in Single Domains, Springer, Berlin, Germany, 2006.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods : Evolution to Complex Geometries and Applications to Fluid Dynamics, Springer, Berlin, Germany, 2006.

Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Ruplishing Company, Boston, Mass, USA, 1966.

L. N. Trefethen and D. B. Bau, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 1997.

P. F. Fischer, “An overlapping Schwarz method for spectral element solution of the incompressible Navier-Stokes equations,” Journal of Computational Physics, vol. 133, no. 1, pp. 84–101, 1997.

S. Kim and S. D. Kim, “Preconditioning on high-order element methods using Chebyshev-Gauss-Lobatto nodes,” Applied Numerical Mathematics, vol. 59, no. 2, pp. 316–333, 2009.

S. D. Kim and S. V. Parter, “Preconditioning Chebyshev spectral collocation method for elliptic partial differential equations,” SIAM Journal on Numerical Analysis, vol. 33, no. 6, pp. 2375–2400, 1996.

S. D. Kim and S. V. Parter, “Preconditioning Chebyshev spectral collocation by finite-difference operators,” SIAM Journal on Numerical Analysis, vol. 34, no. 3, pp. 939–958, 1997.

S. V. Parter and E. E. Rothman, “Preconditioning Legendre spectral collocation approximations to elliptic problems,” SIAM Journal on Numerical Analysis, vol. 32, no. 2, pp. 333–385, 1995.

S. D. Kim, “Piecewise bilinear preconditioning of high-order finite element methods,” Electronic Transactions on Numerical Analysis, vol. 26, pp. 228–242, 2007.

R. A. Horn and C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, Cambridge, UK, 1991.

M. O. Deville, P. F. Fischer, and E. H. Mund, High-Order Methods for Incompressible Fluid Flow, vol. 9 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 2002.