# Finite Element Preconditioning on Spectral Element Discretizations for Coupled Elliptic Equations

Authors: JongKyum Kwon; Soorok Ryu; Philsu Kim; Sang Dong Kim
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

The uniform bounds on eigenvalues ofB^h2−1A^N2 are shown both analytically and numerically by
the 𝒫1 finite element preconditioner B^h2−1 for the Legendre spectral element system A^N2u¯=f¯ 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.


---

## Body

## 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.1)-Δu+b⋅∇u+u+1δv=0inΩ,-Δv-∇⋅(vb)+v-u=-ûinΩ,u=v=0on⁡∂Ω,
which is the result of Lagrange multiplier rule applied to a L2 optimal control problem subject to an elliptic equation (see [1]). Applying the 𝒫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 Section2, 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 functionu̲=[u,v]T, where T denotes the transpose, are defined by(2.1)Δu̲∶=[uxx+uyyvxx+vyy],∇u̲∶=[uxvxuyvy],∇⋅u̲∶=[ux+uyvx+vy].
With the usual L2 inner product 〈·,·〉 and its norm ∥·∥, for vector functions u̲∶=[u,v]T and w̲∶=[w,z]T, define(2.2)〈u̲,w̲〉∶=〈u,w〉+〈v,z〉,‖u̲‖2∶=‖u‖2+‖v‖2
and, for matrix functions U and V, define(2.3)〈U,V〉=∑k=14〈uk,vk〉,‖U‖2∶=∑k=14‖uk‖2,whereU=[u1u2u3u4],V=[v1v2v3v4].We use the standard Sobolev spaces likeH1(Ω) and H01(Ω) 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 by1/δ, the system (1.1) can be expressed as(2.4)Au̲∶=-AΔu̲+B(b⋅∇u̲)T+(A+C)u̲=f̲inΩ,
where(2.5)A=[1001δ],B=[100-1δ],C=[01δ-1δ0],f̲=[0-1δû],
with the zero boundary condition u̲=0̲on∂Ω. Let ℬ be another decoupled uniformly elliptic operator such that(2.6)Bu̲∶=-AΔu̲+Au̲inΩ
with the zero boundary condition.


### 2.2. LGL Nodes, Weights, and Function Spaces

