Finite Element Preconditioning on Spectral Element Discretizations for Coupled Elliptic Equations

Abstract

The uniform bounds on eigenvalues of 𝖡_ℎ⁻¹₂𝖠_𝑁2 are shown both analytically and numerically by the 𝒫₁ finite element preconditioner 𝖡_ℎ⁻¹₂ for the Legendre spectral element system 𝖠_𝑁2𝐮=𝐟 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:1−Δ𝑢+𝐛⋅∇𝑢+𝑢+𝛿()𝑣=0inΩ,−Δ𝑣−∇⋅𝑣𝐛+𝑣−𝑢=−𝑢inΩ,𝑢=𝑣=0on𝜕Ω, (1.1) 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Δ𝐮𝑢∶=^𝑥𝑥+𝑢_𝑦𝑦𝑣_𝑥𝑥+𝑣_𝑦𝑦,∇𝐮𝑢∶=^𝑥𝑣_𝑥𝑢_𝑦𝑣_𝑦,∇⋅𝐮𝑢∶=^𝑥+𝑢_𝑦𝑣_𝑥+𝑣_𝑦. (2.1) With the usual 𝐿² inner product ⟨⋅,⋅⟩ and its norm ‖⋅‖, for vector functions 𝐮∶=[𝑢,𝑣]^𝑇 and 𝐰∶=[𝑤,𝑧]^𝑇, define𝐮,𝐰∶=⟨𝑢,𝑤⟩+⟨𝑣,𝑧⟩,𝐮²∶=‖𝑢‖²+‖𝑣‖² (2.2) and, for matrix functions 𝑈 and 𝑉, define⟨𝑈,𝑉⟩=⁴_𝑘=1𝑢_𝑘,𝑣_𝑘,‖𝑈‖²∶=⁴_𝑘=1𝑢_𝑘²𝑢,where𝑈=¹𝑢₂𝑢₃𝑢₄𝑣,𝑉=¹𝑣₂𝑣₃𝑣₄. (2.3)

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 1/𝛿, the system (1.1) can be expressed as𝒜𝐮∶=−𝖠Δ𝐮+𝖡𝐛⋅∇𝐮^𝑇+(𝖠+𝖢)𝐮=𝐟inΩ, (2.4) where0𝖠=110𝛿,𝖡=1100−𝛿01,𝖢=𝛿−1𝛿0,𝐟=0−1𝛿𝑢, (2.5) with the zero boundary condition 𝐮=𝟎on𝜕Ω. Let ℬ be another decoupled uniformly elliptic operator such thatℬ𝐮∶=−𝖠Δ𝐮+𝖠𝐮inΩ (2.6) with the zero boundary condition.

2.2. LGL Nodes, Weights, and Function Spaces

Let {𝜂_𝑘}^𝑁_𝑘=0 and {𝜔_𝑘}^𝑁_𝑘=0 be the reference LGL nodes and its corresponding LGL weights in 𝐼=[−1,1], respectively, arranged by −1=∶𝜂₀<𝜂₁<⋯<𝜂_𝑁−1<𝜂_𝑁∶=1. We use {𝑡_𝑗}^𝐸_𝑗=0 as the set of knots in the interval 𝐼 such that −1=∶𝑡₀<𝑡₁<⋯<𝑡_𝐸−1<𝑡_𝐸∶=1. Here 𝐸 denotes the number of subintervals of 𝐼. Denote 𝑁 by the degree of a polynomial on each subinterval 𝐼_𝑗∶=[𝑡_𝑗−1,𝑡_𝑗] and 𝒢∶={𝜉_𝑗,𝑘}^𝐸,𝑁_{𝑗=1,𝑘=0} by the set of 𝑘th-LGL nodes 𝜉_𝑗,𝑘 in each subinterval 𝐼_𝑗 (𝑗=1,…,𝐸) arranged by𝜉_𝑗,0∶=𝑡_𝑗−1<𝜉_𝑗,1<⋯<𝜉_{𝑗,𝑁−1}<𝑡_𝑗=∶𝜉_𝑗,𝑁, (2.7) where𝜉_𝑗,𝑘=ℎ_𝑗2𝜂_𝑘+12𝑡_𝑗−1+𝑡_𝑗,ℎ_𝑗=𝑡_𝑗−𝑡_𝑗−1, (2.8) and the corresponding LGL weights {𝜌_𝑗,𝑘}^𝐸,𝑁_{𝑗=1,𝑘=0} are given by𝜌_𝑗,𝑘ℎ∶=^𝑗2𝜔_𝑁+ℎ_𝑗+12𝜔₀𝜌,𝑘=𝑁,𝑗=1,…,𝐸−1,_{𝑗−1,𝑁,}ℎ𝑘=0,𝑗=2,…,𝐸,_𝑗2𝜔_𝑘,otherwise. (2.9) 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 {𝜙_𝑖(𝑡)}^𝑁_𝑖=0 satisfying𝜙_𝑖𝜂_𝑗=𝛿_𝑖𝑗for𝑖,𝑗=0,1,…,𝑁, (2.10) where 𝛿_𝑖𝑗 denotes the Kronecker delta function.

