A posteriori error estimates of mixed discontinuous Galerkin method for a class of Stokes eigenvalue problems

: For a class of Stokes eigenvalue problems including the classical Stokes eigenvalue problem and the magnetohydrodynamic Stokes eigenvalue problem a residual type a posteriori error estimate of the mixed discontinuous Galerkin finite element method using P k − P k − 1 element ( k ≥ 1) is studied in this paper. The a posteriori error estimators for approximate eigenpairs are given. The reliability and e ffi ciency of the posteriori error estimator for the eigenfunction are proved and the reliability of the estimator for the eigenvalue is also analyzed. The numerical results are provided to confirm the theoretical predictions and indicate that the method considered in this paper can reach the optimal convergence order O ( do f − 2 kd ).


Introduction
Stokes eigenvalue problem is of great significance because of its role in the stability analysis of fluid mechanics. Therefore, it is of great interest to study efficient numerical methods for solving this problem (e.g. see [1][2][3][4][5][6][7][8]).
The a posteriori error estimates of the classical Stokes eigenvalue problem based on the velocitypressure formulation received much attention from scholars. For example, the authors in [9][10][11] studied the low-order conforming mixed method and Jia et al. [12] and Sun et al. [6] explored the low-order non-conforming finite elements. Since the accuracy of low-order elements is not high Gedicke et al. [1] used the Arnold-Winther hybrid finite element method to analyze the a posteriori error estimation based on the stress-velocity formulation in R 2 and Gedicke et al. [2] adopted the divergence-conforming discontinuous Galerkin finite element method (DGFEM for short) to discuss the a posteriori error estimate for velocity-pressure formulation on shape-regular rectangular meshes. Lepe et al. [13] proposed a mixed numerical method to study the error estimates for a vorticity-based velocity-stress formulation of the Stokes eigenvalue problem.
In this paper, for a class of Stokes eigenvalue problems (see (2.1)), including the classical Stokes eigenvalue problem in R d (d = 2, 3) and the magnetohydrodynamic (MHD) Stokes eigenvalue problem et al. based on the velocity-pressure formulation we study the residual type a posteriori error estimates of the mixed DGFEM using P k − P k−1 (k ≥ 1) element on shape-regular simplex meshes. For the Stokes equations the DGFEM was researched by [14][15][16][17][18][19] which laid a foundation for us to study further the Stokes eigenvalue problem. Among them, Badia et al. [14] proved the well-posedness of discrete DG formulation and studied the a priori error estimate of DGFEM with P k − P k−1 element without using the discrete inf-sup condition. Our main work is as follows: (1) We present the a posteriori error estimators for approximate eigenpairs. Referring to [20,21], using the enriching operator (see [22,23]) and the lifting operator (see [24,25]) we prove the reliability and efficiency of the estimator for eigenfunctions. We establish an identity (see Lemma 3.8) and use it to analyze the reliability of the estimator for eigenvalues. The characteristic of the adaptive DGFEM discussed in this paper is that for the Stokes eigenvalue problem in two and three-dimensional domains due to the usage of high-order elements it can capture smooth solutions and achieve the optimal convergence order for local less smooth solutions (eigenfunctions that have local singularity or local low smoothness) on adaptive locally refined graded meshes.
(2) We implement adaptive computing and the numerical results confirm our theoretical predictions and show that our method is stable, efficient and can obtain high-accuracy approximate eigenvalues. In the existing literature on the classical Stokes eigenvalue problem, Gedicke et al. [1,2] presented the approximate eigenvalues of 11, 10 and 9 significant digits in the unit square, the L-shaped domain and the slit domain, respectively, which are the most accurate approximations in the existing literature. In this paper, we obtain approximate eigenvalues that have the same accuracy as those in [1,2].
Note that C in different positions in this article represents different positive constants which is independent of mesh size h. We use a ≲ b to represent a ≤ Cb and a ⋍ b to represent a ≲ b and b ≲ a.

