Spectral-Galerkin Approximation and Optimal Error Estimate for Stokes Eigenvalue Problems in Polar Geometries This work is supported in part by the National Natural Science Foundation of China grants No. 11661022, 91130014, 11471312, 91430216, 11471031, and U1530401; and by the US National Science Foundation grant DMS-1419040.

# Spectral-Galerkin Approximation and Optimal Error Estimate for Stokes Eigenvalue Problems in Polar Geometries††thanks: This work is supported in part by the National Natural Science Foundation of China grants No. 11661022, 91130014, 11471312, 91430216, 11471031, and U1530401; and by the US National Science Foundation grant DMS-1419040.

Jing An Beijing Computational Science Research Center, Beijing 100193, China (anjing@csrc.ac.cn); and School of Mathematical Sciences, Guizhou Normal University, Guiyang 550025, China (aj154@163.com).    Huiyuan Li State Key Laboratory of Computer Science/Laboratory of Parallel Computing, Institute of Software, Chinese Academy of Sciences, Beijing 100190, China (huiyuan@iscas.ac.cn).    Zhimin Zhang Beijing Computational Science Research Center, Beijing 100193, China (zmzhang@csrc.ac.cn); and Department of Mathematics, Wayne State University, Detroit, MI 48202, USA (zzhang@math.wayne.edu).
###### Abstract

In this paper we propose and analyze spectral-Galerkin methods for the Stokes eigenvalue problem based on the stream function formulation in polar geometries. We first analyze the stream function formulated fourth-order equation under the polar coordinates, then we derive the pole condition and reduce the problem on a circular disk to a sequence of equivalent one-dimensional eigenvalue problems that can be solved in parallel. The novelty of our approach lies in the construction of suitably weighted Sobolev spaces according to the pole conditions, based on which, the optimal error estimate for approximated eigenvalue of each one dimensional problem can be obtained. Further, we extend our method to the non-separable Stokes eigenvalue problem in an elliptic domain and establish the optimal error bounds. Finally, we provide some numerical experiments to validate our theoretical results and algorithms.

Keywords: Stokes eigenvalue problem, polar geometry, pole condition, spectral-Galerkin approximation, optimal error analysis

## 1 Introduction

We consider in this paper the Stokes eigenvalue problem which arises in stability analysis of the stationary solution of the Navier-Stokes equations [20]:

 −Δu+∇p=λu, in  Ω, (1.1) ∇⋅u=0, in  Ω, (1.2) u=0, on  ∂Ω, (1.3)

where is the flow velocity, is the pressure, is the Laplacian operator, is the flow domain and denotes the boundary of the flow domain .

Let us introduce the stream function such that . Then we derive an alternative formulation for (1.1)-(1.3):

 −Δ2ψ=λΔψ, in  Ω, (1.4) ψ=∂ψ∂n=0, on  ∂Ω, (1.5)

where is the unit outward normal to the boundary . (1.4) is also referred to as the biharmonic eigenvalue problem for plate buckling. The naturally equivalent weak form of (1.4)-(1.5) reads: Find such that

 A(ψ,ϕ)=λB(ψ,ϕ),ϕ∈H20(Ω), (1.6)

where the bilinear forms and are defined by

 A(ψ,ϕ)=(Δψ,Δϕ)=∫ΩΔψΔ¯¯¯ϕdxdy, B(ψ,ϕ)=(∇ψ,∇ϕ)=∫Ω∇ψ⋅∇¯¯¯ϕdxdy.

There are various numerical approaches to solving (1.4)-(1.5). Mixed finite element methods introduce the auxiliary function to reduce the fourth-order equation to a saddle point problem and then discretize the reduced second order equations with (-) continuous finite elements[8, 22, 10, 29]. However, spurious solutions may occur in some situations. The conforming finite element methods including Argyris elements [2] and the partition of unity finite elements [11], require globally continuously differentiable finite element spaces, which are difficult to construct and implement. The third type of approaches use non-conforming finite element methods, such as Adini elements [1], Morley elements [19, 21, 25] and the ordinary -interior penalty Galerkin method [26]. Their disadvantage lies in that such elements do not come in a natural hierarchy. Both the conforming and nonconforming finite element methods are based on the naturally equivalent variational formulation (1.6), and usually involve low order polynomials and guarantee only a low order of convergence.