We define 𝒫^ℎ_𝑁 as the subspace of continuous functions whose basis {𝜙_𝜇}^𝔑+1_𝜇=0 is piecewise continuous Lagrange polynomials of degree 𝑁 on 𝐼_𝑗 with respect to 𝒢 and define 𝒱^ℎ_𝑁 as the space of all piecewise Lagrange linear functions {𝜓_𝜇}^𝔑+1_𝜇=0 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𝜙_𝑝(𝑡)𝜙∶=^𝑝(𝑡),0^𝑇,for𝑝≤𝔑,0,𝜙_𝑝−𝔑(𝑡)^𝑇𝜓,for𝑝>𝑁,^𝑝(𝑡)𝜓∶=^𝑝(𝑡),0^𝑇,for𝑝≤𝔑,0,𝜓_𝑃−𝔑(𝑡)^𝑇,for𝑝>𝑁, (2.11) where 𝑝=1,2,…,2𝔑. 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 {Ξ_𝜇}^𝔑2_𝜇=1∶={(𝜉_𝜇,𝜉_𝜈)}^𝔑,𝔑_{𝜇=1,𝜈=1}, where 𝜇=𝜇+𝔑(𝜈−1) for 𝜇,𝜈=1,2,…,𝔑. 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Φ_𝑝𝜙𝑥,𝑦∶=^𝑝𝑥,𝑦,0^𝑇,for𝑝≤𝔑²,0,𝜙_𝑝−𝔑2𝑥,𝑦^T,for𝑝>𝔑²,Ψ_𝑝𝜓𝑥,𝑦∶=^𝑝𝑥,𝑦,0^𝑇,for𝑝≤𝔑²,0,𝜓_𝑝−𝔑2𝑥,𝑦^𝑇,for𝑝>𝔑², (2.12) where 𝑝=1,…,2𝔑².

2.3. Interpolation Operators

