Pointwise superconvergence of block finite elements for the three-dimensional variable coefficient elliptic equation

1.   Introduction

In the realm of solving differential equations using the finite element method (FEM), it has been observed that the rate of convergence of finite element solutions at specific exceptional points within a domain surpasses the optimal global rate. This phenomenon is known as superconvergence, which has already aroused many scholars' interest. In recent decades, superconvergence has become a significant topic in the research field of the Galerkin FEM. With the advancement of research technologies, numerous superconvergence results have been obtained, and theoretical frameworks for superconvergence have also been established. Currently, several important works related to FEM superconvergence are cited in Ref. ^[1,2,3,4,5]. Depending on the partition types within a domain, the commonly used three-dimensional finite elements mainly include tetrahedral, pentahedral, and hexahedral elements. Substantial progress has been made in studying the superconvergence of these three-dimensional FEMs with numerous superconvergence results documented in various published papers or reports such as ^{[6,7,8,9,10,11]}. Recently, we also obtained some superconvergence results for the three-dimensional FEM ^{[12,13,14,15]}. This paper focuses on the superconvergence of the block finite element for variable coefficient elliptic equations, which is a topic not explored by others. The objective is to demonstrate that the error convergence rates of the finite element approximation and corresponding interpolant in the $W^{1, \infty}$-semi-norm and $L^\infty$-norm are one order (or nearly one order) higher than those of the finite element approximation and the true solution (referred to as optimal global rates). It is important to note that the results presented here are generalizations of the research discussed in ^[14].

In the paper, the letter $C$ is employed to represent a generic constant, which may vary in different instances. Additionally, standard notations for the Sobolev spaces and their norms are utilized.

The model problem considered is

where $\Omega\subset{\cal R}^3$ is a rectangular block with boundary, $\partial\Omega$, consisting of faces parallel to the $x$-, $y$-, and $z$-axes. For simplicity, we assume $a_{ij}, a_i, a_0$, and $f$ are sufficiently smooth given functions, and write $\partial_1u = \frac{\partial u}{\partial x}$, $\partial_2u = \frac{\partial u}{\partial y}$, and $\partial_3u = \frac{\partial u}{\partial z}$.

Thus, the standard variational formulation of the problem (1.1) is as follows.

where

We also write

The existence and uniqueness of the solution to (1.2) is given by the Lax-Milgram lemma, see Ciarlet [2, Theorem 1.1.31].

In order to discretize the problem (1.2), one proceeds as follows. The domain $\Omega$ is first partitioned into a regular family of rectangular blocks ${\cal T}^h$ with mesh size $h\in (0, \, 1)$ such that $\bar\Omega = \cup_{e\in {\cal T}^h}\bar e$. Then we define the finite dimensional subspace, $S^h_0(\Omega)\subset H^1_0(\Omega)$, as the standard tensor-product $m$-order finite element space over the partition. Thus, the discrete problem of approximating (1.2) is

Obviously, from (1.2) and (1.4), the following Galerkin orthogonality relation holds.

To obtain the desired results, for every $Z\in\Omega$, and any directional unit vector $\ell\in {\cal R}^3$, we also need the discrete Green's function $G^h_Z$ and discrete derivative Green's function $\partial_{Z, \ell}G^h_Z$ defined by

and

Here, $\partial_{Z, \ell}G^h_Z$ and $\partial_\ell v(Z)$ stand for the following one-sided directional derivatives, respectively.

As for $G^h_Z$ and $\partial_{Z, \ell}G^h_Z$, we have ^[3]

where $\|G^h_Z\|^h_{2, 1, \Omega} = \sum_{e\in {\cal T}^h}\|G^h_Z\|_{2, 1, e}$ and $\|\partial_{Z, \ell}G^h_Z\|^h_{2, 1, \Omega} = \sum_{e\in {\cal T}^h}\|\partial_{Z, \ell}G^h_Z\|_{2, 1, e}$.

The rest of the paper is arranged as follows. In Section 2, for the second-order elliptic equation with variable coefficients, we discuss two weak estimates for the finite element, which are crucial in the superconvergence analysis. Combined with (1.8)–(1.10), several superconvergence results of the finite element approximation are given in Section 3.

2.   Weak estimates for the finite element

In this section, using the properties of the interpolation operator of projection type, we derive the weak estimates.

We write an element

Let $\{l_j(x)\}^\infty_{j = 0}$, $\{\tilde{l}_j(y)\}^\infty_{j = 0}$, and $\{\bar{l}_j(z)\}^\infty_{j = 0}$ be the normalized orthogonal Legendre polynomial systems on $L^2(I_1)$, $L^2(I_2)$, and $L^2(I_3)$, respectively. It is easy to check that $\{l_i(x)\tilde{l}_j(y)\bar{l}_k(z)\}^\infty_{i, j, k = 0}$ is the normalized orthogonal polynomial system on $L^2(e)$. Set