In contrast, it is observed in [31] that the spectral method, whenever it is applicable, has tremendous advantage over the traditional -version methods. In particular, spectral and spectral element methods using high order orthogonal polynomials for fourth-order equations result in an exponential order of convergence for smooth solutions [23, 6, 5, 13, 30, 14, 9]. In analogy to the Argyris finite element methods, the conforming spectral element method requires globally continuously differentiable element spaces, which are extremely difficult to construct and implement on unstructured (triangular or quadrilateral) meshes. This is exactly the reason why -conforming spectral elements are rarely reported in literature except those on rectangular meshes [30]. Hence, the spectral methods using globally smooth basis functions are naturally suitable choices in practice for (1.6) on some fundamental regions including rectangles, triangles and polar geometries.

To the best of our knowledge there are few reports on spectral-Galerkin approximation for the Stokes eigenvalue problem by the stream function formulation in polar geometries. The polar transformation introduces polar singularities and variable coefficients of the form in polar coordinates [23, 4], which involves intricate pole conditions thus brings forth severe difficulties in both the design of approximation schemes and the corresponding error analysis. The aim of the current paper is to propose and analyze an efficient spectral-Galerkin approximation for the stream function formulation of the Stokes eigenvalue problem in polar geometries. As the first step, we use the separation of variables in polar coordinates to reduce the original problem in the unit disk to equivalent infinite sequence of one-dimensional eigenvalue problems which can be solved individually in parallel. Rigorous pole conditions involved are prerequisite for the equivalence of the original problem and the sequence of the one-dimensional eigenvalue problems, and thus play a fundamental role in our further study. It is worthy to note, however, that the pole conditions derived for the fourth-order source problems in open literature (such as [23, 4]) are inadequate for our eigenvalue problems since they would inevitably induce improper/spurious computational results.

Based on the pole condition, suitable approximation spaces are introduced and spectral-Galerkin schemes are proposed. A rigorous analysis on the optimal error estimate in certain properly introduced weighted Sobolev spaces is made for each one dimensional eigenvalue problem by using the minimax principle. Finally, we extend our spectral-Galerkin method to solving the stream function formulation of the Stokes eigenvalue problem in an elliptic region. Owing to its non-separable property, this problem is actually another challenge both in computation and analysis. A brief explanation on the implementation of the approximation scheme is first given, and an optimal error estimate is then presented in the Cartesian coordinates under the framework of Babǔska and Osborn [3].

The rest of this paper is organized as follows. In the next section, dimension reduction scheme of the Stokes eigenvalue problem is presented. In §3, we derive the weak formulation and prove the error estimation for a sequence of equivalent one-dimensional eigenvalue problems. Also, we describe the details for an efficient implementation of the algorithm. In §4, we extend our algorithm to the case of elliptic region. We present several numerical experiments in §5 to demonstrate the accuracy and efficiency of our method. Finally, in §6 we give some concluding remarks.

## 2 Dimensionality reduction and pole conditions

Before coming to the main body of this section, we would like to introduce some notations and conventions which will be used throughout the paper. Let be a generic positive weight function on a bounded domain , which is not necessarily in . Denote by the inner product of whose norm is denoted by . We use and to denote the usual weighted Sobolev spaces, whose norm is denoted by . In cases where no confusion would arise, (if ) and may be dropped from the notation. Let (resp. ) be the collection of nonnegative integers (resp. integers). For , we denote by the collection of all algebraic polynomials on with the total degree no greater than . We denote by a generic positive constant independent of any function and of any discretization parameters. We use the expression to mean that .

In the current section, we restrict our attention to the unit disk . We shall employ a classical technique, separation of variables, to reduce the problem to a sequence of equivalent one-dimensional problems.

Throughout this paper, we shall use the polar coordinates for points in the disk such that . We associate any function in Cartesian coordinates with its partner in polar coordinates. If no confusion would arise, we shall use the same notation for and . We now recall that, under the polar coordinates,

 Δ=1r∂∂r(r∂∂r)+1r2∂2∂θ2,∇=(cosθ∂∂r−sinθr∂∂θ,sinθ∂∂r+cosθr∂∂θ)t. (2.1)