We denote 𝐶(Ω) as the set of continuous functions in Ω∶=𝐼×𝐼. Let ℐ_𝑁2∶𝐶(Ω)→[[𝒫_𝑁]]∶=𝒫_𝑁⊗𝒫_𝑁 be the usual reference interpolation operator such that (ℐ_𝑁2𝑢)(𝜂_𝑖,𝜂_𝑗)=𝑢(𝜂_𝑖,𝜂_𝑗) for 𝑢∈𝐶(Ω) (see, e.g., [17]). The global interpolation operator ℐ^ℎ_𝑁2∶𝐶(Ω)→[[𝒫⁰_𝑁,ℎ]] is given by ℐ^ℎ_𝑁2𝑢(𝑥,𝑦)=∑^𝔑2_𝜇=1𝑢_𝜇𝜙_𝜇(𝑥,𝑦), where 𝑢_𝜇=𝑢(𝜉_𝜇,𝜉_𝜈). Hence, it follows that (ℐ^ℎ_𝑁2𝑢)(𝜉_𝜇,𝜉_𝜈)=𝑢(𝜉_𝜇,𝜉_𝜈) for 𝑢∈𝐶(Ω). With this interpolation operator ℐ^ℎ_𝑁2, let us define the vector interpolation operator 𝐈^ℎ_𝑁2∶𝐶(Ω)×𝐶(Ω)→[[𝐏⁰_𝑁,ℎ]] such that, for 𝐮∶=[𝑢,𝑣]^𝑇∈𝐶(Ω)×𝐶(Ω),𝐈^ℎ_𝑁2𝐮𝜉_𝜇,𝜉_𝜈∶=ℐ^ℎ_𝑁2𝑢𝜉_𝜇,𝜉_𝜈,ℐ^ℎ_𝑁2𝑣𝜉_𝜇,𝜉_𝜈^𝑇=𝐮𝜉_𝜇,𝜉_𝜈. (2.13) Let Φ_̃𝑖(𝑥,𝑦)=𝜙_𝑖𝑗(𝑥,𝑦)=𝜙_𝑖(𝑥)𝜙_𝑗(𝑦) and Ψ_̃𝑖(𝑥,𝑦)=𝜓_𝑖𝑗(𝑥,𝑦)=𝜓_𝑖(𝑥)𝜓_𝑗(𝑦), where 𝑖̃=𝑖+(𝑁+1)𝑗 and 𝑖,𝑗=0,…,𝑁, be the basis of [[𝒫_𝑁]] and [[𝒱_𝑁]], respectively. Let us denote 𝖬_𝑁2 and 𝖬_ℎ2 by the mass matrices such that𝖬_𝑁2𝑖,̃𝑗̃=Φ_̃𝑖,Φ_̃𝑗,𝖬_ℎ2𝑖,̃𝑗̃=Ψ_̃𝑖,Ψ_̃𝑗 (2.14) and denote 𝖲_𝑁2 and 𝖲_ℎ2 by the stiffness matrices such that𝖲_𝑁2𝑖,̃𝑗̃=∇Φ_̃𝑖,∇Φ_̃𝑗,𝖲_ℎ2𝑖,̃𝑗̃=∇Ψ_̃𝑖,∇Ψ_̃𝑗, (2.15) where 𝑖,̃𝑗̃=1,…,(𝑁+1)².

According to Theorems 5.4 and 5.5 in [17], there are two absolute positive constants 𝑐₀ and 𝑐₁ such that for any 𝑈=[𝑢₁,…,𝑢_(𝑁+1)²]^𝑇,𝑐₀𝖬_𝑁2𝑈,𝑈≤𝖬_ℎ2𝑈,𝑈≤𝑐₁𝖬_𝑁2,𝑐𝑈,𝑈₀𝖲_𝑁2𝑈,𝑈≤𝖲_ℎ2𝑈,𝑈≤𝑐₁𝖲_𝑁2𝑈,𝑈, (2.16) and for all 𝑢∈[[𝒱_𝑁]],𝑐₀‖𝑢‖≤ℐ_𝑁2𝑢≤𝑐₁‖𝑢‖,𝑐₀‖𝑢‖₁≤ℐ_𝑁2𝑢₁≤𝑐₁‖𝑢‖₁. (2.17) The extension of (2.17) to the interpolation operator ℐ^ℎ_𝑁2 leads to𝑐₀‖𝑢‖≤ℐ^ℎ_𝑁2𝑢≤𝑐₁‖𝑢‖,𝑐₀‖𝑢‖₁≤ℐ^ℎ_𝑁2𝑢₁≤𝑐₁‖𝑢‖₁ (2.18) 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 𝑐₀𝐮≤𝐈^ℎ_𝑁2𝐮≤𝑐₁𝐮,𝑐₀𝐮₁≤𝐈^ℎ_𝑁2𝐮₁≤𝑐₁𝐮₁. (2.19)