which are called Lobatto functions. Suppose $u\in H^3(e)$. Then, we have the following expansion (see ^[14]):

where the coefficients $\beta_{ijk}$ can be seen in ^[3,14] and satisfy, for $i, j, k\geq1$,

We introduce the standard tensor-product $m$-order polynomial spaces denoted by $T_m$, namely,

Define the interpolation operator of projection type by $\Pi_m^e$: $H^3(e)\rightarrow T_m(e)$ such that

Thus, combining (2.2) and (2.4) yields

Further, we may define the global tensor-product $m$-order interpolation operator of projection type ^[14]

where $(\Pi_m u)|_e = \Pi_m^eu$.

Theorem 2.1. Suppose $\{{\cal T}^h\}$ is a regular family of rectangular partitions of $\Omega$, and $u\in W^{m+2, \, \infty}(\Omega)\cap H^1_0(\Omega)$. Then, for all $v\in S^h_0(\Omega)$, the interpolation operator $\Pi_m$, defined by (2.6), satisfies the following weak estimates:

where $\|v\|^h_{2, \, 1, \, \Omega} = \sum_{e\in {\cal T}^h}\|v\|_{2, \, 1, \, e}$.

Proof. Note that the term $u-\Pi_m^eu$ can be written as follows.

Now, we first bound $(a_{11}\partial_1(u-\Pi_mu), \, \partial_1v) = \sum_e(a_{11}\partial_1(u-\Pi_m^eu), \, \partial_1v)_e$. For $\lambda_{(m+1)00}$, by the orthogonality of Legendre functions, we have

Set $h_e = \hbox{diam}(e)$. Thus,

By (2.3), (2.10), (2.11), and the properties of Legendre functions, we get

It is easy to see

Consider the main parts of $R$ in (2.9), i.e.,

As the symmetry, only need to discuss $\lambda_{01(m+1)}, \lambda_{1(m+1)0}, \lambda_{(m+1)10}$. Obviously,

For $\lambda_{1(m+1)0}$,

By (2.3), (2.15), and the properties of Legendre and Lobatto functions, we have

Similarly,

From (2.9), (2.12)–(2.14), (2.16), and (2.17),

Summing over all elements yields

Similarly to the arguments of the result (2.18),

Now, we bound the terms $(a_{ij}\partial_i(u-\Pi_mu), \, \partial_jv), i\neq j$. Without loss of generality, we consider $(a_{12}\partial_1(u-\Pi_mu), \, \partial_2v) = \sum_e(a_{12}\partial_1(u-\Pi_m^eu), \, \partial_2v)_e$. Nevertheless, from (2.9),

Integration by parts results in

Combined with (2.3), we get