Then the bilinear forms and in (1.6) become

 A(ψ,ϕ)=∫10rdr∫2π0[∂2ψ∂r2+1r∂ψ∂r+1r2∂2ψ∂θ2][∂2¯¯¯ϕ∂r2+1r∂¯¯¯ϕ∂r+1r2∂2¯¯¯ϕ∂θ2]dθ, B(ψ,ϕ)=∫10rdr∫2π0[∂ψ∂r∂¯¯¯ϕ∂r+1r2∂ψ∂θ∂¯¯¯ϕ∂θ]dθ.

Denote and define the bilinear forms for functions on ,

 Bm(u,v)=∫10(ru′¯¯¯v′+m2ru¯¯¯v)dr.

Further let us assume

 ψ=∑m∈Zψm(r)eimθ,ϕ=∑m∈Zϕm(r)eimθ. (2.2)

By the orthogonality of the Fourier system , one finds that

 A(ψ,ϕ)=2π∑m∈ZAm(ψm,ϕm),B(ψ,ϕ)=2π∑m∈ZBm(ψm,ϕm).

For the well-posedness of and , the following pole conditions for (and the same type of pole conditions for ) should be imposed,

 mψm(0)=0,limr→0+[ψ′m(r)−m2rψm(r)]=(1−m2)ψ′m(0)=0, (2.3)

which can be further simplified into the following three categories,

 (1).ψ′m(0)=0, m=0; (2.4) (2).ψm(0)=0, |m|=1; (2.5) (3).ψm(0)=ψ′m(0)=0, |m|≥2. (2.6)

It is worthy to note that our pole condition (2.5) for is a revision of the pole condition in (4.8) of [23]. A concrete example to support the absence of reads,

 ψ=ψ±1(r)e±iθ∈H20(D),ψ±1(r)=(1−r)2r.

Also, this absence of in (2.5) is also confirmed by [7].

The boundary conditions on states for all integer . Meanwhile, together with implies . It is then easy to verify that (resp.  ) induces a Sobolev norm for any function on which satisfies the boundary condition (resp. ) and the pole condition (resp. ).

We now introduce two non-uniformly weighted Sobolev spaces on ,

 \lx@overaccentset∘H1m(I):={u: Bm(u,u)<∞, mu(0)=u(1)=0}, (2.7) \lx@overaccentset∘H2m(I):={u: Am(u,u)<∞, mu(0)=(1−m2)u′(0)=u(1)=u′(1)=0}, (2.8)

which are endowed with energy norms

 (2.9)

In the sequel, (1.6) is reduced to a system of infinite one-dimensional eigen problems: to find such that and

 Am(ψm,ϕm)=λmBm(ψm,ϕm),ϕm∈\lx@overaccentset∘H2m(I),m∈Z. (2.10)

We now conclude this section with the following lemma on and .

###### Lemma 2.1

For ,

 Bm(u,v)=∫10(u′±mru)(¯¯¯v′±mr¯¯¯v)rdr, (2.11) Am(u,v)=∫10[r(u′∓mru)′(¯¯¯v′∓mr¯¯¯v)′+(1±m)2r(u′∓mru)(¯¯¯v′∓mr¯¯¯v)]dr. (2.12)

Proof. By integration by parts and the pole condition (2.3), one verifies that

 ∫10(u′± mru)(¯¯¯v′±mr¯¯¯v)rdr=∫10(ru′¯¯¯v′+m2ru¯¯¯v′)dr±m∫10(u¯¯¯v)′dr = ∫10(ru′¯¯¯v′+m2ru¯¯¯v′)dr,

which gives (2.11).

As a result,

 Am(u,v)= ∫10[(u′∓mru)′+1±mr(u′∓mru)][(¯¯¯v′∓mr¯¯¯v)′+1±mr(¯¯¯v′∓mr¯¯¯v)]rdr = ∫10[r(u′∓mru)′(¯¯¯v′∓mr¯¯¯v)′+(1±m)2r(u′∓mru)(¯¯¯v′∓mr¯¯¯v)]dr