Let{ηk}k=0N and {ωk}k=0N be the reference LGL nodes and its corresponding LGL weights in I=[-1,1], respectively, arranged by -1=:η0<η1<⋯<ηN-1<ηN∶=1. We use {tj}j=0E as the set of knots in the interval I such that -1=:t0<t1<⋯<tE-1<tE∶=1. Here E denotes the number of subintervals of I. Denote N by the degree of a polynomial on each subinterval Ij∶=[tj-1,tj] and 𝒢∶={ξj,k}j=1,k=0E,N by the set of kth-LGL nodes ξj,k in each subinterval Ij (j=1,…,E) arranged by(2.7)ξj,0∶=tj-1<ξj,1<⋯<ξj,N-1<tj=:ξj,N,
where(2.8)ξj,k=hj2ηk+12(tj-1+tj),hj=tj-tj-1,
and the corresponding LGL weights {ρj,k}j=1,k=0E,N are given by(2.9)ρj,k∶={hj2ωN+hj+12ω0,k=N,j=1,…,E-1,ρj-1,N,k=0,j=2,…,E,hj2ωk,otherwise.
Let 𝒫k be the space of all polynomials defined on I whose degrees are less than or equal to k. The Lagrange basis for 𝒫N on I is given by {ϕ̂i(t)}i=0N satisfying(2.10)ϕ̂i(ηj)=δijfori,j=0,1,…,N,
where δij denotes the Kronecker delta function.We define𝒫Nh as the subspace of continuous functions whose basis {ϕμ}μ=0𝔑+1 is piecewise continuous Lagrange polynomials of degree N on Ij with respect to 𝒢 and define 𝒱Nh as the space of all piecewise Lagrange linear functions {ψμ}μ=0𝔑+1 with respect to 𝒢. Note that these basis functions have a proper support. See [9, 18] for detail. Let 𝒫N,h0∶=H01(I)∩𝒫Nh and 𝒱N,h0∶=H01(I)∩𝒱Nh. With the notation 𝒱×𝒱∶={[u,v]T|u,v∈𝒱}, define P̲N,h0 and V̲N,h0 as subspaces of 𝒫N,h0×𝒫N,h0 and 𝒱N,h0×𝒱N,h0, respectively. The basis functions for P̲N,h0 and V̲N,h0 are given respectively by(2.11)ϕ̲p(t)∶={[ϕp(t),0]T,forp≤N,[0,ϕp-N(t)]T,forp>N,ψ̲p(t)∶={[ψp(t),0]T,forp≤N,[0,ψP-N(t)]T,forp>N,
where p=1,2,…,2𝔑. For two dimensional (2D) case, let [[𝒫N,h0]]∶=𝒫N,h0⊗𝒫N,h0 and [[𝒱N,h0]]∶=𝒱N,h0⊗𝒱N,h0 be tensor product function spaces of one-dimensional function spaces and let [[P̲N,h0]] and [[V̲N,h0]] be subspaces of [[𝒫N,h0]]×[[𝒫N,h0]] and [[𝒱N,h0]]×[[𝒱N,h0]], respectively. Now let us order the interior LGL points in Ω by horizontal lines as {Ξμ̃}μ̃=1𝔑2∶={(ξμ,ξν)}μ=1,ν=1𝔑,𝔑, where μ̃=μ+𝔑(ν-1) for μ,ν=1,2,…,𝔑. Accordingly, the basis functions for [[P̲N,h0]] and [[V̲N,h0]] 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̲N,h0]] and [[V̲N,h0]] are given, respectively, by(2.12)Φ̲p̃(x,y)∶={[ϕp̃(x,y),0]T,forp̃≤N2,[0,ϕp̃-N2(x,y)]T,forp̃>N2,Ψ̲p̃(x,y)∶={[ψp̃(x,y),0]T,forp̃≤N2,[0,ψp̃-N2(x,y)]T,forp̃>N2,
where p̃=1,…,2𝔑2.


### 2.3. Interpolation Operators

We denoteC(Ω) as the set of continuous functions in Ω∶=I×I. Let ℐN2:C(Ω)→[[𝒫N]]∶=𝒫N⊗𝒫N be the usual reference interpolation operator such that (ℐN2u)(ηi,ηj)=u(ηi,ηj) for u∈C(Ω) (see, e.g., [17]). The global interpolation operator ℐN2h:C(Ω)→[[𝒫N,h0]] is given by ℐN2hu(x,y)=∑μ̃=1𝔑2uμ̃ϕμ̃(x,y), where uμ̃=u(ξμ,ξν). Hence, it follows that (ℐN2hu)(ξμ,ξν)=u(ξμ,ξν) for u∈C(Ω). With this interpolation operator ℐN2h, let us define the vector interpolation operator I̲N2h:C(Ω)×C(Ω)→[[P̲N,h0]] such that, for u̲∶=[u,v]T∈C(Ω)×C(Ω),(2.13)(I̲N2hu̲)(ξμ,ξν)∶=[IN2hu(ξμ,ξν),IN2hv(ξμ,ξν)]T=u̲(ξμ,ξν).
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 [[𝒫N]] and [[𝒱N]], respectively. Let us denote MN2 and Mh2 by the mass matrices such that(2.14)MN2(ĩ,j̃)=〈Φ̂ĩ,Φ̂j̃〉,Mh2(ĩ,j̃)=〈Ψ̂ĩ,Ψ̂j̃〉
and denote SN2 and Sh2 by the stiffness matrices such that(2.15)SN2(ĩ,j̃)=〈∇Φ̂ĩ,∇Φ̂j̃〉,Sh2(ĩ,j̃)=〈∇Ψ̂ĩ,∇Ψ̂j̃〉,
where ĩ,j̃=1,…,(N+1)2.According to Theorems 5.4 and 5.5 in [17], there are two absolute positive constants c0 and c1 such that for any U=[u1,…,u(N+1)2]T,(2.16)c0〈MN2U,U〉≤〈Mh2U,U〉≤c1〈MN2U,U〉,c0〈SN2U,U〉≤〈Sh2U,U〉≤c1〈SN2U,U〉,
and for all u∈[[𝒱N]],(2.17)c0‖u‖≤‖IN2u‖≤c1‖u‖,c0‖u‖1≤‖IN2u‖1≤c1‖u‖1.
The extension of (2.17) to the interpolation operator ℐN2h leads to(2.18)c0‖u‖≤‖IN2hu‖≤c1‖u‖,c0‖u‖1≤‖IN2hu‖1≤c1‖u‖1
for all u∈[[𝒱N,h0]], where the constants c0 and c1 are positive constants independent of E and N (see [18]).Theorem 2.1.
For allu̲∈[[V̲N,h0]], there are positive constants c0 and c1 independent of E and N such that
(2.19)c0‖u̲‖≤‖I̲N2hu̲‖≤c1‖u̲‖,c0‖u̲‖1≤‖I̲N2hu̲‖1≤c1‖u̲‖1.Proof.
By the definitions of the interpolation operatorI̲N2h and the norms, we have
(2.20)‖I̲N2hu̲‖2=‖IN2hu‖2+‖IN2hv‖2,‖I̲N2hu̲‖12=‖IN2hu‖12+‖IN2hv‖12,
which completes the proof because of (2.18).


## 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 functionu̲=[u,v]T, where T denotes the transpose, are defined by(2.1)Δu̲∶=[uxx+uyyvxx+vyy],∇u̲∶=[uxvxuyvy],∇⋅u̲∶=[ux+uyvx+vy].
With the usual L2 inner product 〈·,·〉 and its norm ∥·∥, for vector functions u̲∶=[u,v]T and w̲∶=[w,z]T, define(2.2)〈u̲,w̲〉∶=〈u,w〉+〈v,z〉,‖u̲‖2∶=‖u‖2+‖v‖2
and, for matrix functions U and V, define(2.3)〈U,V〉=∑k=14〈uk,vk〉,‖U‖2∶=∑k=14‖uk‖2,whereU=[u1u2u3u4],V=[v1v2v3v4].We use the standard Sobolev spaces likeH1(Ω) and H01(Ω) 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 by1/δ, the system (1.1) can be expressed as(2.4)Au̲∶=-AΔu̲+B(b⋅∇u̲)T+(A+C)u̲=f̲inΩ,
where(2.5)A=[1001δ],B=[100-1δ],C=[01δ-1δ0],f̲=[0-1δû],
with the zero boundary condition u̲=0̲on∂Ω. Let ℬ be another decoupled uniformly elliptic operator such that(2.6)Bu̲∶=-AΔu̲+Au̲inΩ
with the zero boundary condition.


## 2.2. LGL Nodes, Weights, and Function Spaces

Let{ηk}k=0N and {ωk}k=0N be the reference LGL nodes and its corresponding LGL weights in I=[-1,1], respectively, arranged by -1=:η0<η1<⋯<ηN-1<ηN∶=1. We use {tj}j=0E as the set of knots in the interval I such that -1=:t0<t1<⋯<tE-1<tE∶=1. Here E denotes the number of subintervals of I. Denote N by the degree of a polynomial on each subinterval Ij∶=[tj-1,tj] and 𝒢∶={ξj,k}j=1,k=0E,N by the set of kth-LGL nodes ξj,k in each subinterval Ij (j=1,…,E) arranged by(2.7)ξj,0∶=tj-1<ξj,1<⋯<ξj,N-1<tj=:ξj,N,
where(2.8)ξj,k=hj2ηk+12(tj-1+tj),hj=tj-tj-1,
and the corresponding LGL weights {ρj,k}j=1,k=0E,N are given by(2.9)ρj,k∶={hj2ωN+hj+12ω0,k=N,j=1,…,E-1,ρj-1,N,k=0,j=2,…,E,hj2ωk,otherwise.
Let 𝒫k be the space of all polynomials defined on I whose degrees are less than or equal to k. The Lagrange basis for 𝒫N on I is given by {ϕ̂i(t)}i=0N satisfying(2.10)ϕ̂i(ηj)=δijfori,j=0,1,…,N,
where δij denotes the Kronecker delta function.We define𝒫Nh as the subspace of continuous functions whose basis {ϕμ}μ=0𝔑+1 is piecewise continuous Lagrange polynomials of degree N on Ij with respect to 𝒢 and define 𝒱Nh as the space of all piecewise Lagrange linear functions {ψμ}μ=0𝔑+1 with respect to 𝒢. Note that these basis functions have a proper support. See [9, 18] for detail. Let 𝒫N,h0∶=H01(I)∩𝒫Nh and 𝒱N,h0∶=H01(I)∩𝒱Nh. With the notation 𝒱×𝒱∶={[u,v]T|u,v∈𝒱}, define P̲N,h0 and V̲N,h0 as subspaces of 𝒫N,h0×𝒫N,h0 and 𝒱N,h0×𝒱N,h0, respectively. The basis functions for P̲N,h0 and V̲N,h0 are given respectively by(2.11)ϕ̲p(t)∶={[ϕp(t),0]T,forp≤N,[0,ϕp-N(t)]T,forp>N,ψ̲p(t)∶={[ψp(t),0]T,forp≤N,[0,ψP-N(t)]T,forp>N,
where p=1,2,…,2𝔑. For two dimensional (2D) case, let [[𝒫N,h0]]∶=𝒫N,h0⊗𝒫N,h0 and [[𝒱N,h0]]∶=𝒱N,h0⊗𝒱N,h0 be tensor product function spaces of one-dimensional function spaces and let [[P̲N,h0]] and [[V̲N,h0]] be subspaces of [[𝒫N,h0]]×[[𝒫N,h0]] and [[𝒱N,h0]]×[[𝒱N,h0]], respectively. Now let us order the interior LGL points in Ω by horizontal lines as {Ξμ̃}μ̃=1𝔑2∶={(ξμ,ξν)}μ=1,ν=1𝔑,𝔑, where μ̃=μ+𝔑(ν-1) for μ,ν=1,2,…,𝔑. Accordingly, the basis functions for [[P̲N,h0]] and [[V̲N,h0]] 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̲N,h0]] and [[V̲N,h0]] are given, respectively, by(2.12)Φ̲p̃(x,y)∶={[ϕp̃(x,y),0]T,forp̃≤N2,[0,ϕp̃-N2(x,y)]T,forp̃>N2,Ψ̲p̃(x,y)∶={[ψp̃(x,y),0]T,forp̃≤N2,[0,ψp̃-N2(x,y)]T,forp̃>N2,
where p̃=1,…,2𝔑2.


## 2.3. Interpolation Operators

We denoteC(Ω) as the set of continuous functions in Ω∶=I×I. Let ℐN2:C(Ω)→[[𝒫N]]∶=𝒫N⊗𝒫N be the usual reference interpolation operator such that (ℐN2u)(ηi,ηj)=u(ηi,ηj) for u∈C(Ω) (see, e.g., [17]). The global interpolation operator ℐN2h:C(Ω)→[[𝒫N,h0]] is given by ℐN2hu(x,y)=∑μ̃=1𝔑2uμ̃ϕμ̃(x,y), where uμ̃=u(ξμ,ξν). Hence, it follows that (ℐN2hu)(ξμ,ξν)=u(ξμ,ξν) for u∈C(Ω). With this interpolation operator ℐN2h, let us define the vector interpolation operator I̲N2h:C(Ω)×C(Ω)→[[P̲N,h0]] such that, for u̲∶=[u,v]T∈C(Ω)×C(Ω),(2.13)(I̲N2hu̲)(ξμ,ξν)∶=[IN2hu(ξμ,ξν),IN2hv(ξμ,ξν)]T=u̲(ξμ,ξν).
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 [[𝒫N]] and [[𝒱N]], respectively. Let us denote MN2 and Mh2 by the mass matrices such that(2.14)MN2(ĩ,j̃)=〈Φ̂ĩ,Φ̂j̃〉,Mh2(ĩ,j̃)=〈Ψ̂ĩ,Ψ̂j̃〉
and denote SN2 and Sh2 by the stiffness matrices such that(2.15)SN2(ĩ,j̃)=〈∇Φ̂ĩ,∇Φ̂j̃〉,Sh2(ĩ,j̃)=〈∇Ψ̂ĩ,∇Ψ̂j̃〉,
where ĩ,j̃=1,…,(N+1)2.According to Theorems 5.4 and 5.5 in [17], there are two absolute positive constants c0 and c1 such that for any U=[u1,…,u(N+1)2]T,(2.16)c0〈MN2U,U〉≤〈Mh2U,U〉≤c1〈MN2U,U〉,c0〈SN2U,U〉≤〈Sh2U,U〉≤c1〈SN2U,U〉,
and for all u∈[[𝒱N]],(2.17)c0‖u‖≤‖IN2u‖≤c1‖u‖,c0‖u‖1≤‖IN2u‖1≤c1‖u‖1.
The extension of (2.17) to the interpolation operator ℐN2h leads to(2.18)c0‖u‖≤‖IN2hu‖≤c1‖u‖,c0‖u‖1≤‖IN2hu‖1≤c1‖u‖1
for all u∈[[𝒱N,h0]], where the constants c0 and c1 are positive constants independent of E and N (see [18]).Theorem 2.1.
For allu̲∈[[V̲N,h0]], there are positive constants c0 and c1 independent of E and N such that
(2.19)c0‖u̲‖≤‖I̲N2hu̲‖≤c1‖u̲‖,c0‖u̲‖1≤‖I̲N2hu̲‖1≤c1‖u̲‖1.Proof.
By the definitions of the interpolation operatorI̲N2h and the norms, we have
(2.20)‖I̲N2hu̲‖2=‖IN2hu‖2+‖IN2hv‖2,‖I̲N2hu̲‖12=‖IN2hu‖12+‖IN2hv‖12,
which completes the proof because of (2.18).


## 3. Analysis on𝒫1 Finite Element Preconditioner

The bilinear forms corresponding to (2.4) and (2.6) are given by(3.1)α(u̲,w̲)=〈A∇u̲T,∇w̲T〉+〈B(b⋅∇u̲)T,w̲〉+〈(A+C)u̲,w̲〉=〈f̲,w̲〉,(3.2)β(u̲,w̲)=〈A∇u̲T,∇w̲T〉+〈Au̲,w̲〉,
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 H01(Ω)×H01(Ω) for the variational problem (3.1).Proposition 3.1.
For a real valued vector functionu̲(x,y)=[u,v]T∈H01(Ω)×H01(Ω), we have
(3.3)‖u̲‖1≤α(u̲,u̲)≤1δ‖u̲‖12.Proof.
Sinceb is a constant vector, [u,v]T∈H01(Ω)×H01(Ω), and u̲ is a real valued vector function, we have
(3.4)α(u̲,u̲)=β(u̲,u̲)+〈B(b⋅∇u̲)T,u̲〉+〈Cu̲,u̲〉=β(u̲,u̲)+〈b⋅∇u,u〉-1δ〈b⋅∇v,v〉+1δ(〈v,u〉-〈u,v〉)=β(u̲,u̲)=‖u‖12+1δ‖v‖12.
Hence, (3.3) is proved because 0<δ≤1.Lemma 3.2.
Ifu̲=[u,v]T is a complex valued function in H01(Ω)×H01(Ω), then we have the following estimates:
(3.5)Re(α(u̲,u̲))=β(u̲,u̲),|α(u̲,u̲)|≤(1+12|b|+1δ)β(u̲,u̲),
where |b|=b12+b22.Proof.
Let us decomposeu 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 i2=-1. Then we have
(3.6)〈b⋅∇u,u〉=i∫Ωb1(qxp-pxq)+b2(pyq-qyp)dΩ,(3.7)-1δ〈b⋅∇u,u〉=-iδ∫Ωb1(sxr-rxs)+b2(rys-syr)dΩ,(3.8)1δ(〈v,u〉-〈u,v〉)=i(2δ∫Ωsp-qrdΩ).
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 0<δ≤1, we have
(3.9)|〈b⋅∇u,u〉-1δ〈b⋅∇v,v〉+1δ(〈v,u〉-〈u,v〉)|≤|b|‖∇u‖‖u‖+1δ|b|‖∇v‖‖v‖+1δ(2‖v‖‖u‖)≤12|b|(‖∇u‖2+‖u‖2+1δ(‖∇v‖2+‖v‖2))+1δ(‖u‖2+‖v‖2)≤(12|b|+1δ)β(u̲,u̲).
Hence (3.5) is proved.Letσ(A) and ℱ(A) be the spectrum (or set of eigenvalues) and field of values of the square matrix A, respectively. Let ÂN2 and B̂h2 be the two dimensional stiffness matrices on the spaces [[P̲N,h0]] and [[V̲N,h0]] induced by (3.1) and (3.2), respectively. Then, we have the following.Lemma 3.3.
ForU(≠0)∈ℂ2𝔑2, one has
(3.10)σ(B̂h2-1ÂN2)⊂W,
where
(3.11)W∶={〈ÂN2U,U〉〈B̂h2U,U〉∣U≠0}.Proof.
SinceB̂h2 is symmetric positive definite, there exists a unique positive definite square root B̂h21/2 of B̂h2. So, we have
(3.12)〈ÂN2U,U〉〈B̂h2U,U〉=〈ÂN2U,U〉〈B̂h21/2U,B̂h21/2U〉=〈B̂h2-1/2ÂN2B̂h2-1/2V,V〉〈V,V〉,whereV=B̂h21/2U,
for U(≠0)∈ℂ2𝔑2. Now let B̂h2-1/2 be a short-hand notation for (B̂h21/2)-1. Therefore, from the relation of spectrum and field of values (see [19] or [11]), it follows that
(3.13)σ(B̂h2-1ÂN2)=σ(B̂h2-1/2ÂN2B̂h2-1/2)⊂F(B̂h2-1/2ÂN2B̂h2-1/2)=W,
which completes the proof.The following theorem shows the uniform bounds of eigenvalues which is independent of bothN and E for our preconditioned system(3.14)B̂h2-1ÂN2u̲N=B̂h2-1f̲N.Theorem 3.4.
Let{λp}p=12𝔑2 be the set of eigenvalues of
(3.15)B̂h2-1ÂN2,
then, there are constants c0,C0, and Λδ independent of E and N, such that
(3.16)0<c0<Re(λp)<C0,|λp|≤Λδ,
where Λδ=C(1+∥b∥+1/δ).Proof.
Letu̲(x,y)∈[[V̲N,h0]] be represented as u̲(x,y)=∑p=12𝔑2upΨ̲p(x,y). Then, its piecewise polynomial interpolation can be written as
(3.17)(I̲N2hu̲)(x,y)=∑p=12N2upΦ̲p(x,y).
Let U=[u1,…,u2𝔑2]T. From the definitions of the bilinear forms, we have
(3.18)〈ÂN2U,U〉=α(I̲N2hu̲,I̲N2hu̲),〈B̂h2U,U〉=β(u̲,u̲).
This implies that
(3.19)wλ∶=〈ÂN2U,U〉〈B̂h2U,U〉=α(I̲N2hu̲,I̲N2hu̲)β(u̲,u̲).
By Theorem 2.1 and Lemma 3.2, the real part of wλ satisfies
(3.20)Re(wλ)=β(I̲N2hu̲,I̲N2hu̲)β(u̲,u̲)~1,
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)|wλ|=|α(I̲N2hu̲,I̲N2hu̲)|β(u̲,u̲)≤(1+(1/2)|b|+(1/δ))β(I̲N2hu̲,I̲N2hu̲)β(u̲,u̲)≤C(1+|b|+1δ).
Therefore, from Lemma 3.3, the real parts and the absolute values of eigenvalues of B̂h2-1ÂN2 satisfy (3.16).Remark 3.5.
Let the one dimensional (1D) bilinear forms foru̲,w̲∈H01(-1,1)×H01(-1,1) be
(3.22)α1(u̲,w̲)=〈Au̲′,w̲′〉+〈Bbu̲′,w̲〉+〈(A+C)u̲,w̲〉,(3.23)β1(u̲,w̲)=〈Au̲′,w̲′〉+〈Au̲,w̲〉,
and denote ÂN and B̂h by the 1D stiffness matrices on the spaces P̲N,h0 and V̲N,h0 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)SN=〈ϕμ′(t),ϕν′(t)〉,RN=〈ϕμ′(t),ϕν(t)〉,MN=〈ϕμ(t),ϕν(t)〉
and 1D finite element matrices(3.25)Sh=〈ψμ′(t),ψν′(t)〉,Rh=〈ψμ′(t),ψν(t)〉,Mh=〈ψμ(t),ψν(t)〉,
where {ϕμ}μ=1𝔑 and {ψν}ν=1𝔑 are the Lagrange basis of the spaces 𝒫N,h0 and 𝒱N,h0, respectively. For actual computations for SN, RN, and MN, the inner product 〈·,·〉 on the space 𝒫N,h0 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 SN, RN, and MN and the exact matrices SN, RN, and MN are equivalent, respectively, because of the equivalence of LGL quadrature on the polynomial space we used. For example, the mass matrix MN can be computed using the LGL weights {ρj,k} 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 ÂN and B̂h corresponding to (3.22) and (3.23) using the matrices in (3.24). One may have(4.1)ÂN=[SN+bRN+MN1δMN-1δMN1δ(SN-bRN+MN)]2N×2N,B̂h=[Sh+Mh001δ(Sh+Mh)]2N×2N.
Let b=[b1,b2]T be a constant vector in (1.1), then the 2D matrices ÂN2 and B̂h2 can be expressed as(4.2)ÂN2=[SN2+RN2+MN21δMN2-1δMN21δ(SN2-RN2+MN2)]2N2×2N2,B̂h2=[Sh2+Mh2001δ(Sh2+Mh2)]2N2×2N2,
where(4.3)Sk2=Mk(y)⊗Sk(x)+Sk(y)⊗Mk(x),Mk2=Mk(y)⊗Mk(x)(k=Norh),RN2=b1MN(y)⊗RN(x)+b2RN(y)⊗MN(x).


### 4.2. Numerical Analysis on Eigenvalues

The linear systemÂN2u̲N=f̲N and the preconditioned linear system B̂h2-1ÂN2u̲N=B̂h2-1f̲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̂h2-1ÂN2 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̂N2 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.The condition numbers of the matrixÂN2 when b=[1,1]T.
(a)(b)Figures2 and 4 show the spectra of the resulting preconditioned operator B̂h2-1ÂN2 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̂h2-1ÂN2 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̂h2-1ÂN2 are independent of N and E, but they depend still on the penalty parameter δ.The eigenvalues ofB̂h2-1ÂN2 for N=4,6,8,10 and E=4 when δ=1,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for E=4,6,8,10 and N=4 when δ=1,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for N=4,6,8,10 and E=4 when δ=10-4,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for E=4,6,8,10 and N=4 when δ=10-4,b=[1,1]T.
(a)(b)(c)(d)By choosing the convection coefficientb=[10,10]T, in Figure 6, the distributions of eigenvalues of ÂN2(left) and B̂h2-1ÂN2(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̂N2 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̂h2-1ÂN2 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.The eigenvalues ofÂN2(left) and B̂h2-1ÂN2(right) for δ=1 (top) and δ=10-3 (bottom) when b=[10,10]T and N=12,E=4.
(a)(b)(c)(d)


## 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 ÂN and B̂h corresponding to (3.22) and (3.23) using the matrices in (3.24). One may have(4.1)ÂN=[SN+bRN+MN1δMN-1δMN1δ(SN-bRN+MN)]2N×2N,B̂h=[Sh+Mh001δ(Sh+Mh)]2N×2N.
Let b=[b1,b2]T be a constant vector in (1.1), then the 2D matrices ÂN2 and B̂h2 can be expressed as(4.2)ÂN2=[SN2+RN2+MN21δMN2-1δMN21δ(SN2-RN2+MN2)]2N2×2N2,B̂h2=[Sh2+Mh2001δ(Sh2+Mh2)]2N2×2N2,
where(4.3)Sk2=Mk(y)⊗Sk(x)+Sk(y)⊗Mk(x),Mk2=Mk(y)⊗Mk(x)(k=Norh),RN2=b1MN(y)⊗RN(x)+b2RN(y)⊗MN(x).


## 4.2. Numerical Analysis on Eigenvalues

The linear systemÂN2u̲N=f̲N and the preconditioned linear system B̂h2-1ÂN2u̲N=B̂h2-1f̲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̂h2-1ÂN2 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̂N2 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.The condition numbers of the matrixÂN2 when b=[1,1]T.
(a)(b)Figures2 and 4 show the spectra of the resulting preconditioned operator B̂h2-1ÂN2 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̂h2-1ÂN2 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̂h2-1ÂN2 are independent of N and E, but they depend still on the penalty parameter δ.The eigenvalues ofB̂h2-1ÂN2 for N=4,6,8,10 and E=4 when δ=1,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for E=4,6,8,10 and N=4 when δ=1,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for N=4,6,8,10 and E=4 when δ=10-4,b=[1,1]T.
(a)(b)(c)(d)The eigenvalues ofB̂h2-1ÂN2 for E=4,6,8,10 and N=4 when δ=10-4,b=[1,1]T.
(a)(b)(c)(d)By choosing the convection coefficientb=[10,10]T, in Figure 6, the distributions of eigenvalues of ÂN2(left) and B̂h2-1ÂN2(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̂N2 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̂h2-1ÂN2 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.The eigenvalues ofÂN2(left) and B̂h2-1ÂN2(right) for δ=1 (top) and δ=10-3 (bottom) when b=[10,10]T and N=12,E=4.
(a)(b)(c)(d)


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



---

*Source: 245051-2012-02-12.xml*