Proof. By the definitions of the interpolation operator 𝐈^ℎ_𝑁2 and the norms, we have 𝐈^ℎ_𝑁2𝐮²=ℐ^ℎ_𝑁2𝑢²+ℐ^ℎ_𝑁2𝑣²,𝐈^ℎ_𝑁2𝐮²₁=ℐ^ℎ_𝑁2𝑢²₁+ℐ^ℎ_𝑁2𝑣²₁, (2.20) 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𝛼𝐮,𝐰=𝖠∇𝐮^𝑇,∇𝐰^𝑇+𝖡𝐛⋅∇𝐮^𝑇,𝐰+(𝖠+𝐶)𝐮,𝐰=𝐟,𝐰,𝛽𝐮,𝐰=𝖠∇𝐮^𝑇,∇𝐰^𝑇+𝖠𝐮,𝐰, (3.1)(3.2) where 𝐮∶=[𝑢(𝑥,𝑦),𝑣(𝑥,𝑦)]^𝑇, 𝐰∶=[𝑤(𝑥,𝑦),𝑧(𝑥,𝑦)]^𝑇 and 𝖠,𝖡,𝖢 are the same matrices in (2.4) and 𝐟=[0,𝑢(𝑥,𝑦)]^𝑇. 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 𝐮₁≤𝛼𝐮,𝐮≤1𝛿𝐮²₁. (3.3)

Proof. Since 𝐛 is a constant vector, [𝑢,𝑣]^𝑇∈𝐻¹₀(Ω)×𝐻¹₀(Ω), and 𝐮 is a real valued vector function, we have 𝛼𝐮,𝐮=𝛽𝐮,𝐮+𝖡𝐛⋅∇𝐮^𝑇,𝐮+𝖢𝐮,𝐮=𝛽𝐮,𝐮+⟨⟩−1𝐛⋅∇𝑢,𝑢𝛿⟨⟩+1𝐛⋅∇𝑣,𝑣𝛿⟩−⟨𝐮(⟨𝑣,𝑢𝑢,𝑣⟩)=𝛽,𝐮=‖𝑢‖²₁+1𝛿‖𝑣‖²₁. (3.4) Hence, (3.3) is proved because 0<𝛿≤1.

Lemma 3.2. If 𝐮=[𝑢,𝑣]^𝑇 is a complex valued function in 𝐻¹₀(Ω)×𝐻¹₀(Ω), then we have the following estimates: Re𝛼𝐮,𝐮=𝛽𝐮,𝐮,𝛼𝐮,𝐮≤11+2|𝐛|+1𝛿𝛽𝐮,𝐮, (3.5) where |𝐛|=𝑏²₁+𝑏²₂.

Proof. Let us decompose 𝑢 and 𝑣 in 𝐮 as 𝑢(𝑥,𝑦)=𝑝(𝑥,𝑦)+𝑖𝑞(𝑥,𝑦) and 𝑣(𝑥,𝑦)=𝑟(𝑥,𝑦)+𝑖𝑠(𝑥,𝑦), where 𝑝,𝑞,𝑟, and 𝑠 are real functions and 𝑖²=−1. Then we have ⟨𝐛⋅∇𝑢,𝑢⟩=𝑖_Ω𝑏₁𝑞_𝑥𝑝−𝑝_𝑥𝑞+𝑏₂𝑝_𝑦𝑞−𝑞_𝑦𝑝−1𝑑Ω,𝛿⟨⟩𝑖𝐛⋅∇𝑢,𝑢=−𝛿_Ω𝑏₁𝑠_𝑥𝑟−𝑟_𝑥𝑠+𝑏₂𝑟_𝑦𝑠−𝑠_𝑦𝑟1𝑑Ω,𝛿⟩−⟨2(⟨𝑣,𝑢𝑢,𝑣⟩)=𝑖𝛿∫_Ω𝑠𝑝−𝑞𝑟𝑑Ω. (3.6)(3.7)(3.8) 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 0<𝛿≤1, we have ⟩−1⟨𝐛⋅∇𝑢,𝑢𝛿⟨⟩+1𝐛⋅∇𝑣,𝑣𝛿⟩−⟨≤|𝐛𝑢‖+1(⟨𝑣,𝑢𝑢,𝑣⟩)|‖∇𝑢‖‖𝛿|𝐛𝑣‖+1|‖∇𝑣‖‖𝛿(2‖𝑣𝑢≤1‖‖‖)2|𝐛|‖∇𝑢‖²+‖𝑢‖²+1𝛿‖∇𝑣‖²+‖𝑣‖²+1𝛿‖𝑢‖²+‖𝑣‖²≤12|𝐛|+1𝛿𝛽𝐮,𝐮. (3.9) Hence (3.5) is proved.