Meanwhile, the pole conditions (2.4)-(2.6) states that both and vanish at the two endpoints of . Thus the last integral above is zero, and (2.12) is now proved.

## 3 Spectral Galerkin approximation and its error estimates

Let be the space of polynomials of degree less than or equal to on , and setting . Then the spectral Galerkin approximation scheme to (2.10) is: Find such that and

 Am(ψmN,vN)=λmNBm(ψmN,vN),∀vN∈XmN. (3.1)

Due to the symmetry properties and , we shall only consider from now on in this section.

### 3.1 Mini-max principle

To give the error analysis, we will use extensively the minimax principle.

###### Lemma 3.1

Let denote the eigenvalues of (2.10) and be any -dimensional subspace of . Then, for , there holds

 λlm=minVl⊂\lx@overaccentset∘H2m(I)maxv∈VlAm(v,v)Bm(v,v). (3.2)

Proof. See Theorem 3.1 in [18].

###### Lemma 3.2

Let denote the eigenvalues of (2.10) and be arranged in an ascending order, and define

 Ei,j=\rm span{ψim,⋯,ψjm},

where is the eigenfunction corresponding to the eigenvalue . Then we have

 λlm=maxv∈Ek,lAm(v,v)Bm(v,v) k≤l, (3.3) λlm=minv∈El,mAm(v,v)Bm(v,v) l≤m. (3.4)

Proof. See Lemma 3.2 in [18].

It is true that the minimax principle is also valid for the discrete formulation (3.1) (see [18]).

###### Lemma 3.3

Let denote the eigenvalues of (3.1), and be any -dimensional subspace of . Then, for , there holds

 λlmN=minVl⊂XmNmaxv∈VlAm(v,v)Bm(v,v). (3.5)

Define the orthogonal projection such that

 Am(ψm−Π2,mNψm,v)=0,∀v∈XmN. (3.6)
###### Theorem 3.1

Let be obtained by solving (3.1) as an approximation of , an eigenvalue of (2.10). Then, we have

 0<λlm≤λlmN≤λlmmaxv∈E1,lBm(v,v)Bm(Π2,mNv,Π2,mNv). (3.7)

Proof. According to the coerciveness of and we easily derive . Since , from (3.2) and (3.5) we can obtain . Let denote the space spanned by . It is obvious that is a -dimensional subspace of . From the minimax principle, we have

Since from and the non-negativity of , we have

 Am(Π2,mNv,Π2,mNv)≤Am(v,v).

Thus, we have

 λlmN ≤ maxv∈E1,lAm(v,v)Bm(Π2,mNv,Π2,mNv) = maxv∈E1,lAm(v,v)Bm(v,v)Bm(v,v)Bm(Π2,mNv,Π2,mNv) ≤ λlmmaxv∈E1,lBm(v,v)Bm(Π2,mNv,Π2,mNv).

The proof of Theorem 3.1 is completed.

### 3.2 Error estimates

Denote by the Jacobi weight function of index , which is not necessarily in . Define the -orthogonal projection such that

 (π0,0Nu−u,v)I=0,v∈PN(I).

Further, for , define recursively the -orthogonal projections such that

 [π−k,−kNu](r)=∫r0[π1−k,1−kN−1u′](t)dt+u(0).

Next, for any nonnegative integers , define the Sobolev space

 Hs,k(I)={u∈Hk(I):s∑l=0∥∂lru∥ωmax(l−k,0),max(l−k,0),I<∞}.

Now we have the following error estimate on .

###### Lemma 3.4 ([15, Theorem 3.1.4])

is a Legendre tau approximation of such that

 ∂lr[π−k,−kNu](0)=∂lru(0),∂lr[π−k,−kNu](1)=∂lru(1),0≤l≤k−1, (3.8) (π−k,−kNu−u,v)=0,v∈PN−2k. (3.9)

Further suppose with . Then for ,

 ∥∂lr(π−k,−kNu−u)∥ωl−k,l−k,I≲Nl−s∥∂sru∥ωs−k,s−k,I,0≤l≤k≤s. (3.10)
###### Theorem 3.2

Suppose and with and . Then for ,

 (3.11)