Preliminaries
Consider the following class of Stokes eigenvalue problems: in Ω, in Ω, u = 0, on ∂Ω, where Ω ⊂ R d (d = 2, 3) is a bounded polyhedral domain, u = (u 1 , ..., u d ) T is the velocity of the flow, p is the pressure, µ > 0 is the kinematic viscosity parameter of the fluid, λ is the eigenvalue of the problem (2.1) and A is a d × d symmetric semi-definite matrix whose elements belong to L ∞ (Ω). The problem (2.1) includes the classical Stokes eigenvalue problem in R d (d = 2, 3) and the MHD Stokes eigenvalue problem et al. When A is a zero matrix (2.1) is the classical Stokes eigenvalue problem. In the case of d = 2 when the problem (2.1) is the MHD Stokes eigenvalue problem (see [7,8]) where H 0 is the intensity of the externally applied magnetic field on the vertical direction, i.e., magnetic field H = (0, H 0 , 0), and Ha is the Hartmann number (see [7,8]) and when the problem (2.1) is the MHD Stokes eigenvalue problem while the magnetic field is applied horizontally.
For the sake of narrative simplicity we take µ = 1 in this paper. Let H ρ (Ω) be the Sobolev space on Ω of order ρ ≥ 0 equipped with the norm ∥ · ∥ ρ,Ω (denoted by where A(u, z) = (∇u, ∇z) + (Au, z), Let T h = {τ} be a regular simplex partition of Ω with the mesh diameter h = max τ∈T h h τ where h τ is the diameter of element τ. We use E i h and E b h to denote the set of interior faces (edges) and the set of faces (edges) on ∂Ω, respectively.
We use h F to denote the measure of F ∈ E h . We denote by (·, ·) τ and (·, ·) F the inner product in L 2 (τ) and L 2 (F), respectively. We use ω(τ) to represent the union of all elements which share at least one edge (face) with τ and use ω(F) to represent the union of the elements having in common with F.
The broken Sobolev space is defined by For any F ∈ E i h , there are two simplices τ + and τ − such that F = τ + ∩ τ − (e.g. see [14]). Let n + be the unit normal of F pointing from τ + to τ − and let n − = −n + .
The discrete velocity and pressure spaces are defined as follows (see [7]): where P k (τ) is the space of polynomials of degree less than or equal to k ≥ 1 on τ.
The DGFEM for the problem Here, γ/h F is the interior penalty parameter. We select γ = C 1 k 2 with C 1 = 10 in the actual numerical implementations from Remark 2.1 in [26].
Define the DG-norm as follows: Note that ∥ · ∥ h is equivalent to ||| · ||| on X h . From [27] we know that the continuity and coercivity properties hold: We consider the boundary problem : Given g ∈ (L 2 (Ω)) d , in Ω, u g = 0, on ∂Ω.
(2.10) From the Lax-Milgram theorem we have the existence and uniqueness of the velocity u in the space Z = {z ∈ X : b(z, ϱ) = 0, ∀ϱ ∈ W}. From the well-known inf-sup condition (see [28]): the stability of the pressure holds. The weak formulation of (2.10) reads: Find (u g , p g ) ∈ X × W such that and its DGFEM form reads: (2.14) We assume that the following regularity is valid: For any g ∈ (L 2 (Ω) where C is a positive constant independent of g. From Lemma 6.5 in [27] we can obtain the consistency of the DGFEM that is to say when (u g , p g ) is the solution of the boundary problem (2.10), there hold the following equations: Badia et al. [14], Cockburn et al. [29], Hansbo et al. [30] and Schötzau et al. [25] proved that (2.13) and (2.14) are well defined and gave the a priori error estimate. From [14] we obtain the following lemma.
We introduce the following auxiliary problem before estimating the error of velocity in the sense of L 2 norm: (2.25) Using the continuity of A h (·, ·) and the approximation property of I h ω we obtain Using the approximation property of ϑ h ϱ we can get Substituting E 1 , E 2 , E 3 into (2.25) and using (2.23) we obtain (2.24). The proof is completed. From the inf-sup condition and [14] we know that (2.11) and (2.12) are uniquely solvable and stable. Then, we define and it is valid that From [14] we also know that (2.13) and (2.14) are uniquely solvable and stable and we define Hence, Thus, (2.2), (2.3) and (2.4), (2.5) have the following equivalent operator forms, respectively: (2.29) Next, we will derive the error estimates for the eigenvalue problem.

From (2.24), (2.20) and (2.15) we have
Thus, we can obtain the following Lemma 2.3 (see Lemma 2.3 in [31]) from the Babuška-Osborn spectral approximation theory [32,33]. From (2.8) and (2.9) we know that ||| · ||| is a norm stronger than ∥ · ∥ h , i.e., ∥z∥ h ≲ |||z|||. Additionally, we have To show the validity of (2.31), using the trace theorem on the reference element and the scaling argument we have for any τ ∈ T h that and from the inverse inequality and the interpolation estimate and taking By the above inequality and (2.9) we obtain (2.31).

The a posteriori error estimates
Let (λ h , u h , p h ) ∈ R + × X h × W h be an approximate eigenpair. First, for each element τ ∈ T h we introduce the residuals: 3) identity matrix. Next, we introduce the following estimator η J τ to measure the jump of the approximate solution u h : The local error indictor is defined as Then, the global a posteriori error estimator is defined as We denote θ τ =int{ τ i ∩τ ∅τ i , τ i ∈ T h } for τ ∈ T h and use θ F to represent the set of all elements which share at least one node with face F. We denote by z I the Scott-Zhang interpolation function [34], then z I ∈ X ∩ X h and The lifting operator L : and has the following property (see [24,25]): Using the lifting operator, we define the following form: Lemma 3.1. Let (u g , p g ) and (u g h , p g h ) be the solutions of (2.11), (2.12) and (2.13), (2.14), respectively. Then, .

From (3.2) we deduce
For the third term, by the properties of the interpolation function z I we know [[z − z I ]] = 0. Therefore, from the definition of lifting operation L we have By the Cauchy-Schwartz inequality, (3.4) and (3.1) we obtain For the last term of (3.19) we get Substituting B 1 -B 5 into (3.19) results in (3.17).
In [22,23], the authors constructed the enriching operator E h : X h → X h ∩ X by averaging and proved the following lemma. Lemma 3.3. It is valid the following estimate: Theorem 3.2. Suppose that the conditions of Theorem 2.1 hold. Then, For any vector-valued polynomial function σ on F it is valid that Furthermore, for each b F σ there exists an extension σ b ∈ H 1 0 (ω(F)) satisfying σ b | F = b F σ and Using the standard arguments (see, e.g., Lemma 3.13 in [36]) and Lemmas 7 and 8 in [2], we can deduce the following local bounds.

31)
Proof. For any τ ∈ T h define the function R and K locally by From (3.22) and using λu + ∆u − Au − ∇p = 0, we have Using integration by parts and K| ∂τ = 0, we obtain Applying the Cauchy-Schwartz inequality yields
For any F ∈ E i h (Ω), [[u]] = 0 and for any F ∈ E h ∩ ∂Ω, u ⊗ n = 0. Therefore, we obtain (3.31) and finish the proof. Theorem 3.3. Suppose that the conditions of Theorem 2.1 hold. Then, the a posteriori error estimator η h is efficient: Proof. By using (2.16) and (2.17) we get We complete the proof.
Proof. Theorem 2.1 indicates that ∥u − u h ∥ 0 is a term of higher order than |||u − u h ||| + ∥p − p h ∥ 0 . Hence, from (3.36) and (3.21), we obtain  [17], for the problem (2.1) all analysis and conclusions in this paper are valid for the mixed DGFEM using the Q k − Q k−1 element.

