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

1.   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 velocity-pressure 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 $\mathbb{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 $\mathbb{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 $\mathbb{P}_{k}-\mathbb{P}_{k-1}$ ($k\geq 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 $\mathbb{P}_{k}-\mathbb{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\lesssim b$ to represent $a\leq C b$ and $a\backsimeq b$ to represent $a\lesssim b$ and $b\lesssim a$.

2.   Preliminaries

Consider the following class of Stokes eigenvalue problems:

where $\Omega\subset \mathbb{R}^{d}\; (d = 2, 3)$ is a bounded polyhedral domain, $\mathbf{u} = (u_{1}, ..., u_{d})^{T}$ is the velocity of the flow, $p$ is the pressure, $\mu > 0$ is the kinematic viscosity parameter of the fluid, $\lambda$ is the eigenvalue of the problem (2.1) and $A$ is a $d\times d$ symmetric semi-definite matrix whose elements belong to $L^{\infty}(\Omega)$.

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 $\mathbf{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 $\mu = 1$ in this paper.

Let $H^{\rho}(\Omega)$ be the Sobolev space on $\Omega$ of order $\rho\geq 0$ equipped with the norm $\|\cdot\|_{\rho, \Omega}$ (denoted by $\|\cdot\|_{\rho}$ for simplicity). $H_{0}^{1}(\Omega) = \{z\in H^{1}(\Omega), z|_{\partial\Omega} = 0\}$. We denote $\|\mathbf{z}\|_{\rho} = \sum\limits_{i = 1}^{d}\|z_{i}\|_{\rho}$ for $\mathbf{z} = (z_{1}, \cdots, z_{d})\in H^{\rho}(\Omega)^{d}$. We denote by $(\cdot, \cdot)$ the inner product in $L^{2}(\Omega)^{d}$ which is given by $(u, z) = \int_{\Omega}uz dx$ ($d = 1$) and $(\mathbf{u}, \mathbf{z}) = \int_{\Omega}\mathbf{u}\cdot\mathbf{z} dx$ ($d = 2, 3$). Define $\mathbb{X} = H_{0}^{1}(\Omega)^{d}$ with the norm $\|\mathbf{z}\|_\mathbb{X} = (\nabla\mathbf{z}, \nabla\mathbf{z})^{\frac{1}{2}}$ and define $\mathbb{W} = L_{0}^{2}(\Omega) = \{\varrho\in L^{2}(\Omega):(\varrho, 1) = 0\}$.

The weak formulation of (2.1) is given by: Find $(\lambda, \mathbf{u}, p)\in \mathbb{R}\times \mathbb{X} \times \mathbb{W}$, $\|\mathbf{u}\|_{0} = 1$ such that

where

Let $\mathcal{T}_{h} = \{\tau\}$ be a regular simplex partition of $\Omega$ with the mesh diameter $h = \max\limits_{\tau\in \mathcal{T}_{h}} h_{\tau}$ where $h_{\tau}$ is the diameter of element $\tau$. We use $\mathcal{E}_{h}^{i}$ and $\mathcal{E}_{h}^{b}$ to denote the set of interior faces (edges) and the set of faces (edges) on $\partial\Omega$, respectively. $\mathcal{E}_{h} = \mathcal{E}_{h}^{i}\cup\mathcal{E}_{h}^{b}$. We use $h_{F}$ to denote the measure of $F\in \mathcal{E}_{h}$. We denote by $(\cdot, \cdot)_{\tau}$ and $(\cdot, \cdot)_{F}$ the inner product in $L^{2}(\tau)$ and $L^{2}(F)$, respectively. We use $\omega(\tau)$ to represent the union of all elements which share at least one edge (face) with $\tau$ and use $\omega(F)$ to represent the union of the elements having in common with $F$.

The broken Sobolev space is defined by

For any $F\in\mathcal{E}_{h}^{i}$, there are two simplices $\tau^{+}$ and $\tau^{-}$ such that $F = \tau^{+}\cap\tau^{-}$ (e.g. see ^[14]). Let $\mathbf{n}^{+}$ be the unit normal of $F$ pointing from $\tau^{+}$ to $\tau^{-}$ and let $\mathbf{n}^{-} = -\mathbf{n}^{+}$.

For any $\varphi\in H^{1}(\Omega, \mathcal{T}_{h})$ we denote its jump and mean on $F\in\mathcal{E}_{h}^{i}$ by $[\![\varphi]\!] = \varphi^{+}\mathbf{n}^{+}+\varphi^{-}\mathbf{n}^{-}$ and $\mathbf{\{} \pmb{\mathsf{ φ}}\mathbf{\}} = \frac{1}{2}(\varphi^{+}+\varphi^{-})$, respectively, where $\varphi^{\pm} = \varphi|_{\tau^{\pm}}$. For $\boldsymbol{\varphi}\in H^{1}(\Omega, \mathcal{T}_{h})^{d}$ we denote by $[\![\boldsymbol{\varphi}]\!] = \boldsymbol{\varphi}^{+}\cdot\mathbf{n}^{+}+\boldsymbol{\varphi}^{-}\cdot\mathbf{n}^{-}$ the jump and $\{\boldsymbol{\varphi}\} = \frac{1}{2}(\boldsymbol{\varphi}^{+}+\boldsymbol{\varphi}^{-})$ the mean of $\boldsymbol{\varphi}$ on $F\in\mathcal{E}_{h}^{i}$. We also denote by $[\![\underline{\boldsymbol{\varphi}}]\!] = \boldsymbol{\varphi}^{+}\otimes \mathbf{n}^{+}+\boldsymbol{\varphi}^{-}\otimes \mathbf{n}^{-}$ the full jump of $\boldsymbol{\varphi}$ on $F\in\mathcal{E}_{h}^{i}$, where $\boldsymbol{\varphi}\otimes \mathbf{n} = [\varphi_{i}n_{j}]_{1\leq i, j\leq d}$, $\boldsymbol{\varphi} = (\varphi_{i}), \mathbf{n} = (n_{j})$. For tensors $\boldsymbol{\chi}\in H^{1}(\Omega, \mathcal{T}_{h})^{d\times d}$ we denote by $[\![\boldsymbol{\chi}]\!] = \boldsymbol{\chi}^{+}\mathbf{n}^{+}+\boldsymbol{\chi}^{-}\mathbf{n}^{-}$ the jump and $\{\boldsymbol{\chi}\} = \frac{1}{2}(\boldsymbol{\chi}^{+}+\boldsymbol{\chi}^{-})$ the mean on $F\in\mathcal{E}_{h}^{i}$.

For the sake of simplicity, when $F\in\mathcal{E}_{h}^{b}$ by modifying the above definitions appropriately, we obtain the jump and the mean on $\partial\Omega$. That is to say, we modify $\varphi^{-} = 0$ (similarly, $\boldsymbol{\varphi}^{-} = 0$ and $\boldsymbol{\chi}^{-} = 0$) to obtain the definition of jump on $\partial\Omega$ and modify $\varphi^{-} = \varphi^{+}$(similarly, $\boldsymbol{\varphi}^{-} = \boldsymbol{\varphi}^{+}$ and $\boldsymbol{\chi}^{-} = \boldsymbol{\chi}^{+}$) to obtain the definition of mean on $\partial\Omega$.

The discrete velocity and pressure spaces are defined as follows (see ^[7]):

where $\mathbb{P}_{k}(\tau)$ is the space of polynomials of degree less than or equal to $k\geq 1$ on $\tau$.

The DGFEM for the problem (2.1) is to find $(\lambda_{h}, \mathbf{u}_{h}, p_{h})\in \mathbb{R}\times\mathbb{X}_{h}\times \mathbb{W}_{h}$, $\|\mathbf{u}_{h}\|_{0} = 1$ such that

where

Here, $\gamma/h_{F}$ is the interior penalty parameter. We select $\gamma = 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 $\|\cdot\|_{h}$ is equivalent to $|||\cdot|||$ on $\mathbb{X}_{h}$.

From ^[27] we know that the continuity and coercivity properties hold:

We consider the boundary problem : Given $\boldsymbol{g}\in (L^{2}(\Omega))^{d}$,

From the Lax-Milgram theorem we have the existence and uniqueness of the velocity $\mathbf{u}$ in the space $\mathbf{Z} = \{\mathbf{z\in\mathbb{X}}: b(\mathbf{z}, \varrho) = 0, \forall \varrho\in \mathbb{W}\}$. From the well-known inf-sup condition (see ^[28]):

the stability of the pressure holds.

The weak formulation of (2.10) reads: Find $(\mathbf{u}^{g}, p^{g})\in \mathbb{X}\times \mathbb{W}$ such that

and its DGFEM form reads: Find $(\mathbf{u}^{g}_{h}, p^{g}_{h})\in \mathbb{X}_{h}\times \mathbb{W}_{h}$ such that

We assume that the following regularity is valid: For any $\boldsymbol{g}\in (L^{2}(\Omega))^{d}\; (d = 2, 3)$ there exists $(\mathbf{u}^{g}, p^{g})\in (H^{1+r}(\Omega)^{d}\times H^{r}(\Omega))\cap (W^{2, p}(\Omega)^{d}\times W^{1, p}(\Omega))\; (\frac{1}{2} < r\leq 1, \; p > \frac{2d}{d+1})$ satisfying (2.10) and

where $C$ is a positive constant independent of $\boldsymbol{g}$.

From Lemma 6.5 in ^[27] we can obtain the consistency of the DGFEM that is to say when $(\mathbf{u}^{g}, p^{g})$ is the solution of the boundary problem (2.10), there hold the following equations:

From (2.13), (2.14), (2.16) and (2.17) we have

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.

Lemma 2.1. Assume that $(\mathbf{u}^{g}, p^{g})\in H^{1+s}(\Omega)^{d}\times H^{s}(\Omega)$ ($r\leq s\leq k$) with $\boldsymbol{g}\in H^{l}(\Omega)^{d}$ ($0\leq l\leq k+1$). Then,

Proof. Since the bilinear form $\mathbb{A}(\cdot, \cdot)$ is coercive on $\mathbb{X}$, $\mathbb{A}_{h}(\cdot, \cdot)$ is also coercive on $\mathbb{X}_{h}$. Using the proof method of Theorem 4.1 in ^[14] we can obtain (2.20).

Let $I_{h}: \mathbb{X}\cap C^{0}(\overline{\Omega})^{d}\to \mathbb{X}_{h}\cap\mathbb{X}$ be the conforming element interpolation operator and let $\vartheta_{h}:\; H^{s}(\Omega)\rightarrow \mathbb{W}_{h}$ be the local $L^{2}$ projection operator satisfying $\vartheta_{h}p|_{\tau}\in \mathbb{P}_{k-1}(\tau)$ and

We introduce the following auxiliary problem before estimating the error of velocity in the sense of $L^{2}$ norm:

Using (2.15), we have

From Theorem 6.12 in ^[27] using the Nitsche's technique can prove the following lemma.

Lemma 2.2. Suppose that the conditions of Lemma 2.1 and (2.15) hold. Then,

Proof. From (2.18) and (2.19) we can derive

Using the continuity of $\mathbb{A}_{h}(\cdot, \cdot)$ and the approximation property of $I_{h}\boldsymbol{\omega}$ we obtain

Since $div \boldsymbol{\omega} = 0$, $\boldsymbol{\omega}\in[H_{0}^{1}(\Omega)]^{d}$, $[\boldsymbol{\omega}] = 0$ and $[I_{h}\boldsymbol{\omega}] = 0, \forall F\in \mathcal{E}_{h}$, we have

Using the approximation property of $\vartheta_{h}\varrho$ 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:

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 $|||\cdot|||$ is a norm stronger than $\|\cdot\|_{h}$, i.e., $\|z\|_{h}\lesssim|||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 $\tau\in \mathcal{T}_{h}$ that

and from the inverse inequality and the interpolation estimate and taking $w = \nabla (\mathbf{u}-I_{h}\mathbf{u})$ in (2.32) we deduce

By the above inequality and (2.9) we obtain (2.31).

Theorem 2.1. Let $(\lambda, \mathbf{u}, p)$ and $(\lambda_{h}, \mathbf{u}_{h}, p_{h})$ be the $j$th eigenpair of (2.2), (2.3) and (2.4), (2.5), respectively. Assume that the regularity estimate (2.15) is valid and $(\mathbf{u}, p)\in H^{1+s}(\Omega)^{d}\times H^{s}(\Omega)$ for some $s\in[r, k]$. Then,

Proof. In (2.11)–(2.14), we take $\boldsymbol{g} = \lambda\mathbf{u}$ then we obtain $\mathbf{u}^{g} = \lambda \mathbb{T}\mathbf{u}$, $\mathbf{u}^{g}_{h} = \lambda \mathbb{T}_{h}\mathbf{u}$, $p^{g} = \lambda \mathbb{S}\mathbf{u}$ and $p^{g}_{h} = \lambda \mathbb{S}_{h}\mathbf{u}$. Hence, using (2.20) we deduce

By using (2.16), (2.18), (2.19) and (2.36) we obtain

From Lemma 2.3 in ^[31] we know that

Substituting (2.37) and (2.24) into (2.39) we obtain (2.34).

Applying the triangle inequality and (2.27) we get

From (2.24), (2.38) and (2.39) we deduce

Then, from (2.40)–(2.42) we obtain

Thus, we get (2.33).

Combining (2.43) with (2.36) we get (2.35).

3.   The a posteriori error estimates

Let $(\lambda_{h}, \mathbf{u}_{h}, p_{h})\in \mathbb{R}^{+}\times\mathbb{X}_{h}\times \mathbb{W}_{h}$ be an approximate eigenpair. First, for each element $\tau\in\mathcal{T}_{h}$ we introduce the residuals:

where $\mathbf{I}$ denotes the $d\times d$ $(d = 2, 3)$ identity matrix. Next, we introduce the following estimator $\eta_{J_{\tau}}$ to measure the jump of the approximate solution $\mathbf{u}_{h}$:

The local error indictor is defined as

Then, the global a posteriori error estimator is defined as

We denote $\theta_{\tau} =$int$\{\bigcup\limits_{\overline{\tau}_{i}\cap \overline{\tau}\not = \emptyset}\bar{\tau}_{i}, \tau_{i}\in\mathcal{T}_{h}\}$ for $\tau\in\mathcal{T}_{h}$ and use $\theta_{F}$ to represent the set of all elements which share at least one node with face $F$. We denote by $\mathbf{z}^{I}$ the Scott-Zhang interpolation function ^[34], then $\mathbf{z}^{I}\in \mathbb{X}\cap\mathbb{X}_{h}$ and

Denote

The lifting operator $\mathcal{L}:\mathbb{X}(h)\rightarrow \underline{\sum}_{h}$ is defined by

and has the following property (see ^[24,25]):

Using the lifting operator, we define the following form:

by

Note that $\widetilde{\mathbb{A}}_{h} = \mathbb{A}_{h}$ on $\mathbb{X}_{h}\times \mathbb{X}_{h}$ and $\mathbb{A} = \widetilde{\mathbb{A}}_{h}$ on $\mathbb{X}\times \mathbb{X}$. The DGFEM presented in (2.4) and (2.5) is equivalent to finding $(\lambda_{h}, \mathbf{u}_{h}, p_{h})\in \mathbb{R}^{+}\times \mathbb{X}_{h}\times \mathbb{W}_{h}$ and satisfying

Lemma 3.1. Let $(\mathbf{u}^{g}, p^{g})$ and $(\mathbf{u}^{g}_{h}, p^{g}_{h})$ be the solutions of (2.11), (2.12) and (2.13), (2.14), respectively. Then,

Proof. For $\forall \mathbf{\bar{u}}\in \mathbb{X}$, from (2.11) we have

For $\forall \mathbf{\bar{u}}\in \mathbb{X}$, $\overline{p}\in \mathbb{W}$ we have

Combining the above two equations and taking $\mathbf{z} = \mathbf{u}^{g}-\mathbf{\bar{u}}$ we obtain

From the well-known inf-sup condition we obtain

Dividing both sides of (3.9) by $\|\mathbf{z}\|_{h}$ and taking supremum for $\mathbf{z}\in \mathbb{X}$ we get

Using the triangle inequality we obtain

Since $(\mathbf{\bar{u}}, \bar{p})$ is arbitrary and $\inf\limits_{\bar{p}\in \mathbb{W}}\|\bar{p}- p^{g}_{h}\|_{0} = 0$, the part $\lesssim$ in (3.8) is valid. The other part $\gtrsim$ in (3.8) is obvious.

Lemma 3.1 can be extended to the eigenvalue problem.

Theorem 3.1. Let $(\lambda, \mathbf{u}, p)$ and $(\lambda_{h}, \mathbf{u}_{h}, p_{h})$ be the $j$th eigenpair of (2.2), (2.3) and (2.4), (2.5), respectively. Then,

Proof. Using (2.26) and (2.27) we can obtain

In (2.11)–(2.14) we take $\boldsymbol{g} = \lambda_{h}\mathbf{u}_{h}$ and obtain $\mathbf{u}^{g} = \lambda_{h} \mathbb{T}\mathbf{u}_{h }$, $\mathbf{u}^{g}_{h} = \lambda_{h}\mathbb{T}_{h}\mathbf{u}_{h}$, $p^{g} = \lambda_{h}\mathbb{S}\mathbf{u}_{h }$ and $p^{g}_{h} = \lambda_{h}\mathbb{S}_{h}\mathbf{u}_{h}$. Therefore, from (3.8) we have

From (2.11) with $\boldsymbol{g} = \lambda_{h}\mathbf{u}_{h}$, (2.26) and (2.27) we deduce

Substituting (3.15) into (3.14) gives us

Theorem 2.1 indicates that $\|\lambda_{h}\mathbf{u}_{h}-\lambda\mathbf{u}\|_{0}$ is a small quantity of higher order compared with $\|\mathbf{u}-\mathbf{u}_{h}\|_{h}+\|p-p_{h}\|_{0}$. From (3.16) and (3.14) the side $\lesssim$ in (3.12) is true. The other side $\gtrsim$ in (3.12) is obvious.

Lemma 3.2. Under the conditions of Theorem 2.1,

Proof. Using (2.7), (3.6), (3.7) and the Green's formula we deduce that

By $\mathbf{z}^{I}\in \mathbb{X}\cap\mathbb{X}_{h}$ and (2.2)–(2.5) we obtain

Using (3.6), the Cauchy-Schwartz inequality, (3.1) and (3.2), (3.18) can be written as follows:

Next, we will analyze each item on the right-hand side of (3.19). Using the Cauchy-Schwartz inequality and the approximation property (3.1) and (3.2), we have

From (3.2) we deduce

For the third term, by the properties of the interpolation function $\mathbf{z}^{I}$ we know $[\![\mathbf{z}-\mathbf{z}^{I}]\!] = 0$. Therefore, from the definition of lifting operation $\mathcal{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}:\mathbb{X}_{h}\to \mathbb{X}_{h}\cap\mathbb{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,

Proof. Substituting (3.17) and (3.20) into (3.12), we obtain (3.21).

Let $b_{\tau}$ and $b_{F}$ be the standard bubble function on element $\tau$ and face $F$ ($d = 3$) or edge $F$ ($d = 2$) of $\tau$, respectively. Then, from ^[20,21,35] we obtain the following lemma.

Lemma 3.4. For any vector-valued polynomial function $\mathbf{z}_{h}$ on $\tau$,

For any vector-valued polynomial function $\sigma$ on $F$ it is valid that

Furthermore, for each $b_{F}\sigma$ there exists an extension $\sigma_{b}\in H_{0}^{1}(\omega(F))$ satisfying $\sigma_{b}|_{F} = b_{F}\sigma$ 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.

Lemma 3.5. Under the conditions of Theorem 2.1,

Proof. For any $\tau\in \mathcal{T}_{h}$ define the function $R$ and $K$ locally by

From (3.22) and using $\lambda\mathbf{u}+\Delta \mathbf{u}-A\mathbf{u}-\nabla p = 0$, we have

Using integration by parts and $K|_{\partial\tau} = 0$, we obtain

Applying the Cauchy-Schwartz inequality yields

From (3.23) and (3.24) we get

Dividing (3.32) by $h_{\tau}\|R\|_{0, \tau}$ and noting $\|\nabla\cdot \mathbf{u}_{h}\|_{0} = \|\nabla\cdot (\mathbf{u}_{h}-\mathbf{u})\|_{0}$, we obtain (3.29).

For any interior edge $F\in\mathcal{E}_{h}^{i}$ let the functions $R$ and $\Theta$ be such that

Using (3.26) and $[\![p\mathbf{I}-\nabla\mathbf{u}]\!]|_{F} = 0$ we get

Applying the Green's formula over each element of $\omega(F)$ we derive

Using $\lambda\mathbf{u}+\Delta\mathbf{u}-A\mathbf{u}-\nabla p = 0$ we deduce

Using the Cauchy-Schwartz inequality, (3.27) and (3.28) yieids

Combing the above estimates of $T_{1}$, $T_{2}$, $T_{3}$ and $T_{4}$, dividing (3.33) by $h_{F}^{1/2}\|R\|_{0, F}$ and summing over all interior edges of $\tau$ gives us (3.30).

For any $F\in \mathcal{E}_{h}^{i}(\Omega)$, $[\![\underline{\mathbf{u}}]\!] = 0$ and for any $F\in \mathcal{E}_{h}\cap \partial\Omega$, $\mathbf{u}\otimes\mathbf{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 $\eta_{h}$ is efficient:

Lemma 3.6. Let $(\lambda, \mathbf{u}, p)$ and $(\lambda_{h}, \mathbf{u}_{h}, p_{h})$ be the eigenpairs of (2.2), (2.3) and (2.4), (2.5), respectively. Then,

Proof. By using (2.16) and (2.17) we get

From (2.2) and (2.3) with $(\mathbf{z}, \varrho) = (\mathbf{u}, p)$, (2.4) and (2.5) with $(\mathbf{z}_{h}, \varrho_{h}) = (\mathbf{u_{h}}, p_{h})$ and (3.37), (3.38) we deduce

We complete the proof.

Theorem 3.4. Under the conditions of Theorem 2.1,

Proof. Theorem 2.1 indicates that $\|\mathbf{u}-\mathbf{u}_{h}\|_{0}$ is a term of higher order than $|||\mathbf{u}-\mathbf{u}_{h}|||+\|p-p_{h}\|_{0}$. Hence, from (3.36) and (3.21), we obtain

Thus, from (2.39) and (3.21) we obtain (3.39).

Remark 3.1. Theorem 2.1 indicates that $\|\lambda_{h}\mathbf{u}_{h}-\lambda \mathbf{u}\|_{0}$ is a small quantity of higher order than $\|\mathbf{u}-\mathbf{u}_{h}\|_{h}+\|p-p_{h}\|_{0}$. Theorems 3.2 and 3.3 show that the estimator $\eta_{h}$ for the eigenfunction error $\|\mathbf{u}-\mathbf{u}_{h}\|_{h}+\|p-p_{h}\|_{0}$ is reliable and efficient up to data oscillation. Therefore, a good graded mesh is generated by the adaptive algorithm for the estimator, which makes the eigenfunction error $\|\mathbf{u}-\mathbf{u}_{h}\|_{h}+\|p-p_{h}\|_{0}$ reach the optimal convergence rate $O(dof^{-\frac{2k}{d}})$. Hence, from ^[37,38,39] we can look forward to getting $\sum\limits_{\tau\in\mathcal{T}_{h}}h_{\tau}^{2r}(|\mathbf{u}-I_{h}\mathbf{u}|_{1+r, \tau}^{2} +\|p-\vartheta_{h}p\|_{r}^{2})\lesssim dof^{-\frac{2k}{d}}$. Thereby from (3.39) we have $|\lambda-\lambda_{h}|\lesssim dof^{-\frac{2k}{d}}$. Therefore, we think $\eta_{h}^{2}$ can be regarded as the error estimator of $\lambda_{h}$.

Remark 3.2. Based on ^[17], for the problem (2.1) all analysis and conclusions in this paper are valid for the mixed DGFEM using the $\mathbb{Q}_{k}-\mathbb{Q}_{k-1}$ element.

4.   Numerical experiments

When $\lambda$ 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 $l$th iteration.

$\lambda_{k, h_{l}}$: the $k$th approximate eigenvalue at the $l$th iteration.

$dof$: the degrees of freedom at the $l$th iteration.

4.1.   The numerical results for the classical Stokes eigenvalue problem

4.1.1.   The results for two-dimensional domains

The experiment is conducted in three two-dimensional domains the slit domain $\Omega_{slit} =$ $(-1, 1)^{2}\setminus \{{0\leq x\leq 1, y = 0}\}$, the L-shaped domain $\Omega_{L} =$$(-1, 1)^{2}\setminus [0, 1]\times [-1, 0]$ and the unit square domain $\Omega_{square} =$ $(0, 1)^{2}$. In the step mark we select the parameter $\theta = 0.5$, and the initial mesh $\pi_{h_{0}}$ with $h_{0} = \frac{\sqrt{2}}{16}$ for the above three two-dimensional domains.

The reference values for the first eigenvalue of the classical Stokes eigenvalue problem are $\lambda_{1, slit} = 29.9168629$, $\lambda_{1, L} = 32.13269465$ and $\lambda_{1, square} = 52.344691168$ for $\Omega_{slit}$, $\Omega_{L}$ and $\Omega_{square}$, respectively (see ^[1,2]) and the reference value for the fourth eigenvalue reads $\lambda_{4, square} = 128.209584313$ in $\Omega_{square}$ (see ^[2]). We choose the values $\lambda_{4, slit} = 40.1527333966$ and $\lambda_{4, L} = 48.9835839778$ as the reference values for the $\Omega_{slit}$ and $\Omega_{L}$ respectively, which are obtained by adaptive procedure using $\mathbb{P}_{3}-\mathbb{P}_{2}$ element with as much degrees of freedom as possible.

The error curves for the first eigenvalue of the classical Stokes eigenvalue problem are shown in Figures 1–3 and the fourth eigenvalue are shown in Figures 4–6. The adaptive refined meshes for the first eigenvalue of the classical Stokes eigenvalue problem by the ADGFEM are shown in Figure 7.

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 $dof$ the approximations obtained by the ADGFEM are more accurate than those computed on uniform meshes.

The approximations of the first eigenvalue obtained by $\mathbb{P}_{3}-\mathbb{P}_{2}$ element in $\Omega_{slit}$, $\Omega_{L}$ and $\Omega_{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 $\Omega_{slit}$, $\Omega_{L}$ and $\Omega_{square}$, respectively. Furthermore, it shows that our method is effective. The approximations of the fourth eigenvalue obtained by $\mathbb{P}_{3}-\mathbb{P}_{2}$ element in $\Omega_{slit}$, $\Omega_{L}$ and $\Omega_{square}$ are listed in Tables 4–6, respectively.

4.1.2.   The results for three-dimensional domains

The experiment is conducted in two three-dimensional domains: $\Omega_{1} = (0, 1)^{3}\setminus \{0\leq x\leq 0.5, 0\leq y\leq 0.5, 0.5\leq z\leq 1\}$ and $\Omega_{2} = (0, 1)^{3}$. In computation we select the initial mesh $\pi_{h_{0}}$ with $h_{0} = \frac{\sqrt{3}}{8}$ and $\theta = 0.25$.

The reference values for the first eigenvalue of the classical Stokes eigenvalue problem are $\lambda_{\Omega_{1}} = 70.98560$ and $\lambda_{\Omega_{2}} = 62.17341$ for the domains $\Omega_{1}$ and $\Omega_{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.

4.2.   The numerical results for the MHD Stokes eigenvalue problem

We conduct experiments in $\Omega_{L}$ and $\Omega_{square}$. We select $\theta = 0.5$ and the initial mesh $\pi_{h_{0}}$ with $h_{0} = \frac{\sqrt{2}}{16}$ for the above two two-dimensional domains.

For the MHD Stokes eigenvalue problem with $Ha = 5$, we choose the values $\lambda_{1, square} = 64.68920947$ and $\lambda_{1, L} = 40.2764915$ as the reference values for $\Omega_{square}$ and $\Omega_{L}$, respectively, and while $Ha = 30$, we choose the values $\lambda_{1, square} = 234.34458093$ and $\lambda_{1, L} = 125.24247135$ as the reference values for $\Omega_{square}$ and $\Omega_{L}$ respectively. These reference values are obtained by the ADGFEM using $\mathbb{P}_{3}-\mathbb{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 $\Omega_{square}$ using $\mathbb{P}_{3}-\mathbb{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 $\mathbb{P}_{k}-\mathbb{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.

5.   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 $\mathbb{P}_{k}-\mathbb{P}_{k-1} \; (k\geq1)$ 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(dof^{-\frac{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.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 11561014, 11761022), and Science and Technology Planning Project of Guizhou Province (Guizhou Kehe fundamental research-ZK[2022] No.324).

Conflict of interest

The authors declare that this work does not have any conflicts of interest.