Solving the Stokes Equations via a Least Squares Weak Galerkin Method
July 12, 2026
Abstract
We present a least-squares weak Galerkin (LS-WG) finite element method for solving the Stokes equations on arbitrary polygonal and polyhedral meshes. By utilizing discrete weak derivatives on discontinuous polynomial spaces, the proposed framework naturally accommodates complex domain geometries and general partitions. Crucially, this least-squares formulation bypasses the traditional inf-sup (LBB) compatibility condition, transforming the standard indefinite saddle-point problem into an inherently symmetric and positive definite (SPD) discrete linear system. We establish the well-posedness of the numerical scheme and rigorously derive optimal-order error estimates in a custom discrete energy norm. Specifically, we prove convergence rates of \(\mathcal{O}(h^k)\) for the discrete projection error and \(\mathcal{O}(h^{k-1})\) for the global approximation error when employing polynomials of degree \(k \ge 1\) for the velocity field and \(k-1\) for the pressure. Extensive numerical experiments confirm these theoretical convergence rates, demonstrating the method’s robustness, geometric flexibility, and overall efficiency.
1 Introduction↩︎
The Stokes equations are fundamental to fluid mechanics, governing the behavior of incompressible viscous fluids in the creeping flow regime where the Reynolds number is very low. In this regime, viscous forces dominate inertial effects, making the model essential for simulating diverse phenomena across science and engineering. Key applications range from microfluidics and lab-on-a-chip technologies to biological systems—such as the swimming mechanisms of microorganisms—and large-scale geophysical processes, including mantle convection and glacier dynamics. Furthermore, the Stokes problem provides the underlying linear structure for the full, nonlinear Navier-Stokes equations. Consequently, the development of robust and geometrically flexible numerical solvers for the Stokes system remains a critical area of computational research.
In this paper, we introduce a least squares weak Galerkin (WG) finite element method for the Stokes equations. We seek a velocity field \({\boldsymbol{u}}\) and a pressure \(p\) that satisfy the following system: \[\begin{align} -\Delta{\boldsymbol{u}}+ \nabla p &= {\boldsymbol{f}}\quad \text{in } \Omega, \tag{1} \\ \nabla \cdot {\boldsymbol{u}}&= 0 \quad \text{in } \Omega, \tag{2} \\ {\boldsymbol{u}}&= {\boldsymbol{g}}\quad \text{on } \partial\Omega, \tag{3} \end{align}\] where \(\Omega \subset \mathbb{R}^d\) (\(d=2,3\)) is a polygonal or polyhedral domain, \({\boldsymbol{f}}\) is the body force, and \({\boldsymbol{g}}\) represents the Dirichlet boundary data.
The numerical approximation of the Stokes problem presents significant challenges, primarily due to the indefinite saddle-point nature of the system. Standard mixed finite element methods (FEM) require that the velocity and pressure spaces satisfy the restrictive discrete inf-sup (Ladyzhenskaya-Babuška-Brezzi or LBB) condition to ensure stability and well-posedness [1]. Constructing inf-sup stable element pairs on complex or non-matching grids is often computationally intensive and geometrically limiting.
To overcome these constraints, discontinuous Galerkin (DG) methods were developed to relax inter-element continuity, thereby enhancing geometric flexibility and enabling divergence-free velocity approximations [2]. More recently, the weak Galerkin (WG) finite element method has emerged as a powerful alternative. The WG framework [3]–[27] utilizes discrete weak derivatives and enforces continuity across interfaces through specially designed stabilizers. First applied to the Stokes problem in [28], the WG method naturally accommodates arbitrary polygonal and polyhedral meshes. Recent refinements have further introduced lifting operators and advanced stabilization to ensure pressure robustness across complex polytopal domains [29]–[31].
A notable evolution within this field is the Primal-Dual Weak Galerkin (PDWG) method [32]–[46], which formulates numerical approximations as constrained minimization problems. By enforcing the governing equations as weak constraints through Lagrange multipliers, PDWG offers favorable stability and symmetry for various problems, including non-self-adjoint linear transport equations [32], [37], [43], [47].
Building upon these developments, this work proposes a novel least-squares weak Galerkin finite element method for the Stokes problem 1 –3 . Unlike traditional WG formulations that may yield non-symmetric systems, or primal-dual WG (PDWG) methods that increase the total number of unknowns via dual variables, the proposed LS-WG method minimizes a least-squares functional in a weak sense. This approach offers several distinct advantages:
Symmetry and Positivity: The formulation transforms the standard indefinite saddle-point problem into a discrete linear system that is inherently symmetric and positive definite (SPD). This allows for the use of highly efficient iterative solvers, such as the conjugate gradient (CG) method.
Mesh Flexibility: The method inherits the standard WG capability of handling arbitrary polytopal meshes with pendant nodes, eliminating the need for matching grids or specialized continuous basis functions.
No Inf-Sup Requirement: The least-squares framework bypasses the need to satisfy the traditional LBB compatibility condition, greatly simplifying the choice of velocity and pressure finite element spaces.
This paper provides a rigorous theoretical and numerical foundation for the LS-WG scheme. We establish the existence and uniqueness of the discrete solution and derive optimal-order error estimates in a carefully defined discrete energy norm. Specifically, we prove convergence rates of \(\mathcal{O}(h^k)\) for the discrete projection error and \(\mathcal{O}(h^{k-1})\) for the global approximation error when employing polynomials of degree \(k \ge 1\) for the velocity field and \(k-1\) for the pressure. Finally, we present extensive numerical experiments to validate these theoretical convergence rates and demonstrate the robustness and geometric flexibility of the algorithm.
The remainder of this paper is organized as follows. In Section 2, we define the fundamental discrete weak differential operators: the discrete weak Laplacian, the discrete weak gradient, and the discrete weak divergence. Section 3 presents the development of the least-squares weak Galerkin finite element scheme for the Stokes equations. The error equations for the LS-WG approximation are derived in Section 4, followed by a rigorous derivation of optimal-order error estimates, specifically establishing the \(\mathcal{O}(h^k)\) and \(\mathcal{O}(h^{k-1})\) convergence rates in the discrete energy norm, in Section 5. Finally, Section 6 provides a series of numerical experiments to validate the theoretical findings and demonstrate the robustness of the method.
Let \(D \subset \mathbb{R}^d\) (\(d=2, 3\)) be an open, bounded domain with a Lipschitz continuous boundary. We employ standard notation for Sobolev spaces. For any \(s \ge 0\), \(H^s(D)\) denotes the Sobolev space with the associated inner product \((\cdot,\cdot)_{s,D}\), norm \(\|\cdot\|_{s,D}\), and seminorm \(|\cdot|_{s,D}\). The space \(H^0(D)\) coincides with \(L^2(D)\), for which the inner product and norm are simplified to \((\cdot,\cdot)_{D}\) and \(\|\cdot \|_{D}\), respectively. In cases where \(D=\Omega\), the subscript \(D\) will be omitted for brevity.
For vector-valued functions, we define the space \(H(\text{div}; D)\) as the set of functions in \([L^2(D)]^d\) whose divergence also belongs to \(L^2(D)\): \[H(\text{div}, D) = \left\{ {\boldsymbol{v}}\in [L^2(D)]^d : \nabla \cdot {\boldsymbol{v}}\in L^2(D) \right\}.\]
2 Discrete Weak Differential Operators↩︎
The hallmark of the weak Galerkin method is the replacement of classical (strong) derivatives with weak derivatives. These weak operators are defined on a class of weak functions that are discontinuous across element boundaries. For the least-squares WG formulation of the Stokes problem 1 –3 , we require three specific weak differential operators: the weak Laplacian, the weak divergence, and the weak gradient. For completeness, this section reviews the formal definitions of these weak differential operators and their discrete counterparts [28], [30], [48].
2.1 Discrete Weak Laplacian↩︎
Let \(T\) be a polygonal or polyhedral domain in \(\mathbb{R}^d\) (\(d=2,3\)) with boundary \(\partial T\). We define a weak function for the velocity field on \(T\) as an ordered triplet \({\boldsymbol{v}}=\{{\boldsymbol{v}}_0, {\boldsymbol{v}}_b, \mathbf{v}_g\}\), where the components are defined as follows: \({\boldsymbol{v}}_0 \in [L^2(T)]^d\) represents the value of \({\boldsymbol{v}}\) in the interior of \(T\); \({\boldsymbol{v}}_b \in [L^2(\partial T)]^d\) represents the value of \({\boldsymbol{v}}\) on the boundary \(\partial T\); \(\mathbf{v}_g \in [L^2(\partial T)]^{d\times d}\) represents the value of the gradient \(\nabla {\boldsymbol{v}}\) on the boundary \(\partial T\). Crucially, \({\boldsymbol{v}}_b\) and \(\mathbf{v}_g\) are defined independently and are not required to be the traces of \({\boldsymbol{v}}_0\) or \(\nabla {\boldsymbol{v}}_0\) on \(\partial T\). We denote the space of such weak functions on \(T\) by: \[\begin{align} \mathcal{W}(T) = \bigg\{ {\boldsymbol{v}}= \{{\boldsymbol{v}}_0, {\boldsymbol{v}}_b, \mathbf{v}_g\} & : {\boldsymbol{v}}_0 \in [L^2(T)]^d,\\ &\quad {\boldsymbol{v}}_b \in [L^2(\partial T)]^d, \;\mathbf{v}_g \in [L^2(\partial T)]^{d\times d} \bigg\}. \end{align}\]
The weak Laplacian of a function \({\boldsymbol{v}}\in \mathcal{W}(T)\), denoted by \(\Delta_w {\boldsymbol{v}}\), is a linear functional in the dual space of \([H^2(T)]^d\) defined by: \[\begin{align} (\Delta_w {\boldsymbol{v}}, {\boldsymbol{w}})_T &=({\boldsymbol{v}}_0, \Delta{\boldsymbol{w}})_T - \langle {\boldsymbol{v}}_b, \nabla {\boldsymbol{w}}\cdot \mathbf{n}\rangle_{\partial T} \\ & \quad \; + \langle \mathbf{v}_g \cdot \mathbf{n}, {\boldsymbol{w}}\rangle_{\partial T}, \quad \forall {\boldsymbol{w}}\in [H^2(T)]^d, \end{align}\] where \(\mathbf{n}\) is the outward unit normal vector on \(\partial T\).
Let \(P_r(T)\) denote the space of polynomials of degree at most \(r\) on \(T\). The discrete weak Laplacian, \(\Delta_{w,r,T} {\boldsymbol{v}}\), is the unique polynomial in \([P_r(T)]^d\) satisfying: \[\begin{align} \label{dislap}\begin{aligned} (\Delta_{w,r,T} {\boldsymbol{v}}, {\boldsymbol{w}})_T & = ({\boldsymbol{v}}_0, \Delta {\boldsymbol{w}})_T - \langle {\boldsymbol{v}}_b, \nabla {\boldsymbol{w}}\cdot \mathbf{n} \rangle_{\partial T} \\ &\quad \; + \langle \mathbf{v}_g \cdot \mathbf{n}, {\boldsymbol{w}}\rangle_{\partial T}, \quad \forall {\boldsymbol{w}}\in [P_r(T)]^d. \end{aligned} \end{align}\tag{4}\] If \({\boldsymbol{v}}_0\) possesses sufficient regularity (e.g., \({\boldsymbol{v}}_0 \in [H^2(T)]^d\)), the discrete weak Laplacian can be equivalently expressed through integration by parts: \[\begin{align} \label{dislap2} \begin{aligned} (\Delta_{w,r,T} {\boldsymbol{v}}, {\boldsymbol{w}})_T & = (\Delta {\boldsymbol{v}}_0, {\boldsymbol{w}})_T + \langle {\boldsymbol{v}}_0 - {\boldsymbol{v}}_b, \nabla {\boldsymbol{w}}\cdot \mathbf{n} \rangle_{\partial T} \\ &\quad \;+ \langle (\mathbf{v}_g - \nabla {\boldsymbol{v}}_0) \cdot \mathbf{n}, {\boldsymbol{w}}\rangle_{\partial T}, \quad \forall {\boldsymbol{w}}\in [P_r(T)]^d. \end{aligned} \end{align}\tag{5}\]
2.2 Discrete Weak Divergence↩︎
For any \({\boldsymbol{v}}\in \mathcal{W}(T)\), we define the weak divergence \(\nabla_w \cdot {\boldsymbol{v}}\) as a linear functional in the dual space of \(H^1(T)\) such that for each \(\varphi \in H^1(T)\): \[\label{wd} (\nabla_w \cdot {\boldsymbol{v}}, \varphi)_T = -({\boldsymbol{v}}_0, \nabla\varphi)_T + \langle {\boldsymbol{v}}_b \cdot \mathbf{n}, \varphi\rangle_{\partial T}.\tag{6}\] The discrete weak divergence operator \(\nabla_{w,r,T} \cdot {\boldsymbol{v}}\in P_r(T)\) is the unique polynomial satisfying: \[\label{d-d} (\nabla_{w,r,T} \cdot {\boldsymbol{v}}, \varphi)_T = -({\boldsymbol{v}}_0, \nabla\varphi)_T + \langle {\boldsymbol{v}}_b \cdot \mathbf{n}, \varphi\rangle_{\partial T}, \quad \forall \varphi \in P_r(T).\tag{7}\] By integration by parts, for \({\boldsymbol{v}}_0 \in H(\text{div}, T)\), we have the equivalent form: \[\label{div2} (\nabla_{w,r,T} \cdot {\boldsymbol{v}}, \varphi)_T = (\nabla \cdot {\boldsymbol{v}}_0, \varphi)_T + \langle ({\boldsymbol{v}}_b - {\boldsymbol{v}}_0) \cdot \mathbf{n}, \varphi\rangle_{\partial T}, \quad \forall \varphi \in P_r(T).\tag{8}\]
2.3 Discrete Weak Gradient↩︎
For the pressure variable, we define a weak scalar function space \(\mathcal{V}(T)\) as: \[\mathcal{V}(T) = \{q = \{q_0, q_b\} : q_0 \in L^2(T), \;q_b \in L^2(\partial T)\}.\] The weak gradient of \(q \in \mathcal{V}(T)\), denoted \(\nabla_w q\), is a linear functional in the dual space of \([H^1(T)]^d\) such that for each \({\boldsymbol{w}}\in [H^1(T)]^d\): \[\label{wg} (\nabla_w q, {\boldsymbol{w}})_T = -(q_0, \nabla \cdot {\boldsymbol{w}})_T + \langle q_b, {\boldsymbol{w}}\cdot \mathbf{n}\rangle_{\partial T}.\tag{9}\] The discrete weak gradient \(\nabla_{w,r, T} q \in [P_r(T)]^d\) is defined as the unique polynomial vector satisfying: \[\label{d-g} (\nabla_{w,r, T} q, {\boldsymbol{w}})_T = -(q_0, \nabla \cdot {\boldsymbol{w}})_T + \langle q_b, {\boldsymbol{w}}\cdot \mathbf{n}\rangle_{\partial T}, \quad \forall {\boldsymbol{w}}\in [P_r(T)]^d.\tag{10}\] For \(q_0 \in H^1(T)\), this is equivalently given by: \[\label{nabla2} (\nabla_{w,r, T} q, {\boldsymbol{w}})_T = (\nabla q_0, {\boldsymbol{w}})_T + \langle q_b - q_0, {\boldsymbol{w}}\cdot \mathbf{n}\rangle_{\partial T}, \quad \forall {\boldsymbol{w}}\in [P_r(T)]^d.\tag{11}\]
3 A Least Squares Weak Galerkin Finite Element Scheme↩︎
Let \(\mathcal{T}_h\) be a partition of the domain \(\Omega\) with mesh size \(h\), consisting of shape regular polygons or polyhedra. Denote by \(\mathcal{E}_h\) the set of all edges (in 2D) or flat faces (in 3D) in \(\mathcal{T}_h\), and let \(\mathcal{E}_h^0 = \mathcal{E}_h \setminus \partial\Omega\) be the set of all interior edges or faces.
For any integer \(k \ge 2\), we define the weak Galerkin finite element space for the velocity variable as follows: \[\begin{align} \label{W-h} \begin{aligned} \mathcal{W}_h &= \{ {\boldsymbol{v}}=\{{\boldsymbol{v}}_0, {\boldsymbol{v}}_b, \mathbf{v}_g\} : {\boldsymbol{v}}_0 \in [P_{k}(T)]^d, \;{\boldsymbol{v}}_b \in [P_{k}(e)]^d, \\ &\qquad \; \mathbf{v}_g \in [P_{k-1}(e)]^{d\times d}, \;\forall T \in \mathcal{T}_h, \;\forall e \subset \partial T \}. \end{aligned} \end{align}\tag{12}\] We emphasize that there is a single, uniquely defined value for \({\boldsymbol{v}}_b\) on each edge \(e \in \mathcal{E}_h^0\). For the pressure variable, the associated weak finite element space is given by: \[\begin{align} \label{V-h} \begin{aligned} \mathcal{V}_h & = \bigg\{ q=\{q_0, q_b\} : q_0 \in P_{k-1}(T) \cap L_0^2(\Omega), \\ &\qquad \; q_b \in P_{k-1}(e), \;\forall T \in \mathcal{T}_h, \;\forall e \subset \partial T \bigg\}. \end{aligned} \end{align}\tag{13}\] Let \(\mathcal{W}_h^0\) denote the subspace of \(\mathcal{W}_h\) consisting of discrete weak functions with vanishing boundary values; that is, \[\mathcal{W}_h^0 = \left\{ {\boldsymbol{v}}=\{{\boldsymbol{v}}_0, {\boldsymbol{v}}_b, \mathbf{v}_g\} \in \mathcal{W}_h : {\boldsymbol{v}}_b = \mathbf{0} \text{ on } \partial\Omega \right\}.\]
The discrete weak Laplacian \(\Delta_{w,k-2}\) and the discrete weak divergence \((\nabla_{w,k-1}\cdot)\) on the space \(\mathcal{W}_h\) are computed locally on each element \(T\) using 4 and 7 , respectively. More precisely, they are defined by: \[\begin{align} (\Delta_{w,k-2}{\boldsymbol{v}})|_T &= \Delta_{w,k-2, T} ({\boldsymbol{v}}|_T), \quad \forall {\boldsymbol{v}}\in \mathcal{W}_h,\\ (\nabla_{w,k-1}\cdot{\boldsymbol{v}})|_T &= \nabla_{w,k-1, T}\cdot ({\boldsymbol{v}}|_T), \quad \forall {\boldsymbol{v}}\in \mathcal{W}_h. \end{align}\] For notational simplicity, we shall hereafter drop the subscripts \(k-2\) and \(k-1\) and denote the discrete weak Laplacian and discrete weak divergence simply by \(\Delta_w\) and \(\nabla_w\cdot\), respectively.
Similarly, the discrete weak gradient \((\nabla_{w,k-2})\) on the space \(\mathcal{V}_h\) is computed using 10 on each element \(T\): \[(\nabla_{w,k-2} q)|_T = \nabla_{w,k-2, T} (q|_T), \quad \forall q \in \mathcal{V}_h.\] Again, we will drop the subscript \(k-2\) and denote the discrete weak gradient simply by \(\nabla_w\).
Let \(Q_{0}^k\), \(Q_{0}^{k-1}\), and \(Q_{0}^{k-2}\) denote the local \(L^2\) projection operators from \(L^2(T)\) onto \(P_k(T)\), \(P_{k-1}(T)\), and \(P_{k-2}(T)\), respectively. For each face \(e \in \mathcal{E}_h\), let \(Q_{b}^k\) and \(Q_{b}^{k-1}\) denote the \(L^2\) projections from \(L^2(e)\) onto \(P_{k}(e)\) and \(P_{k-1}(e)\), respectively.
We are now in a position to present the least squares weak Galerkin finite element scheme for the Stokes equations 1 –3 . First, we introduce the following three bilinear forms: for \({\boldsymbol{u}}, {\boldsymbol{v}}, {\boldsymbol{w}}\in {\mathcal{W}}_h\) and \(p, q\in {\mathcal{V}}_h\), \[\begin{align} s_1({\boldsymbol{v}}, {\boldsymbol{w}}) &= \sum_{T\in \mathcal{T}_h} h_T^{-1}\langle {\boldsymbol{v}}_0-{\boldsymbol{v}}_b, {\boldsymbol{w}}_0-{\boldsymbol{w}}_b\rangle_{\partial T} \\ &\qquad\qquad + h_T \langle \nabla {\boldsymbol{v}}_0-\mathbf{v}_g, \nabla{\boldsymbol{w}}_0-\mathbf{w}_g\rangle_{\partial T}, \\ s_2(p, q) &= \sum_{T\in \mathcal{T}_h} h_T\langle p_0-p_b, q_0-q_b\rangle_{\partial T}, \\ a(({\boldsymbol{u}},p),({\boldsymbol{v}}, q)) &= \sum_{T\in \mathcal{T}_h} (-\Delta_w{\boldsymbol{u}}+\nabla_w p, -\Delta_w{\boldsymbol{v}}+\nabla_w q)_T + (\nabla_w\cdot{\boldsymbol{u}}, \nabla_w\cdot{\boldsymbol{v}})_T. \end{align}\]
Figure 1: 
.
Theorem 1. The Least Squares Weak Galerkin method 1 possesses a unique solution.
Proof. It suffices to show that the solution to 1 is identically zero when the data vanishes (i.e., \({\boldsymbol{f}}=\mathbf{0}\) and \({\boldsymbol{g}}=\mathbf{0}\)). Assuming \({\boldsymbol{f}}=\mathbf{0}\) and \({\boldsymbol{g}}=\mathbf{0}\), we choose the test functions \({\boldsymbol{v}}={\boldsymbol{u}}_h\) and \(q=p_h\) in 1 to obtain: \[a(({\boldsymbol{u}}_h,p_h),({\boldsymbol{u}}_h, p_h)) + s_1({\boldsymbol{u}}_h,{\boldsymbol{u}}_h) + s_2(p_h,p_h) = 0.\] Due to the positivity of these bilinear forms, this implies that \(-\Delta_w{\boldsymbol{u}}_h+\nabla_w p_h = \mathbf{0}\) and \(\nabla_w \cdot {\boldsymbol{u}}_h = 0\) on each element \(T\in \mathcal{T}_h\). Furthermore, the stabilizers dictate that \({\boldsymbol{u}}_0={\boldsymbol{u}}_b\), \(\nabla{\boldsymbol{u}}_0=\mathbf{u}_g\), and \(p_0=p_b\) on each \(\partial T\).
Using 5 alongside the facts that \({\boldsymbol{u}}_0={\boldsymbol{u}}_b\) and \(\nabla{\boldsymbol{u}}_0=\mathbf{u}_g\) on each \(\partial T\), we have \[(\Delta_{w} {\boldsymbol{u}}_h, {\boldsymbol{w}})_T = (\Delta {\boldsymbol{u}}_0, {\boldsymbol{w}})_T, \quad \forall {\boldsymbol{w}}\in [P_{k-2}(T)]^d,\] which yields \(\Delta_{w} {\boldsymbol{u}}_h = \Delta {\boldsymbol{u}}_0\).
Similarly, using 8 and \({\boldsymbol{u}}_0={\boldsymbol{u}}_b\) on each \(\partial T\) gives \[(\nabla_w \cdot {\boldsymbol{u}}_h, \phi)_T = (\nabla \cdot {\boldsymbol{u}}_0, \phi)_T, \quad \forall \phi \in P_{k-1}(T),\] yielding \(\nabla_w \cdot {\boldsymbol{u}}_h = \nabla \cdot {\boldsymbol{u}}_0\).
Applying 11 and \(p_0=p_b\) on each \(\partial T\) yields \[(\nabla_w p_h, {\boldsymbol{w}})_T = (\nabla p_0, {\boldsymbol{w}})_T, \quad \forall {\boldsymbol{w}}\in [P_{k-2}(T)]^d,\] which implies \(\nabla_w p_h = \nabla p_0\).
Because \({\boldsymbol{u}}_0={\boldsymbol{u}}_b\) and \(\nabla{\boldsymbol{u}}_0=\mathbf{u}_g\) across all element boundaries, it follows that \({\boldsymbol{u}}_0\) and its gradient are globally continuous, implying \(\nabla{\boldsymbol{u}}_0 \in [H(\text{div},\Omega)]^{d \times d}\). Likewise, \(p_0=p_b\) implies that \(p_0\) is continuous over the entire domain \(\Omega\). Substituting \(\Delta_{w} {\boldsymbol{u}}_h = \Delta {\boldsymbol{u}}_0\), \(\nabla_w \cdot {\boldsymbol{u}}_h = \nabla \cdot {\boldsymbol{u}}_0\), and \(\nabla_w p_h = \nabla p_0\) into the identities \(-\Delta_w{\boldsymbol{u}}_h+\nabla_w p_h=\mathbf{0}\) and \(\nabla_w\cdot{\boldsymbol{u}}_h=0\) yields the strong classical system: \[\begin{align} -\Delta{\boldsymbol{u}}_0+\nabla p_0 &= \mathbf{0} \quad \text{in } \Omega, \tag{14}\\ \nabla\cdot{\boldsymbol{u}}_0 &= 0 \quad \text{in } \Omega, \tag{15}\\ {\boldsymbol{u}}_0 &= \mathbf{0} \quad \text{on } \partial\Omega. \tag{16} \end{align}\] Recall that \(p_0 \in L_0^2(\Omega)\). Therefore, the continuous system 14 –16 admits only the trivial solution \({\boldsymbol{u}}_0 \equiv \mathbf{0}\) and \(p_0 \equiv 0\) in \(\Omega\). This forces the boundary terms \({\boldsymbol{u}}_b \equiv \mathbf{0}\), \({\boldsymbol{u}}_g \equiv \mathbf{0}\) and \(p_b \equiv 0\), which ultimately guarantees that \({\boldsymbol{u}}_h \equiv \mathbf{0}\) and \(p_h \equiv 0\) throughout \(\Omega\).
This completes the proof. ◻
We define the energy norm on the finite element space \(\mathcal{W}_h^0 \times \mathcal{V}_h\) by \[{|||} ({\boldsymbol{u}}_h, p_h){|||}^2 = a(({\boldsymbol{u}}_h, p_h), ({\boldsymbol{u}}_h, p_h)) + s_1({\boldsymbol{u}}_h, {\boldsymbol{u}}_h) + s_2(p_h, p_h).\] Following the mathematical logic established in Theorem 1, it is straightforward to verify that \({|||} \cdot {|||}\) satisfies all the properties of a norm on \(\mathcal{W}_h^0 \times \mathcal{V}_h\), providing a rigorous framework for the subsequent optimal-order error analysis.
4 Error Equations↩︎
In this section, we derive the error equation that governs the relationship between the exact solutions \(({\boldsymbol{u}}, p)\) and their LS-WG approximations \(({\boldsymbol{u}}_h, p_h)\). Let \(({\boldsymbol{u}}, p)\) be the exact solutions to the Stokes problem 1 –3 , and let \({\boldsymbol{u}}_h \in \mathcal{W}_h\) and \(p_h\in \mathcal{V}_h\) be the solutions to the discrete problem 1. We define the projection operators \(Q_h^1 {\boldsymbol{u}}\) and \(Q_h^2 p\) into the weak spaces, and the corresponding error functions, as follows: \[\begin{align} e_{{\boldsymbol{u}}_h} &= {\boldsymbol{u}}_h - Q_h^1 {\boldsymbol{u}}= \{{\boldsymbol{u}}_0 - Q_0^k {\boldsymbol{u}}, \;{\boldsymbol{u}}_b - Q_b^k {\boldsymbol{u}}, \;\mathbf{u}_g - Q_b^{k-1}(\nabla {\boldsymbol{u}})\}, \\ e_{p_h} &= p_h - Q_h^2 p = \{p_0 - Q_0^{k-1}p, \;p_b - Q_b^{k-1}p\}. \end{align}\]
Lemma 1 (Commutative Properties). The projection operators \(Q_h^1\) and \(Q_h^2\) satisfy the following commutative properties: \[\begin{align} \Delta_w(Q_h^1 {\boldsymbol{w}}) &= Q_0^{k-2}(\Delta {\boldsymbol{w}}), \quad \forall {\boldsymbol{w}}\in [H^2(\Omega)]^d, \label{eq:commute1} \\ \nabla_w\cdot(Q_h^1 {\boldsymbol{w}}) &= Q_0^{k-1}(\nabla \cdot {\boldsymbol{w}}), \quad \forall {\boldsymbol{w}}\in [H^1(\Omega)]^d, \label{eq:commute2} \\ \nabla_w(Q_h^2 q) &= Q_0^{k-2}(\nabla q), \quad \forall q \in H^1(\Omega). \label{eq:commute3} \end{align}\] {#eq: sublabel=eq:eq:commute1,eq:eq:commute2,eq:eq:commute3}
Proof. For any test function \({\boldsymbol{v}}\in [P_{k-2}(T)]^d\), the definition of the discrete weak Laplacian 4 gives: \[\begin{align} (\Delta_w (Q_h^1 {\boldsymbol{w}}), {\boldsymbol{v}})_T &= (Q_0^k {\boldsymbol{w}}, \Delta {\boldsymbol{v}})_T - \langle Q_b^k {\boldsymbol{w}}, \nabla {\boldsymbol{v}}\cdot \mathbf{n} \rangle_{\partial T} + \langle Q_b^{k-1}(\nabla {\boldsymbol{w}}) \cdot \mathbf{n}, {\boldsymbol{v}}\rangle_{\partial T} \\ &= ({\boldsymbol{w}}, \Delta {\boldsymbol{v}})_T - \langle {\boldsymbol{w}}, \nabla {\boldsymbol{v}}\cdot \mathbf{n} \rangle_{\partial T} + \langle \nabla {\boldsymbol{w}}\cdot \mathbf{n}, {\boldsymbol{v}}\rangle_{\partial T} \\ &= (\Delta {\boldsymbol{w}}, {\boldsymbol{v}})_T = (Q_0^{k-2} (\Delta {\boldsymbol{w}}), {\boldsymbol{v}})_T. \end{align}\] This confirms ?? .
For any scalar function \(\varphi \in P_{k-1}(T)\), the definition of the discrete weak divergence 7 implies: \[\begin{align} (\nabla_w \cdot (Q_h^1 {\boldsymbol{w}}), \varphi)_T &= -(Q_0^{k} {\boldsymbol{w}}, \nabla \varphi)_T + \langle Q_b^k {\boldsymbol{w}}\cdot \mathbf{n}, \varphi \rangle_{\partial T} \\ &= -({\boldsymbol{w}}, \nabla \varphi)_T + \langle {\boldsymbol{w}}\cdot \mathbf{n}, \varphi \rangle_{\partial T} \\ &= (\nabla \cdot {\boldsymbol{w}}, \varphi)_T = (Q_0^{k-1} (\nabla \cdot {\boldsymbol{w}}), \varphi)_T. \end{align}\] This confirms ?? .
Finally, for any vector function \({\boldsymbol{v}}\in [P_{k-2}(T)]^d\), the definition of the discrete weak gradient 10 yields: \[\begin{align} (\nabla_w (Q_h^2 q), {\boldsymbol{v}})_T &= -(Q_0^{k-1} q, \nabla \cdot {\boldsymbol{v}})_T + \langle Q_b^{k-1} q, {\boldsymbol{v}}\cdot \mathbf{n} \rangle_{\partial T} \\ &= -(q, \nabla \cdot {\boldsymbol{v}})_T + \langle q, {\boldsymbol{v}}\cdot \mathbf{n} \rangle_{\partial T} \\ &= (\nabla q, {\boldsymbol{v}})_T = (Q_0^{k-2} (\nabla q), {\boldsymbol{v}})_T. \end{align}\] This confirms ?? and completes the proof. ◻
Lemma 2 (Error Equation). For any test functions \({\boldsymbol{v}}\in \mathcal{W}_h^0\) and \(q \in \mathcal{V}_h\), the error functions \(e_{{\boldsymbol{u}}_h}\) and \(e_{p_h}\) satisfy the following identity: \[\begin{align} \label{erroreqn}\begin{aligned} &\quad \; a( (e_{{\boldsymbol{u}}_h}, e_{p_h}), ({\boldsymbol{v}}, q)) + s_1(e_{{\boldsymbol{u}}_h}, {\boldsymbol{v}}) + s_2(e_{p_h}, q)\\ & = -s_1(Q_h^1 {\boldsymbol{u}}, {\boldsymbol{v}}) - s_2(Q_h^2 p, q). \end{aligned} \end{align}\qquad{(1)}\]
Proof. Testing the exact momentum equation 1 with \(-\Delta_w {\boldsymbol{v}}+ \nabla_w q\) and the continuity equation 2 with \(\nabla_w \cdot {\boldsymbol{v}}\) on each element \(T \in \mathcal{T}_h\), and summing over all elements, we obtain: \[\begin{align} \label{eq:err95proof1} \begin{aligned} &\quad \; \sum_{T \in \mathcal{T}_h} (- \Delta {\boldsymbol{u}}+ \nabla p, - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T + (\nabla\cdot {\boldsymbol{u}}, \nabla_w\cdot {\boldsymbol{v}})_T \\ & = \sum_{T \in \mathcal{T}_h} ({\boldsymbol{f}}, - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T. \end{aligned} \end{align}\tag{17}\] Note that by definition, \(-\Delta_w {\boldsymbol{v}}+ \nabla_w q \in [P_{k-2}(T)]^d\) and \(\nabla_w \cdot {\boldsymbol{v}}\in P_{k-1}(T)\) on each element \(T\). Therefore, we can apply the properties of the \(L^2\) projections alongside the commutative properties ?? –?? to deduce: \[\begin{align} &\quad \; (- \Delta {\boldsymbol{u}}, - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T \\ &= (-Q_0^{k-2}(\Delta {\boldsymbol{u}}), - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T = (-\Delta_w (Q_h^1 {\boldsymbol{u}}), - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T, \\ &\quad \;(\nabla p, - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T \\ &= (Q_0^{k-2}(\nabla p), - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T = (\nabla_w (Q_h^2 p), - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T, \\ &\quad \;(\nabla\cdot {\boldsymbol{u}}, \nabla_w\cdot {\boldsymbol{v}})_T \\ &= (Q_0^{k-1}(\nabla\cdot {\boldsymbol{u}}), \nabla_w\cdot {\boldsymbol{v}})_T = (\nabla_w\cdot (Q_h^1 {\boldsymbol{u}}), \nabla_w\cdot {\boldsymbol{v}})_T. \end{align}\] Substituting these identities into 17 recovers the bilinear form \(a(\cdot,\cdot)\): \[\label{eq:err95proof2} a((Q_h^1 {\boldsymbol{u}}, Q_h^2 p), ({\boldsymbol{v}}, q)) = \sum_{T \in \mathcal{T}_h} ({\boldsymbol{f}}, - \Delta_w {\boldsymbol{v}}+ \nabla_w q)_T.\tag{18}\] Subtracting 18 from the discrete LS-WG numerical scheme 1 yields: \[a(({\boldsymbol{u}}_h - Q_h^1 {\boldsymbol{u}}, p_h - Q_h^2 p), ({\boldsymbol{v}}, q)) + s_1({\boldsymbol{u}}_h, {\boldsymbol{v}}) + s_2(p_h, q) = 0.\] Using the linearity of the stabilization forms, we rewrite \(s_1({\boldsymbol{u}}_h, {\boldsymbol{v}}) = s_1(e_{{\boldsymbol{u}}_h} + Q_h^1 {\boldsymbol{u}}, {\boldsymbol{v}})\) and \(s_2(p_h, q) = s_2(e_{p_h} + Q_h^2 p, q)\). Substituting these into the above equation directly yields the desired error equation: \[a((e_{{\boldsymbol{u}}_h}, e_{p_h}), ({\boldsymbol{v}}, q)) + s_1(e_{{\boldsymbol{u}}_h}, {\boldsymbol{v}}) + s_2(e_{p_h}, q) = -s_1(Q_h^1 {\boldsymbol{u}}, {\boldsymbol{v}}) - s_2(Q_h^2 p, q).\] This completes the proof. ◻
5 Error Estimates↩︎
In this section, we establish optimal-order error estimates for the LS-WG approximation in the energy norm. Throughout the analysis, we denote by \(C\) a generic positive constant independent of the mesh parameter \(h\).
Lemma 3 (Approximation Properties). Let \(\mathcal{T}_h\) be a shape-regular finite element partition of \(\Omega\). For any \({\boldsymbol{u}}\in [H^{k+1}(\Omega)]^d\) and \(p\in H^{k}(\Omega)\), the following approximation estimates hold: \[\begin{align} \sum_{T \in \mathcal{T}_h} \|{\boldsymbol{u}}- Q_0^k {\boldsymbol{u}}\|_{s,T}^2 &\le C h^{2(k+1-s)} \|{\boldsymbol{u}}\|_{k+1}^2, \quad s=0, 1, 2, \label{error95uh} \\ \sum_{T \in \mathcal{T}_h} \|\Delta{\boldsymbol{u}}- Q_0^{k-2} (\Delta{\boldsymbol{u}})\|_{T}^2 &\le C h^{2(k-1)} \|{\boldsymbol{u}}\|_{k+1}^2, \label{error95uh2} \\ \sum_{T \in \mathcal{T}_h} \|p - Q_0^{k-1} p\|_{s,T}^2 &\le C h^{2(k-s)} \|p\|_{k}^2, \quad s=0, 1, \label{error95ph} \\ \sum_{T \in \mathcal{T}_h} \|\nabla p - Q_0^{k-2} (\nabla p)\|_{T}^2 &\le C h^{2(k-1)} \|p\|_{k}^2. \label{error95ph2} \end{align}\] {#eq: sublabel=eq:error95uh,eq:error95uh2,eq:error95ph,eq:error95ph2} Furthermore, for any \(\phi \in H^1(T)\), the following trace inequality is valid: \[\|\phi\|_{\partial T}^2 \le C \left( h_T^{-1} \|\phi\|_T^2 + h_T \|\nabla \phi\|_T^2 \right). \label{trace95H1}\qquad{(2)}\]
Theorem 2. Let \(({\boldsymbol{u}}, p) \in [H^{k+1}(\Omega)]^d \times H^k(\Omega)\) be the exact solutions of the Stokes problem 1 –3 , and let \(({\boldsymbol{u}}_h, p_h) \in \mathcal{W}_h \times \mathcal{V}_h\) be the LS-WG numerical solution defined by 1. Then, there exists a constant \(C > 0\) such that \[\label{main95estimate} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}) {|||} \le C h^{k} (\|{\boldsymbol{u}}\|_{k+1} + \|p\|_{k}).\qquad{(3)}\]
Proof. By setting the test functions \({\boldsymbol{v}}= e_{{\boldsymbol{u}}_h}\) and \(q = e_{p_h}\) in the error equation ?? , and recalling the definition of the energy norm \({|||} \cdot {|||}\), we have: \[{|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}) {|||}^2 = - s_1(Q_h^1 {\boldsymbol{u}}, e_{{\boldsymbol{u}}_h}) - s_2(Q_h^2 p, e_{p_h}).\]
Bound for \(s_1(Q_h^1 {\boldsymbol{u}}, e_{{\boldsymbol{u}}_h})\): Applying the Cauchy-Schwarz inequality, the trace inequality ?? , and the approximation estimate ?? , we obtain: \[\begin{align} \label{s1} \begin{aligned} &\quad \; s_1(Q_h^1 {\boldsymbol{u}}, e_{{\boldsymbol{u}}_h}) \\ &\leq \left( \sum_{T \in \mathcal{T}_h} h_T^{-1} \|Q_0^k {\boldsymbol{u}}- Q_b^k {\boldsymbol{u}}\|_{\partial T}^2 + h_T \|\nabla Q_0^{k}{\boldsymbol{u}}- Q_b^{k-1}(\nabla {\boldsymbol{u}})\|_{\partial T}^2 \right)^{\frac{1}{2}} \\ &\qquad {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||} \\ &\le C \left( \sum_{T \in \mathcal{T}_h} h_T^{-2} \|Q_0^k {\boldsymbol{u}}- {\boldsymbol{u}}\|_{T}^2 + \|Q_0^k {\boldsymbol{u}}- {\boldsymbol{u}}\|_{1, T}^2 \right. \\ &\quad \left. + \|\nabla Q_0^{k}{\boldsymbol{u}}- \nabla {\boldsymbol{u}}\|_{T}^2 + h_T^2\|\nabla Q_0^{k}{\boldsymbol{u}}- \nabla {\boldsymbol{u}}\|_{1, T}^2 \right)^{\frac{1}{2}} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||} \\ &\le C h^{k} \|{\boldsymbol{u}}\|_{k+1} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||}. \end{aligned} \end{align}\tag{19}\]
Bound for \(s_2(Q_h^2 p, e_{p_h})\): Similarly, applying the Cauchy-Schwarz inequality, the trace inequality ?? and the estimate ?? yields: \[\begin{align} \label{s2} s_2(Q_h^2 p, e_{p_h}) &\leq \left( \sum_{T \in \mathcal{T}_h} h_T \|Q_0^{k-1} p - Q_b^{k-1} p\|_{\partial T}^2 \right)^{\frac{1}{2}} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||} \nonumber \\ &\le C \left( \sum_{T \in \mathcal{T}_h} \|Q_0^{k-1} p - p\|_{T}^2 + h_T^2\|Q_0^{k-1} p - p\|_{1, T}^2 \right)^{\frac{1}{2}} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||} \nonumber \\ &\leq C h^{k} \|p\|_{k} {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||}. \end{align}\tag{20}\] Substituting 19 and 20 back into the energy norm identity and dividing by \({|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}){|||}\) completes the proof. ◻
Theorem 3. Let \(({\boldsymbol{u}}, p) \in [H^{k+1}(\Omega)]^d \times H^{k}(\Omega)\) be the exact solutions of the Stokes problem 1 –3 , and let \(({\boldsymbol{u}}_h, p_h) \in \mathcal{W}_h \times \mathcal{V}_h\) be the LS-WG numerical solution defined by 1. Then, there exists a constant \(C > 0\) such that \[\label{main95estimate2} {|||} ({\boldsymbol{u}}-{\boldsymbol{u}}_h, p - p_h) {|||} \le C h^{k-1} (\|{\boldsymbol{u}}\|_{k+1} + \|p\|_{k}).\qquad{(4)}\]
Proof. By the triangle inequality and the definition of the error functions, we have: \[\begin{align} &\quad \; {|||} ({\boldsymbol{u}}-{\boldsymbol{u}}_h, p- p_h) {|||}^2\\ &\leq {|||} ({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, p- Q_h^2p ) {|||}^2 + {|||} ( Q_h^1{\boldsymbol{u}}-{\boldsymbol{u}}_h, Q_h^2p-p_h ) {|||}^2 \\ &= {|||} ({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, p- Q_h^2p ) {|||}^2 + {|||} (e_{{\boldsymbol{u}}_h}, e_{p_h}) {|||}^2. \end{align}\] Expanding the first term using the definition of the energy norm, and applying the commutative properties ?? –?? , we obtain: \[\begin{align} &\quad \; {|||} ({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, p- Q_h^2p ) {|||}^2 \\ &= a(({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, p- Q_h^2p ), ({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, p- Q_h^2p )) \\ &\quad + s_1({\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}, {\boldsymbol{u}}-Q_h^1{\boldsymbol{u}}) + s_2 (p- Q_h^2p, p- Q_h^2p ) \\ &= \sum_{T\in \mathcal{T}_h} \|-\Delta {\boldsymbol{u}}+ Q_0^{k-2}(\Delta {\boldsymbol{u}}) + \nabla p - Q_0^{k-2}(\nabla p) \|_T^2 \\ &\quad + \sum_{T\in \mathcal{T}_h} \|\nabla \cdot {\boldsymbol{u}}- Q_0^{k-1}(\nabla\cdot {\boldsymbol{u}})\|_T^2 \\ &\quad + s_1(Q_h^1{\boldsymbol{u}}, Q_h^1{\boldsymbol{u}}) + s_2 (Q_h^2p, Q_h^2p). \end{align}\] Note that \(s_1({\boldsymbol{u}}, \cdot) = 0\) and \(s_2(p, \cdot) = 0\) since the continuous solutions have no jumps across element boundaries. Using the approximation estimates ?? and ?? , along with the bounds established for the stabilizers in 19 and 20 , and the discrete error bound ?? , we have: \[\begin{align} &\quad \; {|||} ({\boldsymbol{u}}-{\boldsymbol{u}}_h, p- p_h) {|||}^2 \\ &\leq C \sum_{T\in \mathcal{T}_h} \left( \|\Delta {\boldsymbol{u}}- Q_0^{k-2}(\Delta {\boldsymbol{u}})\|_T^2 + \|\nabla p - Q_0^{k-2}(\nabla p) \|_T^2 \right) \\ &\quad + C \sum_{T\in \mathcal{T}_h} \|\nabla \cdot {\boldsymbol{u}}- Q_0^{k-1}(\nabla\cdot {\boldsymbol{u}})\|_T^2 + C h^{2k} (\|{\boldsymbol{u}}\|^2_{k+1} + \|p\|^2_{k}) \\ &\leq C h^{2k-2} \|{\boldsymbol{u}}\|^2_{k+1} + C h^{2k-2} \|p\|^2_{k} + C h^{2k} (\|{\boldsymbol{u}}\|^2_{k+1} + \|p\|^2_{k}) \\ &\leq C h^{2k-2} (\|{\boldsymbol{u}}\|^2_{k+1} + \|p\|^2_{k}). \end{align}\] Taking the square root of both sides yields the final estimate ?? . ◻
6 Numerical experiment↩︎
In the first numerical test, we solve the Stokes equations 1 –3 , where \(\Omega=(0,1)\times(0,1)\), \[\begin{align} \label{sol1} \begin{aligned} \;{\boldsymbol{u}}&=\begin{pmatrix}8(1-2y)(y-y^2)(x-x^2)^2\\ -8 (1-2x)(x-x^2)(y-y^2)^2\end{pmatrix}, \;\\\ p&=4 (1-2x)(x-x^2)(1-2y)(y-y^2) . \end{aligned} \end{align}\tag{21}\] The computation is done on the triangular grids shown in Figure 2 and on the non-convex pentagonal grids shown in Figure 3, by the weak Galerkin \(P_k\)-\(P_{k-1}^2\) finite elements defined in 12 and 13 , \(k=2,3\) and \(4\). The results are listed in Tables 1-3, where we can see that the optimal orders of convergence are achieved roughly. In these tables, \(G_i\) denotes the \(i\)-th grid, shown in Figure 2 or Figure 3 .
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On triangular grids (Figure [f-g1]) | ||||||
| 1 | 0.603E-02 | — | 0.463E-01 | — | 0.236E-01 | — |
| 2 | 0.308E-02 | 1.0 | 0.170E-01 | 1.4 | 0.147E-01 | 0.7 |
| 3 | 0.814E-03 | 1.9 | 0.889E-01 | 0.0 | 0.910E-02 | 0.7 |
| 4 | 0.193E-03 | 2.1 | 0.486E-01 | 0.9 | 0.258E-02 | 1.8 |
| 5 | 0.490E-04 | 2.0 | 0.265E-01 | 0.9 | 0.653E-03 | 2.0 |
| 6 | 0.127E-04 | 2.0 | 0.137E-01 | 1.0 | 0.163E-03 | 2.0 |
| On pentagonal grids (Figure [f-g2]) | ||||||
| 1 | 0.256E-01 | — | 0.259E+00 | — | 0.598E-01 | — |
| 2 | 0.663E-02 | 1.9 | 0.576E-01 | 2.2 | 0.301E-01 | 1.0 |
| 3 | 0.115E-02 | 2.5 | 0.791E-01 | 0.0 | 0.108E-01 | 1.5 |
| 4 | 0.217E-03 | 2.4 | 0.401E-01 | 1.0 | 0.334E-02 | 1.7 |
| 5 | 0.529E-04 | 2.0 | 0.227E-01 | 0.8 | 0.884E-03 | 1.9 |
| 6 | 0.143E-04 | 1.9 | 0.123E-01 | 0.9 | 0.224E-03 | 2.0 |
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On triangular grids (Figure [f-g1]) | ||||||
| 1 | 0.472E-02 | — | 0.942E-02 | — | 0.297E-01 | — |
| 2 | 0.602E-03 | 3.0 | 0.142E+00 | 0.0 | 0.109E-01 | 1.4 |
| 3 | 0.775E-04 | 3.0 | 0.460E-01 | 1.6 | 0.177E-02 | 2.6 |
| 4 | 0.522E-05 | 3.9 | 0.115E-01 | 2.0 | 0.216E-03 | 3.0 |
| 5 | 0.296E-06 | 4.1 | 0.259E-02 | 2.1 | 0.256E-04 | 3.1 |
| On pentagonal grids (Figure [f-g2]) | ||||||
| 1 | 0.246E-01 | — | 0.866E+00 | — | 0.179E+00 | — |
| 2 | 0.229E-02 | 3.4 | 0.161E+00 | 2.4 | 0.210E-01 | 3.1 |
| 3 | 0.233E-03 | 3.3 | 0.399E-01 | 2.0 | 0.322E-02 | 2.7 |
| 4 | 0.162E-04 | 3.8 | 0.961E-02 | 2.1 | 0.409E-03 | 3.0 |
| 5 | 0.105E-05 | 4.0 | 0.222E-02 | 2.1 | 0.532E-04 | 2.9 |
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On triangular grids (Figure [f-g1]) | ||||||
| 1 | 0.312E-02 | — | 0.343E+00 | — | 0.327E-01 | — |
| 2 | 0.265E-03 | 3.6 | 0.145E+00 | 1.2 | 0.601E-02 | 2.4 |
| 3 | 0.116E-04 | 4.5 | 0.138E-01 | 3.4 | 0.320E-03 | 4.2 |
| 4 | 0.363E-06 | 5.0 | 0.142E-02 | 3.3 | 0.183E-04 | 4.1 |
| On pentagonal grids (Figure [f-g2]) | ||||||
| 1 | 0.298E-01 | — | 0.140E+01 | — | 0.214E+00 | — |
| 2 | 0.186E-02 | 4.0 | 0.183E+00 | 2.9 | 0.142E-01 | 3.9 |
| 3 | 0.741E-04 | 4.6 | 0.252E-01 | 2.9 | 0.110E-02 | 3.7 |
| 4 | 0.250E-05 | 4.9 | 0.320E-02 | 3.0 | 0.715E-04 | 3.9 |
For a 3D numerical test, we solve the Stokes equations 1 –3 , where \(\Omega=(0,1)^3\), \[\begin{align} \label{sol2} \begin{aligned} \;{\boldsymbol{u}}&=\begin{pmatrix}- 2^{10}(x-x^2)^2 (y-y^2)^2 (1-2z)(z-z^2)\\ \;2^{10}(x-x^2)^2 (y-y^2)^2 (1-2z)(z-z^2)\\ 2^{10} (2 x^3 - 2 y^3 - 3 x^2 + 3 y^2 + x - y) (z-z^2)^2\end{pmatrix}, \\ p&=z^3 -\frac{1}{4} . \end{aligned} \end{align}\tag{22}\] We also solve the problem with a smoother solution: \[\begin{align} \label{s3} \begin{aligned} \;{\boldsymbol{u}}&=\begin{pmatrix}e^y \\ e ^ z\\ e^x \end{pmatrix}, & p&=y -\frac{1}{2} . \end{aligned} \end{align}\tag{23}\]
The computation is done on the cuboid grids shown in Figure 4, on the wedge-shape grids shown in Figure 5, and on the tetrahedral grids shown in Figure 6, by the weak Galerkin \(P_k\)-\(P_{k-1}^2\) finite elements defined in 12 and 13 , \(k=2,3\) and \(4\). The results are listed in Tables 4-6, where we can see barely that the optimal orders of convergence are achieved, due to large condition numbers and the computer accuracy.
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On cuboid grids (Figure [f-g3]) for [s3] | ||||||
| 1 | 0.195E-01 | 0.00 | 0.147E+00 | 0.00 | 0.733E-01 | 0.00 |
| 2 | 0.284E-02 | 2.78 | 0.679E-01 | 1.11 | 0.242E-01 | 1.60 |
| 3 | 0.395E-03 | 2.85 | 0.171E-01 | 1.99 | 0.523E-02 | 2.21 |
| 4 | 0.521E-04 | 2.92 | 0.416E-02 | 2.04 | 0.121E-02 | 2.11 |
| 5 | 0.684E-05 | 2.93 | 0.101E-02 | 2.04 | 0.291E-03 | 2.06 |
| On wedge grids (Figure [f-g4]) for [s3] | ||||||
| 1 | 0.120E-01 | 0.00 | 0.237E+00 | 0.00 | 0.507E-01 | 0.00 |
| 2 | 0.191E-02 | 2.65 | 0.115E+00 | 1.05 | 0.228E-01 | 1.15 |
| 3 | 0.269E-03 | 2.83 | 0.285E-01 | 2.00 | 0.469E-02 | 2.28 |
| 4 | 0.365E-04 | 2.88 | 0.101E-01 | 1.50 | 0.104E-02 | 2.17 |
| On tetrahedral grids (Figure [f-g5]) for [sol2] | ||||||
| 1 | 0.116E+00 | 0.00 | 0.243E+01 | 0.00 | 0.726E+00 | 0.00 |
| 2 | 0.836E-01 | 0.47 | 0.255E+01 | 0.00 | 0.492E+00 | 0.56 |
| 3 | 0.239E-01 | 1.81 | 0.207E+01 | 0.30 | 0.179E+00 | 1.46 |
| 4 | 0.420E-02 | 2.51 | 0.943E+00 | 1.14 | 0.508E-01 | 1.82 |
| On cuboid grids (Figure [f-g3]) for [sol2] | ||||||
| 1 | 0.690E+00 | 0.00 | 0.592E-11 | 0.00 | 0.182E+02 | 0.00 |
| 2 | 0.337E+00 | 1.03 | 0.175E+02 | 0.00 | 0.738E+00 | 4.62 |
| 3 | 0.665E-01 | 2.34 | 0.278E+01 | 2.65 | 0.185E+00 | 1.99 |
| 4 | 0.936E-02 | 2.83 | 0.506E+00 | 2.46 | 0.495E-01 | 1.90 |
| 5 | 0.186E-02 | 2.33 | 0.222E+00 | 1.19 | 0.159E-01 | 1.63 |
| On wedge grids (Figure [f-g4]) for [sol2] | ||||||
| 1 | 0.538E+00 | 0.00 | 0.241E+02 | 0.00 | 0.201E+01 | 0.00 |
| 2 | 0.196E+00 | 1.45 | 0.865E+01 | 1.48 | 0.575E+00 | 1.81 |
| 3 | 0.431E-01 | 2.19 | 0.300E+01 | 1.53 | 0.186E+00 | 1.63 |
| 4 | 0.627E-02 | 2.78 | 0.133E+01 | 1.18 | 0.459E-01 | 2.02 |
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On cuboid grids (Figure [f-g3]) for [s3] | ||||||
| 1 | 0.159E-02 | 0.00 | 0.640E-01 | 0.00 | 0.285E-01 | 0.00 |
| 2 | 0.195E-03 | 3.03 | 0.729E-02 | 3.13 | 0.182E-02 | 3.97 |
| 3 | 0.138E-04 | 3.82 | 0.106E-02 | 2.78 | 0.140E-03 | 3.70 |
| 4 | 0.903E-06 | 3.94 | 0.140E-03 | 2.93 | 0.993E-05 | 3.82 |
| On wedge grids (Figure [f-g4]) for [s3] | ||||||
| 1 | 0.106E-02 | 0.00 | 0.661E-01 | 0.00 | 0.178E-01 | 0.00 |
| 2 | 0.147E-03 | 2.85 | 0.110E-01 | 2.59 | 0.161E-02 | 3.47 |
| 3 | 0.102E-04 | 3.85 | 0.196E-02 | 2.49 | 0.119E-03 | 3.76 |
| On cuboid grids (Figure [f-g3]) for [sol2] | ||||||
| 1 | 0.103E+01 | 0.00 | 0.950E+02 | 0.00 | 0.808E+01 | 0.00 |
| 2 | 0.254E+00 | 2.01 | 0.167E+02 | 2.50 | 0.156E+01 | 2.37 |
| 3 | 0.291E-01 | 3.13 | 0.204E+01 | 3.04 | 0.113E+00 | 3.78 |
| 4 | 0.383E-02 | 2.93 | 0.411E+00 | 2.31 | 0.941E-02 | 3.59 |
| On wedge grids (Figure [f-g4]) for [sol2] | ||||||
| 1 | 0.578E+00 | 0.00 | 0.546E+02 | 0.00 | 0.415E+01 | 0.00 |
| 2 | 0.117E+00 | 2.31 | 0.108E+02 | 2.33 | 0.372E+00 | 3.48 |
| 3 | 0.100E-01 | 3.55 | 0.352E+01 | 1.62 | 0.434E-01 | 3.10 |
| \(G_i\) | \(\| u-u_h\|_{0}\) | \(O(h^r)\) | \(\|\Delta_w( u-u_h)\|_0\) | \(O(h^r)\) | \(\| p-p_h\|_{0}\) | \(O(h^r)\) |
|---|---|---|---|---|---|---|
| On cuboid grids (Figure [f-g3]) for [sol2] | ||||||
| 1 | 0.480E+01 | 0.00 | 0.421E+03 | 0.00 | 0.686E+02 | 0.00 |
| 2 | 0.393E+00 | 3.61 | 0.462E+02 | 3.19 | 0.126E+01 | 5.76 |
| 3 | 0.162E-01 | 4.60 | 0.542E+01 | 3.09 | 0.600E-01 | 4.40 |
| On wedge grids (Figure [f-g4]) for [sol2] | ||||||
| 1 | 0.719E+00 | 0.00 | 0.106E+03 | 0.00 | 0.739E+01 | 0.00 |
| 2 | 0.961E-01 | 2.90 | 0.569E+02 | 0.90 | 0.212E+00 | 5.12 |
| 3 | 0.492E-02 | 4.29 | 0.807E+01 | 2.82 | 0.113E-01 | 4.23 |
| On tetrahedral grids (Figure [f-g5]) for [sol2] | ||||||
| 1 | 0.881E-01 | 0.00 | 0.744E+01 | 0.00 | 0.508E+00 | 0.00 |
| 2 | 0.919E-02 | 3.26 | 0.302E+01 | 1.30 | 0.122E+00 | 2.06 |
| 3 | 0.349E-02 | 1.40 | 0.432E+00 | 2.81 | 0.978E-02 | 3.64 |