As for $B_e$, we need to apply the element canceling technique. At the adjacent element $e' = (x_e-n_e, x_e+n_e)\times(y_{e'}-k_{e'}, y_{e'}+k_{e'})\times(z_e-d_e, z_e+d_e)$ of $e$, where $k_{e'}-k_e = O(h_e) = O(h_{e'})$ and $y_{e'} < y_e$,

Obviously, there is the same integration factor between $B_e$ and $B_{e'}$, the summing of which is

Here,

Additionally,

By the inverse estimate, we have

Combined with (2.22) and (2.23), we obtain

From (2.21) and (2.24), summing over all elements yields

As for the main parts of $\partial_1R$ in (2.20), it is easy to see

For $\lambda_{1(m+1)0}, \lambda_{10(m+1)}, \lambda_{(m+1)10}, \lambda_{(m+1)01}$, by (2.3) and the properties of Legendre and Lobatto functions, we immediately obtain

From (2.25)–(2.30), and summing over all elements, we obtain

Hence,

It remains to bound the terms $(a_i\partial_i(u-\Pi_mu), \, v)$ and $(a_0(u-\Pi_mu), \, v)$. Using the integration by parts, the element canceling technique, and the interpolation error estimate, we have

Obviously, by the interpolation error estimate, we immediately get

The desired result (2.7) follows from (2.18), (2.19), and (2.32)–(2.34).

Let's prove (2.8) when $m\geq2$. From (2.20), we first bound $(a_{11}\partial_1\lambda_{(m+1)00}, \, \partial_1v)_e$. By the integration by parts and the orthogonality of Lobatto functions, we have

Combined with (2.3), (2.11), and the properties of Lobatto functions,

Now consider the main parts of $\partial_1R$ in (2.20). As the symmetry, we only need to discuss $\partial_1\lambda_{01(m+1)}, \partial_1\lambda_{1(m+1)0}, \partial_1\lambda_{(m+1)10}$. For $\partial_1\lambda_{1(m+1)0}$, integration by parts yields

By (2.3) and the properties of Legendre and Lobatto functions, we have

For $\partial_1\lambda_{(m+1)10}$, integration by parts yields

By (2.3) and the properties of Legendre and Lobatto functions again, we get

From (2.14), (2.35)–(2.37), and summing over all elements,

Similarly,

Next, we bound the terms $(a_{ij}\partial_i(u-\Pi_mu), \, \partial_jv), \, i\neq j$. Applying integration by parts twice, we have

Combined with (2.3) and the properties of Lobatto functions,

At the adjacent element $e' = (x_e-n_e, x_e+n_e)\times(y_{e'}-k_{e'}, y_{e'}+k_{e'})\times(z_e-d_e, z_e+d_e)$ of $e$, where $k_{e'}-k_e = O(h_e) = O(h_{e'})$ and $y_{e'} < y_e$,

Obviously, there is the same integration factor between $D_e$ and $D_{e'}$, the summing of which is

Similarly to the arguments of (2.24),

From (2.40) and (2.41), summing over all elements yields

Obviously,

Next, we discuss $\partial_1\lambda_{1(m+1)0}, \partial_1\lambda_{10(m+1)}, \partial_1\lambda_{(m+1)10}$, and $\partial_1\lambda_{(m+1)01}$. By the integration by parts and the properties of Legendre and Lobatto functions, we immediately obtain

From (2.20) and (2.42)–(2.47),

Hence,

It remains to bound

As for $(a_1\partial_1(u-\Pi_mu), \, v) = \sum_e(a_1\partial_1(u-\Pi_m^eu), \, v)_e$, we first consider $(a_1\partial_1\lambda_{(m+1)00}, \, v)_e$. Integration by parts twice yields

From (2.3) and the properties of Lobatto functions,

In addition,

Utilizing the integration by parts and the properties of Legendre and Lobatto functions, we get

From (2.20), (2.49)–(2.54), and summing over all elements,

Thus,

For $(a_0(u-\Pi_mu), \, v)$, integration by parts results in

From (2.3) and the properties of Lobatto functions,

Similarly,

Obviously,

From (2.9), (2.56)–(2.59), and summing over all elements,

The desired result (2.8) follows from (2.38), (2.39), (2.48), (2.55), and (2.60).

3.   Superconvergence of the finite element

In this section, we will use the weak estimates and the estimates for discrete Green's function and discrete derivative Green's function to obtain superconvergence estimates of the finite element.

Theorem 3.1. Suppose $\{{\cal T}^h\}$ is a regular family of rectangular partitions of $\Omega$, $u\in W^{m+2, \, \infty}(\Omega)\cap H^1_0(\Omega)$, and $u_h\in S^h_0(\Omega)$ and $\Pi_mu\in S^h_0(\Omega)$ are the finite element approximation and interpolant of projection type to $u$, respectively. Then, we have

Proof. For every $Z\in\Omega$, applying the definitions of $G^h_Z$ and $\partial_{Z, \ell}G^h_Z$ as well as the Galerkin orthogonality relation (1.5), we derive

and

From (1.8), (2.7), and (3.4), the result (3.1) is obtained. Combining (1.9), (2.8), and (3.4) yields the result (3.2). The result (3.3) is immediately proved by using (1.10), (2, 8), and (3.5).

Example 3.1. Consider the following equation:

where $a = e^{x+y+z}$, and the exact solution is $u = \sin{\pi x}\sin{\pi y}\sin{\pi z}.$ Let $u_h$ and $\Pi_2$ be the tensor-product quadratic finite element approximation and the interpolation operator of projection type, respectively. We solve Example 3.1 and obtain the following numerical results (see Table 1):

The numerical results demonstrate our theoretical results.

Comments. The domain treated in the paper is a rectangular block in ${\cal R}^3$, which is also discussed by Goodsell ^[10]. Actually, as the Brandts and M. Křížek discussed ^[7], the results of the paper hold for a bounded polyhedral domain, which is usually presented in engineering problems.

4.   Conclusions

In this paper, we generalized superconvergence results of the FEM from constant coefficient elliptic equations to variable coefficient settings in three dimensions. Applying the properties of interpolation operator of projection type, we obtained the weak estimates. Combined with the estimates for discrete Green's function and discrete derivative Green's function, superconvergence results were derived. Among the arguments, how to deal with the given variable coefficients is a challenging issue. The methods presented in the paper can also be applied to other high-dimensional second-order variable coefficient elliptic equations.

Author contributions

The first author proposed ideas of the study, and gave the theoretical analysis of Section 2 and Section 3. The second author mainly focused on the collection of literatures and the numerical example. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work was supported by the Special Projects in Key Fields of Guangdong Province Ordinary Universities (2022ZDZX3016), the Characteristic Innovation Project of Guangdong Province Ordinary University (2023KTSCX089), and the Projects of Talents Recruitment of Guangdong University of Petrochemical Technology (2021rc003).

Conflict of interest

The authors declare no conflicts of interest.