Proof. Define the differential operator and then set

 uN(r)=−1rm∫1rtm[π−1,−1N−1Dmu](t)dt.

We shall first prove . By (3.9), we find that

 ∫10 tm[π−1,−1N−1Dmu](t)dt=∫10tm[Dmu](t)dt=∫10∂t[tmu(t)]dt=0,N≥m+3,m≠0,

where the last equality sign is derived from the boundary condition . Moreover,

As a result, and

 uN(0)=−limr→01rm∫1rtm[π−1,−1N−1Dmu](t)dt=0,m≠0.

Further, implies

 [Dmu](1)=0,m∈Z;[Dmu](0)=0,m≠1,

which, together with the property (3.8) of , gives

 [π−1,−1N−1Dmu](1)=0,m∈Z;[π−1,−1N−1Dmu](0)=0,m≠1.

In the sequel, we deduce that if and . In summary, we conclude that .

Next by (2.12) and (3.10), we have

 ∥∥uN−u∥∥22,m= ≤ ≲ [N4−2s+(m−1)2N2−2s]∥∥∂s−1rDmu∥∥2ωs−2,s−2.

Finally, (3.11) is an immediate consequence of the projection theorem,

 ∥∥Π2,mNu−u∥∥2,m,I=infv∈XmN∥∥v−u∥∥2,m,I≤∥∥uN−u∥∥2,m,I.

The proof is now completed.

###### Theorem 3.3

Let is the -th approximate eigenvalue of . If with , then we have

 |λlmN−λlm|≲(N2+m2)N2−2smax1≤i≤l∥∥∂s−1r(∂r+mr)ψim∥∥2ωs−2,s−2,I.

Proof. For any , it can be represented by ; we then have

 Bm(v,v)−Bm(Π2,mNv,Π2,mNv)Bm(v,v)≤2|Bm(v,v−Π2,mNv)|Bm(v,v) ≤2∑li,j=1|μi||μj||Bm(ψim−Π2,mNψim,ψjm)|∑li=1|μi|2 ≤2lmaxi,j=1,⋯,l|Bm(ψim−Π2,mNψim,ψjm)|:=ε.

Meanwhile, by the variational form (2.10), the definition (3.6) of , Cauchy-Schwarz inequality and Theorem 3.2, we have

 |Bm(ψim−Π2,mNψim,ψjm)|=1λjm|λjmbm(ψjm,ψim−Π2,mNψim)| =1λjm|Am(ψjm,ψim−Π2,mNψim)|=1λjm|Am(ψjm−Π2,mNψjm,ψim−Π2,mNψim)| ≤1λjm∥ψjm−Π2,mNψjm∥2,m,I∥ψim−Π2,mNψim∥2,m,I ≲(N2+m4)N2−2s∥∂srψjm∥ωs−2,s−2,I∥∂srψim∥ωs−2,s−2,I.

As a result, we have the following estimate for ,

For sufficiently large , . Thus

 0

and we finally deduce from Theorem 3.1 that

 0<λlmN−λlm≤2λlmε≲(N2+m2)N2−2smax1≤i≤l∥∥∂s−1r(∂r+mr)ψim∥∥2ωs−2,s−2,I.

The proof is now completed.

### 3.3 Implementations

We describe in this section how to solve the problems (3.1) efficiently. To this end, we first construct a set of basis functions for . Let

 ϕi(r)=(1−r)2r2J2,1i−4(2r−1),i≥4, (3.12)

where is the Jacobi polynomial of degree .

It is clear that

 XmN=span{ϕmi=ϕi:4≤i≤N},m≥2 X0N=span{ϕ0i=ϕi:4≤i≤N}⊕span{ϕ03(r)=14(1−r)2(2r+1)}, X1N=span{ϕ1i=ϕi:4≤i≤N}⊕span{ϕ13(r)=12(1−r)2r}.

Define if and otherwise. Our basis functions lead to the penta-diagonal matrix and the deca-diagonal mass matrix instead of the hepta- and hendecagon-diagonal ones in [23].

###### Lemma 3.5

For ,

 ϕ′′i(r)=(i−3)(i−1)i2(2i−3)J0,1i−2(2r−1)+2(i−3)2(