Let 𝜎(𝖠) and ℱ(𝖠) be the spectrum (or set of eigenvalues) and field of values of the square matrix 𝖠, respectively. Let 𝖠_𝑁2 and 𝖡_ℎ2 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 𝑈(≠𝟎)∈ℂ^2𝔑2, one has 𝜎𝖡_ℎ2⁻¹𝖠_𝑁2⊂𝒲, (3.10) where 𝖠𝒲∶=^𝑁2𝑈,𝑈𝖡_ℎ2𝑈,𝑈∣𝑈≠𝟎. (3.11)

Proof. Since 𝖡_ℎ2 is symmetric positive definite, there exists a unique positive definite square root 𝖡_ℎ2^1/2 of 𝖡_ℎ2. So, we have 𝖠_𝑁2𝑈,𝑈𝖡_ℎ2=𝖠𝑈,𝑈^𝑁2𝑈,𝑈𝖡_ℎ2^1/2𝑈,𝖡_ℎ2^1/2𝑈=𝖡_ℎ2^−1/2𝖠_𝑁2𝖡_ℎ2^−1/2𝑉,𝑉⟨⟩𝖡𝑉,𝑉,where𝑉=_ℎ2^1/2𝑈, (3.12) for 𝑈(≠𝟎)∈ℂ^2𝔑2. Now let 𝖡_ℎ2^−1/2 be a short-hand notation for (𝖡_ℎ2^1/2)⁻¹. Therefore, from the relation of spectrum and field of values (see [19] or [11]), it follows that 𝜎𝖡_ℎ2⁻¹𝖠_𝑁2=𝜎𝖡_ℎ2^−1/2𝖠_𝑁2𝖡_ℎ2^−1/2⊂ℱ𝖡_ℎ2^−1/2𝖠_𝑁2𝖡_ℎ2^−1/2=𝒲, (3.13) which completes the proof.

The following theorem shows the uniform bounds of eigenvalues which is independent of both 𝑁 and 𝐸 for our preconditioned system𝖡_ℎ2⁻¹𝖠_𝑁2𝐮^𝑁=𝖡_ℎ2⁻¹𝐟^𝑁. (3.14)

Theorem 3.4. Let {𝜆_𝑝}^2𝔑2_𝑝=1 be the set of eigenvalues of 𝖡_ℎ2⁻¹𝖠_𝑁2, (3.15) then, there are constants 𝑐₀,𝐶₀, and Λ_𝛿 independent of 𝐸 and 𝑁, such that 0<𝑐₀<Re𝜆_𝑝<𝐶₀,𝜆_𝑝≤Λ_𝛿, (3.16) where Λ_𝛿=𝐶(1+‖𝐛‖+1/𝛿).