Numerical experiments
When λ is a multiple eigenvalue the exact eigenfunction approximated by the discrete eigenfunction will change with the change of mesh diameter. In order to implement the adaptive algorithm better, we will conduct our numerical experiments on simple eigenpairs (multiplicity 1).
Based on [40][41][42], we design an adaptive DGFEM algorithm (ADGFEM) by adopting the standard adaptive loop with the steps solve, estimate, mark and refine with the a posteriori error estimator given in Section 3. We compile our program with the help of the iFEM package [43] and solve the matrix eigenvalue problem by means of the command ′ eigs ′ in MATLAB.
We adapt the following symbols in our tables: l: the lth iteration. λ k,h l : the kth approximate eigenvalue at the lth iteration. do f : the degrees of freedom at the lth iteration.
The error curves for the first eigenvalue of the classical Stokes eigenvalue problem are shown in  Error of the smallest eigenvalue error indicator error on adaptive mesh error on uniform mesh The line with slope -3 Figure 1. The error curves of the first eigenvalue by the ADGFEM using P 2 −P 1 element (left) and P 3 − P 2 element (right) for the classical Stokes eigenvalue problem in Ω slit .  Error of the smallest eigenvalue error indicator error on adaptive mesh error on uniform mesh The line with slope -3 Figure 3. The error curves of the first eigenvalue by the ADGFEM using P 2 −P 1 element (left) and P 3 − P 2 element (right) for the classical Stokes eigenvalue problem in Ω square .     We observe from Figures 1-6 that the error curves and error estimators curves for ADGFEM are both almost parallel to the straight line with a slope of −k which indicates that the error estimators are reliable and efficient and the adaptive algorithm can achieve the optimal convergence order. This is consistent with our theoretical results. We also observe from the error curves that under the same do f the approximations obtained by the ADGFEM are more accurate than those computed on uniform meshes.
The approximations of the first eigenvalue obtained by P 3 − P 2 element in Ω slit , Ω L and Ω square are listed in Tables 1-3, respectively. These eigenvalues have the same accuracy as those in [1,2] which achieve 9, 10 and 11 significant digits in Ω slit , Ω L and Ω square , respectively. Furthermore, it shows that our method is effective. The approximations of the fourth eigenvalue obtained by P 3 − P 2 element in Ω slit , Ω L and Ω square are listed in Tables 4-6, respectively.   The experiment is conducted in two three-dimensional domains: Ω 1 = (0, 1) 3 \ {0 ≤ x ≤ 0.5, 0 ≤ y ≤ 0.5, 0.5 ≤ z ≤ 1} and Ω 2 = (0, 1) 3 . In computation we select the initial mesh π h 0 with h 0 = Error of the smallest eigenvalue error indicator error on adaptive mesh The line with slope -2 Figure 9. Adaptive mesh after l=5 refinement times (left) and the error curves (right) of the first eigenvalue by the ADGFEM using P 3 − P 2 element for the classical Stokes eigenvalue problem in Ω 2 .
The reference values for the first eigenvalue of the classical Stokes eigenvalue problem are λ Ω 1 = 70.98560 and λ Ω 2 = 62.17341 for the domains Ω 1 and Ω 2 respectively, which are calculated by adaptive procedure with as much degrees of freedom as possible.
The adaptive refined meshes and the error curves are shown in Figures 8 and 9. We observe from Figures 8 and 9 that the error estimators are reliable and efficient and the adaptive algorithm achieve the optimal convergence order.