Proof. Let 𝐮(𝑥,𝑦)∈[[𝐕⁰_𝑁,ℎ]] be represented as 𝐮(𝑥,𝑦)=∑^2𝔑2_𝑝=1𝑢_𝑝Ψ_𝑝(𝑥,𝑦). Then, its piecewise polynomial interpolation can be written as 𝐈^ℎ_𝑁2𝐮𝑥,𝑦=^2𝔑2_𝑝=1𝑢_𝑝Φ_𝑝𝑥,𝑦. (3.17) Let 𝑈=[𝑢₁,…,𝑢_2𝔑2]^𝑇. From the definitions of the bilinear forms, we have 𝖠_𝑁2𝑈,𝑈=𝛼𝐈^ℎ_𝑁2𝐮,𝐈^ℎ_𝑁2𝐮,𝖡_ℎ2𝑈,𝑈=𝛽𝐮,𝐮. (3.18) This implies that 𝑤_𝜆𝖠∶=^𝑁2𝑈,𝑈𝖡_ℎ2=𝛼𝐈𝑈,𝑈^ℎ𝑁2𝐮,𝐈^ℎ_𝑁2𝐮𝛽𝐮,𝐮. (3.19) By Theorem 2.1 and Lemma 3.2, the real part of 𝑤_𝜆 satisfies Re𝑤_𝜆=𝛽𝐈^ℎ_𝑁2𝐮,𝐈^ℎ_𝑁2𝐮𝛽𝐮,𝐮∼1, (3.20) 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 𝑤_𝜆=𝛼𝐈^ℎ_𝑁2𝐮,𝐈^ℎ_𝑁2𝐮𝛽𝐮,𝐮≤(1+(1/2)|𝐛|+(1/𝛿))𝛽𝐈^ℎ_𝑁2𝐮,𝐈^ℎ_𝑁2𝐮𝛽𝐮,𝐮|𝐛|+1≤𝐶1+𝛿. (3.21) Therefore, from Lemma 3.3, the real parts and the absolute values of eigenvalues of 𝖡_ℎ2⁻¹𝖠_𝑁2 satisfy (3.16).