The numerical results for the MHD Stokes eigenvalue problem
We conduct experiments in Ω L and Ω square . We select θ = 0.5 and the initial mesh π h 0 with h 0 = √ 2 16 for the above two two-dimensional domains.
For the MHD Stokes eigenvalue problem with Ha = 5, we choose the values λ 1,square = 64.68920947 and λ 1,L = 40.2764915 as the reference values for Ω square and Ω L , respectively, and while Ha = 30, we choose the values λ 1,square = 234.34458093 and λ 1,L = 125.24247135 as the reference values for Ω square and Ω L respectively. These reference values are obtained by the ADGFEM using P 3 − P 2 element with as much degrees of freedom as possible.
The error curves for the first eigenvalue are shown in Figures 10 and 11 and the adaptive refined meshes for the first eigenvalue of the MHD Stokes eigenvalue problem by the ADGFEM are shown in Figure 12. We observe from Figures 10 and 11 that the error curves and error estimators curves are both approximately parallel to the line with slope −k, which indicates that the error estimators are reliable and efficient and the adaptive algorithm can achieve the optimal convergence order.
The approximations of the first eigenvalue for the MHD Stokes eigenvalue problem in Ω square using P 3 − P 2 element are listed in Tables 7 and 8, from which we can see that the approximate eigenvalues also has high accuracy.
We also use the ADGFEM with P k − P k−1 (k = 1, 2) element to calculate the classical Stokes eigenvalue problem and the MHD Stokes eigenvalue problem. The numerical results indicate that the discrete formulations are stable and effective. Due to article length limitations, these results are not listed in the paper. Error of the smallest eigenvalue error indicator error on adaptive mesh error on uniform mesh The line with slope -3 Figure 10. The error curves of the first eigenvalue by the ADGFEM using P 3 − P 2 element for the MHD Stokes eigenvalue problem with Ha = 5 in Ω square (left) and Ω L (right). Error of the smallest eigenvalue error indicator error on adaptive mesh error on uniform mesh The line with slope -3 Figure 11. The error curves of the first eigenvalue by the ADGFEM using P 3 − P 2 element for the MHD Stokes eigenvalue problem with Ha = 30 in Ω square (left) and Ω L (right). Figure 12. The adaptive meshes for the first eigenvalue of the MHD Stokes eigenvalue problem when Ha = 5 by the ADGFEM at l = 25 refinement times using P 3 − P 2 element in Ω L (left) and at l = 8 refinement times using P 3 − P 2 element in Ω square (right).

Conclusions
In this paper, for a class of Stokes eigenvalue problems including the classical Stokes eigenvalue problem in R d (d = 2, 3) and the MHD Stokes eigenvalue problem et al, based on the velocity-pressure formulation we studied the residual type a posteriori error estimates of the mixed DGFEM using P k − P k−1 (k ≥ 1) element on shape-regular simplex meshes. We proposed the a posteriori error estimator for approximate eigenpairs and proved the reliability and efficiency of the estimator for eigenfunctions and also analyzed their reliability for eigenvalues. The characteristic of the adaptive DGFEM is that it can use high-order elements and capture local low smooth solutions and can achieve the optimal convergence order O(do f − 2k d ) in two and three-dimensional domains. Our method is easy to implement on existing software packages. The numerical results confirmed our theoretical predictions and showed that our method is stable, efficient and can obtain high-accuracy approximate eigenvalues.

Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.