Remark 3.5. Let the one dimensional (1D) bilinear forms for 𝐮,𝐰∈𝐻¹₀(−1,1)×𝐻¹₀(−1,1) be 𝛼₁𝐮,𝐰=𝖠𝐮^′,𝐰^′+𝖡𝑏𝐮^′,𝐰+(𝖠+𝖢)𝐮,𝐰𝛽₁,𝐮,𝐰=𝖠𝐮^′,𝐰^′+𝖠𝐮,𝐰, (3.22)(3.23) 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𝖲_𝑁=𝜙^′_𝜇(𝑡),𝜙^′_𝜈(𝑡),𝖱_𝑁=𝜙^′_𝜇(𝑡),𝜙_𝜈(𝑡),𝖬_𝑁=𝜙_𝜇(𝑡),𝜙_𝜈(𝑡) (3.24) and 1D finite element matrices𝖲_ℎ=𝜓^′_𝜇(𝑡),𝜓^′_𝜈(𝑡),𝖱_ℎ=𝜓^′_𝜇(𝑡),𝜓_𝜈(𝑡),𝖬_ℎ=𝜓_𝜇(𝑡),𝜓_𝜈(𝑡), (3.25) where {𝜙_𝜇}^𝔑_𝜇=1 and {𝜓_𝜈}^𝔑_𝜈=1 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𝖠_𝑁=𝖲_𝑁+𝑏𝖱_𝑁+𝖬_𝑁1𝛿𝖬_𝑁−1𝛿𝖬_𝑁1𝛿𝖲_𝑁−𝑏𝖱_𝑁+𝖬_{𝑁2𝔑×2𝔑},𝖡_ℎ=𝖲_ℎ+𝖬_ℎ𝟎𝟎1𝛿𝖲_ℎ+𝖬_{ℎ2𝔑×2𝔑}. (4.1) Let 𝐛=[𝑏₁,𝑏₂]^𝑇 be a constant vector in (1.1), then the 2D matrices 𝖠_𝑁2 and 𝖡_ℎ2 can be expressed as𝖠_𝑁2=𝖲_𝑁2+𝖱_𝑁2+𝖬_𝑁21𝛿M_𝑁2−1𝛿𝖬_𝑁21𝛿(𝖲_𝑁2−𝖱_𝑁2+𝖬_𝑁2)_{2𝔑2×2𝔑2},𝖡_ℎ2=𝖲_ℎ2+𝖬_ℎ2𝟎𝟎1𝛿𝖲_ℎ2+𝖬_{ℎ22𝔑2×2𝔑2}, (4.2) where𝖲_𝑘2=𝖬_𝑘𝑦⊗𝖲_𝑘(𝑥)+𝖲_𝑘𝑦⊗𝖬_𝑘(𝑥),𝖬_𝑘2=𝖬_𝑘𝑦⊗𝖬_𝑘(𝑥),𝖱)(𝑘=𝑁orℎ_𝑁2=𝑏₁𝖬_𝑁𝑦⊗𝖱_𝑁(𝑥)+𝑏₂𝖱_𝑁𝑦⊗𝖬_𝑁(𝑥). (4.3)

4.2. Numerical Analysis on Eigenvalues

The linear system 𝖠_𝑁2𝐮^𝑁=𝐟^𝑁 and the preconditioned linear system 𝖡_ℎ2⁻¹𝖠_𝑁2𝐮^𝑁=𝖡_ℎ2⁻¹𝐟^𝑁 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 𝖡_ℎ2⁻¹𝖠_𝑁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 𝖠_𝑁2 are presented in Figure 1 for fixed 𝛿=1 (left) and fixed 𝐸=3 (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 𝖠𝑁2 when 𝐛=[1,1]𝑇.

Figures 2 and 4 show the spectra of the resulting preconditioned operator 𝖡_ℎ2⁻¹𝖠_𝑁2 for the polynomials of degrees 𝑁=4,6,8,10 and 𝐸=4 when 𝛿=1 and 𝛿=10⁻⁴, respectively. Also, Figures 3 and 5 show the spectra of 𝖡_ℎ2⁻¹𝖠_𝑁2 for 𝐸=4,6,8,10 and 𝑁=4 for 𝛿=1 and 𝛿=10⁻⁴, respectively. The same axis scales are presented for the same 𝛿 when 𝐛=[1,1]^𝑇. As proven in Theorem 3.4, the eigenvalues of 𝖡_ℎ2⁻¹𝖠_𝑁2 are independent of 𝑁 and 𝐸, but they depend still on the penalty parameter 𝛿.

(a)

(b)

(c)

(d)

Figure 2 

The eigenvalues of 𝖡ℎ2−1𝖠𝑁2 for 𝑁=4,6,8,10 and 𝐸=4 when 𝛿=1,𝐛=[1,1]𝑇.

(a)

(b)

(c)

(d)

Figure 3 

The eigenvalues of 𝖡ℎ2−1𝖠𝑁2 for 𝐸=4,6,8,10 and 𝑁=4 when 𝛿=1,𝐛=[1,1]𝑇.

(a)

(b)

(c)

(d)

Figure 4 

The eigenvalues of 𝖡ℎ2−1𝖠𝑁2 for 𝑁=4,6,8,10 and 𝐸=4 when 𝛿=10−4,𝐛=[1,1]𝑇.

(a)

(b)

(c)

(d)

Figure 5 

The eigenvalues of 𝖡ℎ2−1𝖠𝑁2 for 𝐸=4,6,8,10 and 𝑁=4 when 𝛿=10−4,𝐛=[1,1]𝑇.

By choosing the convection coefficient 𝐛=[10,10]^𝑇, in Figure 6, the distributions of eigenvalues of 𝖠_𝑁2(left) and 𝖡_ℎ2⁻¹𝖠_𝑁2(right) are presented for penalty parameters 𝛿=1,10⁻³ to examine their dependence. In this figure, we see that the distributions of eigenvalues (both real and imaginary part) of 𝖠_𝑁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 𝖡_ℎ2⁻¹𝖠_𝑁2 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 𝖠𝑁2(left) and 𝖡ℎ2−1𝖠𝑁2(right) for 𝛿=1 (top) and 𝛿=10−3 (bottom) when 𝐛=[10,10]𝑇 and 𝑁=12,𝐸=4.

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.