A Posteriori Error Estimate of Two Grid Mixed Finite Element Methods for Semilinear Elliptic Equations

Yiming Wen, Luoping Chen^*, Jiajia Dai

School of Mathematics, Southwest Jiaotong University, Chengdu, China.

DOI: 10.4236/jamp.2023.112020                

Abstract

In this paper, we present the a posteriori error estimate of two-grid mixed finite element methods by averaging techniques for semilinear elliptic equations. We first propose the two-grid algorithms to linearize the mixed method equations. Then, the averaging technique is used to construct the a posteriori error estimates of the two-grid mixed finite element method and theoretical analysis are given for the error estimators. Finally, we give some numerical examples to verify the reliability and efficiency of the a posteriori error estimator.

Keywords

Two-Grid Mixed Finite Element Methods, Posteriori Error Estimates, Semilinear Elliptic Equations, Averaging Technique

1. Introduction

Mixed finite element method is a kind of method in solving partial differential equations (PDEs). Mixed methods are based on writing a higher order differential equation into lower order differential system. The purpose of this article is to study the a posteriori error estimate of the mixed finite element methods for the following semilinear elliptic equation

$-\Delta u=f\left(x\mathrm{,}u\right)\mathrm{,}\text{in}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\Omega \mathrm{,}$ (1.1)

with mixed boundary conditions. $\Omega \subset R^2$ is a convex polygonal domain and $f\left(x\mathrm{,}u\right)$, a real-valued function on $\Omega$, has continuous first and second derivatives to u.

The a posteriori error estimates of mixed finite element method have been studied extensively in the past several decades for solving many differential model problems, for example, the Navier-Stokes equations based on Newton-type linearization by Durango and Novo [1], the linear elliptic problems by Larson and Målqvit [2], the Poisson problem about an error estimate in the $H\left(\text{div}\mathrm{,}\Omega \right)$ norm of the flux by Carstensen [3], the general convex optimal control problems by Chen and Liu [4]. In order to combine the advantage of adaptive mixed finite element method and the efficiency of two-grid finite element method for semilinear elliptic equations, in this study, we proposed the posteriori error estimator for the two-grid mixed finite element methods.

The two grid method is a widely used numerical method in solving nonlinear problems. It was first introduced by Xu [5] [6] to solve the nonsymmetric linear and nonlinear elliptic problems. Many numerical methods combined with two-grid method were used to solve different model problems, for instance, nonlinear reaction-diffusion equations using mixed finite element methods by Chen and Chen [7], nonlinear parabolic equations by Chen and Liu [8], the coupled Stokes-Darcy system by Sun, Shi, et al. [9], two-dimensional semi-linear elliptic interface problems by Chen, Li, et al. [10]. In recent years, the residual-based a posteriori error estimates of two-grid finite element methods and finite volume methods are investigated for nonlinear PDEs [11] [12]. Adaptive two-grid finite element methods based on residual-based a posteriori error estimator are studied in [13].

In order to investigate efficient two-grid adaptive mixed finite element method for semilinear or nonlinear PDEs, in this paper, we study two-grid mixed finite element method and its posteriori error estimates for semilinear elliptic problem (1.1). We first propose two algorithms for the model problem. Then, for both two-grid mixed finite element methods, by using averaging technique, the posteriori error estimators are proposed for the flux error in L²-norm. Theoretical analysis is given to prove the efficiency and reliability of the error estimators. Two numerical examples are given to verify the theoretical results and from the numerical results, we find that the error estimators proposed in this paper are efficient and reliable.

The outline of this paper is organized as follows. In Section 2, we present some notations and weak form of the semilinear elliptic Equations (1.1). Two-grid mixed finite element methods for the model problem are presented in Section 3. In Section 4, we give a theoretical analysis of the reliability and efficiency for the posteriori error estimators. Numerical experiments are given to verify the theoretical results in Section 5.

2. Weak Form and Preliminaries

In this section, we will present some preliminaries and weak form for the semilinear model problem (1.1).

2.1. Preliminary

We first introduce the standard notations used in this paper. We denote $W^{m\mathrm{,}p}\left(\Omega \right)=\left\{v\in L^p\left(\Omega \right)\mathrm{,}{D}^{\alpha }\left(v\right)\in L^p\left(\Omega \right)\mathrm{,}\forall \alpha \in Z_{+}^n\mathrm{,}\left|\alpha \right|\le m\right\}$ as Sobolev spaces with the norm ${\Vert \cdot \Vert }_{m\mathrm{,}p}$, for integer $m\ge 0$ and real number $1\le p\le \infty$, ${\Vert v\Vert }_{m\mathrm{,}p}^p={\displaystyle {\sum }_{\left|\alpha \right|\le m}}{\Vert D^{\alpha }v\Vert }_{L^p\left(\Omega \right)}^p$ and $W_0^{m,p}\left(\Omega \right)=\left\{v\in W^{m,p}\left(\Omega \right):{v|}_{\partial \Omega }=0\right\}$. When $p=2$, we denote $W^{m\mathrm{,2}}\left(\Omega \right)$ by $H^m\left(\Omega \right)$, $W_0^{m\mathrm{,2}}\left(\Omega \right)$ by $H_0^m\left(\Omega \right)$, and we will use ${\Vert \cdot \Vert }_{m\mathrm{,2}}={\Vert \cdot \Vert }_m$ and ${\Vert \cdot \Vert }_{\mathrm{0,2}}=\Vert \cdot \Vert$.

Throughout this paper, we will use letter C to denote a generic positive constant that may represent different values at different places.

2.2. Weak Form

In order to introduce a mixed variational formulation on $\Omega$, we first introduce the following spaces

$V=H\left(\text{div}\mathrm{,}\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{v\in L^2{\left(\Omega \right)}^2\mathrm{,}\text{div}v\in L^2\left(\Omega \right)\right\}\mathrm{,}$

$W=L^2\left(\Omega \right),$

$H_D^1\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{u\in H^1\left(\Omega \right)\mathrm{,}{u|}_{{\Gamma }_D}=0\right\}\mathrm{,}$

with the norms

${\Vert v\Vert }_V={\left({\Vert v\Vert }^2+{\Vert \text{div}v\Vert }^2\right)}^{1/2}\text{and}{\Vert q\Vert }_W=\Vert q\Vert \mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}v\in V\mathrm{,}q\in W\mathrm{.}$

The Lipschitz boundary $\Gamma =\partial \Omega$ of the bounded domain $\Omega$ is split into a closed Dirichlet part ${\Gamma }_D$ and possibly empty Neumann part ${\Gamma }_N$. Set $p=-\nabla u$. Rewrite the problem (1.1), we have

$\{\begin{array}{ll}\text{div}p=f\left(u\right)\mathrm{,}\hfill & \text{in}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\Omega \mathrm{,}\hfill \\ p=-\nabla u\mathrm{,}\hfill & \text{in}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\Omega \mathrm{,}\hfill \\ u=\mathrm{0,}\hfill & \text{on}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\Gamma }_D\mathrm{,}\hfill \\ p\cdot n=g_N\mathrm{,}\hfill & \text{on}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\Gamma }_N\mathrm{.}\hfill \end{array}$ (2.1)

Here $g_N$ is given function.

We define a space $V^0=\left\{v\in V\mathrm{<mo>\; :\; </mo>}v\cdot n=\mathrm{0,}\text{on}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\Gamma }_N\right\}$, the standard mixed variational form of (2.1) is to find $\left(p\mathrm{,}u\right)\in V\times W$, such that

$\left(p\mathrm{,}v\right)-\left(u\mathrm{,}\text{div}v\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}v\in V^0\mathrm{,}$ (2.2)

$\left(\text{div}p\mathrm{,}w\right)=\left(f\left(u\right)\mathrm{,}w\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}w\in W\mathrm{.}$ (2.3)

where $\left(\cdot \mathrm{,}\cdot \right)$ is the inner product of $L^2\left(\Omega \right)$.

Given $\mu \in W$, the linearised form of (2.2) and (2.3) is

$\left(p\mathrm{,}v\right)-\left(u\mathrm{,}\text{div}v\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}v\in V^0\mathrm{,}$ (2.4)

$\left(\text{div}p\mathrm{,}w\right)=\left(f\left(\mu \right)+f^{\prime }\left(\mu \right)\left(u-\mu \right)\mathrm{,}w\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}w\in W\mathrm{.}$ (2.5)

Let ${\mathcal{T}}_h$ denote a regular triangulation of the polygonal domain $\Omega$, $h_T$ denotes the diameter of the element $T\in {\mathcal{T}}_h$ and $h=\max \left\{h_T\right\}$. And for $T\in {\mathcal{T}}_h$, let $P_T^k={\mathcal{P}}_k\left(T\right)$ denote the set of algebraic polynomials in $d\left(d=2\right)$ variables on T of total degree ≤ k.

The set of all nodes and edges appearing in ${\mathcal{T}}_h$ are denoted as $\mathcal{N}$ and $\epsilon$. $E\in {\epsilon }_D$ denotes edges $E\in \epsilon$ on the boundary ${\Gamma }_D$, $E\in {\epsilon }_N$ denotes edges $E\in \epsilon$ on the boundary ${\Gamma }_N$, $E\in {\epsilon }_I$ denotes edges $E\in \epsilon$ but $E\cancel{\subset }\partial \Omega$.

The space ${\mathcal{L}}^k\left({\mathcal{T}}_h\right)$ (possibly discontinuous) of ${\mathcal{T}}_h$ -piecewise polynomials of degree ≤ k is the set of all $U\in L^{\infty }\left(\Omega \right)$ (a set composed of all bounded number columns) with ${U|}_T\in P_T^k$ for all T in ${\mathcal{T}}_h$. Set

$V\left(T\right)={\left[{\mathcal{P}}_0\left(T\right)\right]}^d\oplus \text{span}\left(x{\mathcal{P}}_0\left(T\right)\right)\mathrm{,}$

$W\left(T\right)={\mathcal{P}}_0\left(T\right)\mathrm{,}$

$S^k\left({\mathcal{T}}_h\right)={\mathcal{L}}^k\left({\mathcal{T}}_h\right)\cap C\left(\Omega \right)\text{and}{S}_D^k\left({\mathcal{T}}_h\right)\mathrm{<mo>\; :\; </mo>}=\left\{u_h\in S^k\left({\mathcal{T}}_h\right)\mathrm{<mo>\; :\; </mo>}{u_h|}_{{\Gamma }_D}=0\right\}\mathrm{.}$

Here $C\left(\Omega \right)$ denotes continuous space.

Let

$\begin{array}{l}{V}_h\mathrm{<mo>\; :\; </mo>}=\left\{v\in V\mathrm{<mo>\; :\; </mo>}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}T\in {\mathcal{T}}_h\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>}{v|}_T\in V\left(T\right)\right\},\\ W_h=\left\{w\in W\mathrm{<mo>\; :\; </mo>}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}T\in {\mathcal{T}}_h\mathrm{,}{w|}_T\in W\left(T\right)\right\}\mathrm{.}\end{array}$

In this paper, we mainly study $V_h\times W_h\subset V\times W$ as the lowest order Raviart-Thomas mixed finite element spaces for the discretization of the flux $p$ and u, we define $V_h^0=\left\{v\in V_h\mathrm{<mo>\; :\; </mo>}v\cdot n=\mathrm{0,}\text{on}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{\Gamma }_N\right\}$ ( $V_H^0$ has the same definition as $V_h^0$ ). Therefore, the discretization of mixed finite element method is to find $\left(p_h\mathrm{,}{u}_h\right)\in V_h\times W_h$ such that

$\left(p_h\mathrm{,}{v}_h\right)-\left(u_h\mathrm{,}\text{div}{v}_h\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v}_h\in V_h^0\mathrm{,}$ (2.6)

$\left(\text{div}{p}_h\mathrm{,}{w}_h\right)=\left(f\left(u_h\right)\mathrm{,}{w}_h\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{w}_h\in W_h\mathrm{.}$ (2.7)

2.3. Helmholtz Decomposition and Interpolation Operator

In order to make theoretical analysis, we need to introduce Helmholtz decomposition and the interpolation operator $\mathcal{J}$ [14].

We first define the curl operator as follows [15],

$\text{if}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\tau \mathrm{<mo>\; :\; </mo>}\Omega \to R^3\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{curl}\tau =\nabla \times \tau =\left(\begin{array}{c}\frac{\partial {\tau }_3}{\partial x_2}-\frac{\partial {\tau }_2}{\partial x_3}\\ \frac{\partial {\tau }_1}{\partial x_3}-\frac{\partial {\tau }_3}{\partial x_1}\\ \frac{\partial {\tau }_2}{\partial x_1}-\frac{\partial {\tau }_1}{\partial x_2}\end{array}\right)\mathrm{,}$

$\text{if}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\tau \mathrm{<mo>\; :\; </mo>}\Omega \to R^2\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{curl}\tau =\frac{\partial {\tau }_2}{\partial x_1}-\frac{\partial {\tau }_1}{\partial x_2}\mathrm{,}$

$\text{if}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}v\mathrm{<mo>\; :\; </mo>}\Omega \to R\mathrm{,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}d=\mathrm{2,\; <mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{curl}v=\left(\begin{array}{c}\frac{\partial v}{\partial x_2}\\ -\frac{\partial v}{\partial x_1}\end{array}\right)\mathrm{.}$

Then, we can get $\text{div}\left(\text{curl}\phi \right)=0$ and the Gauss theorem yields the following relation

${\displaystyle {\int }_{\Omega }}\text{curl}\phi \cdot \tau \text{d}x={\displaystyle {\int }_{\Omega }}\phi \cdot \text{curl}\tau \text{d}x+{\displaystyle {\int }_{\partial \Omega }}\left(\tau -\left(\tau \cdot n\right)n\right)\cdot \phi \text{d}s\mathrm{.}$ (2.8)

here $n$ is outer unit normal vector of $\partial \Omega$, and $\tau -\left(\tau \cdot n\right)n$ denotes the tangential component of $\tau$.

We define an approximation operator $\mathcal{J}\mathrm{<mo>\; :\; </mo>}{H}_D^1\left(\Omega \right)\to S_D^1\left({\mathcal{T}}_h\right)$. Let $\left\{{{\phi }_z|}_{z\in \mathcal{N}}\right\}$ denote the nodal basis of $S^1\left({\mathcal{T}}_h\right)$, ${\phi }_z\in S^1\left({\mathcal{T}}_h\right)$ satisfies ${\phi }_z\left(x\right)=0$ if $x\in \mathcal{N}\backslash \left\{z\right\}$ and ${\phi }_z\left(z\right)=1$, the open patches defined by $w_z=\left\{x\in \Omega :0<{\phi }_z\left(x\right)\right\}$.

Then, we modify $\left\{{{\phi }_z|}_{z\in \mathcal{N}}\right\}$ to be a partition of unity $\left({\psi }_z\mathrm{<mo>\; |\; </mo>}z\in \mathcal{K}\right)$ ( $\mathcal{K}=\mathcal{N}\backslash \partial \Omega$ denotes the set of free nodes). Find each fixed node $z\in \mathcal{N}\backslash \mathcal{K}$, we choose a node $\zeta \left(z\right)\in \mathcal{K}$ and let $\zeta \left(z\right)=z$ if $z\in \mathcal{K}$. In this way, we define a partition of $\mathcal{N}$ into $\text{card}\left(\mathcal{K}\right)$ classes $I\left(z\right)=\left\{\tilde{z}\in \mathcal{N}\mathrm{<mo>\; :\; </mo>}\zeta \left(\tilde{z}\right)=z\right\}$, where $z\in \mathcal{K}$. For each $z\in \mathcal{K}$, set

${\psi }_z={\displaystyle \underset{\zeta \in I\left(z\right)}{\sum }}{\phi }_{\zeta }\mathrm{,}$

and notice that $\left({\psi }_z\mathrm{<mo>\; |\; </mo>}z\in \mathcal{K}\right)$ is a partition of unity. It is required that

${\Omega }_z=\left\{x\in \Omega :0<{\psi }_z\left(x\right)\right\},$

is connected.

For $g\in L^1\left(\Omega \right)$ and $z\in \mathcal{K}$, let $g_z\in R$ be

$g_z=\frac{{\displaystyle {\int }_{{\Omega }_z}}g{\psi }_z\text{d}x}{{\displaystyle {\int }_{{\Omega }_z}}{\phi }_z\text{d}x}\mathrm{,}$

and then define

$\mathcal{J}g={\displaystyle \underset{z\in \mathcal{K}}{\sum }}{g}_z{\phi }_z\mathrm{.}$

We also define local mesh-sizes by $h_{\mathcal{T}}$ and $h_{\epsilon }$, where $h_{\mathcal{T}}$ denotes the element-size, ${h_{\mathcal{T}}|}_T=h_T$ for $T\in {\mathcal{T}}_h$, and the edge-size ${h_{\epsilon }|}_E=h_E=\text{diam}\left(E\right)$.

We also use the orthogonal $L^2$ -projection $P_h$ [16] : $W\to W_h$, which satisfies

$\Vert P_{h}w-w\Vert \le Ch{\Vert w\Vert }_1\mathrm{.}$ (2.9)

Lemma 2.1 ( [14] ). There exist $\left(h_{\mathcal{T}}\mathrm{,}{h}_{\epsilon }\right)$ -independent constant C such that for all $g\in H_D^1\left(\Omega \right)$ and $f\in L^2\left(\Omega \right)$, there holds

$\Vert \nabla \mathcal{J}g-\nabla g\Vert \le C\Vert \nabla g\Vert \mathrm{,}$ (2.10)

${\displaystyle {\int }_{\Omega }}f\left(g-\mathcal{J}g\right)\text{d}x\le C\Vert \nabla g\Vert {\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert f-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\mathrm{,}$ (2.11)

$\Vert h_{\mathcal{T}}^{-1}\left(g-\mathcal{J}g\right)\Vert \le C\Vert \nabla g\Vert \mathrm{,}$ (2.12)

${\Vert h_{\epsilon }^{-1/2}\left(g-\mathcal{J}g\right)\Vert }_{L^2\left({\Gamma }_N\right)}\le C\Vert \nabla g\Vert \mathrm{,}$ (2.13)

the constants C only depend on $\Omega$ , ${\Gamma }_D$, ${\Gamma }_N$ and the shape of the elements and patches.

3. Two-Grid Finite Element Methods for Semilinear Problems

In this section, we present two-grid mixed finite element methods for the semilinear elliptic problems and analyze the lower and upper bounds of posteriori error estimates by averaging techniques.

The idea of two-grid methods is to solve the semilinear partial differential equations(PDEs) on the coarse mixed finite element spaces $V_H\times W_H\subset V_h\times W_h$ first and then find the solution ( $p^h\mathrm{,}{u}^h$ ) (or ( ${\tilde{p}}^h\mathrm{,}{\tilde{u}}^h$ )) of a linear PDEs on the finer mixed finite element spaces $V_h\times W_h$. The basic mechanism in these algorithms is to construct two shape-regular subdivision of $\Omega$ as ${\mathcal{T}}_H$ and ${\mathcal{T}}_h$ with different mesh sizes H and h ( $h\ll H$ ).

Two-grid Algorithm 1

Step 1: On the coarse mesh ${\mathcal{T}}_H$, compute $\left(p_H\mathrm{,}{u}_H\right)\in V_H\times W_H$ to satisfy the following original nonlinear system:

$\left(p_H\mathrm{,}{v}_H\right)-\left(u_H\mathrm{,}\text{div}{v}_H\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v}_H\in V_H^0\mathrm{,}$ (3.1)

$\left(\text{div}{p}_H\mathrm{,}{w}_H\right)=\left(f\left(u_H\right)\mathrm{,}{w}_H\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{w}_H\in W_H\mathrm{.}$ (3.2)

Step 2: On the fine grid ${\mathcal{T}}_h$, compute $\left(p^h\mathrm{,}{u}^h\right)\in V_h\times W_h$ to satisfy the following linear system:

$\left(p^h\mathrm{,}{v}_h\right)-\left(u^h\mathrm{,}\text{div}{v}_h\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v}_h\in V_h^0\mathrm{,}$ (3.3)

$\left(\text{div}{p}^h\mathrm{,}{w}_h\right)=\left(f\left(u_H\right)\mathrm{,}{w}_h\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{w}_h\in W_h\mathrm{.}$ (3.4)

The second two-grid algorithm introduces the Newton linearized procedure on the fine mesh to linearize the semilinear system.

Two-grid Algorithm 2

Step 1: On the coarse grid ${\mathcal{T}}_H$, compute $\left(p_H\mathrm{,}{u}_H\right)\in V_H\times W_H$ to satisfy the following original nonlinear system:

$\left(p_H\mathrm{,}{v}_H\right)-\left(u_H\mathrm{,}\text{div}{v}_H\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v}_H\in V_H^0\mathrm{,}$ (3.5)

$\left(\text{div}{p}_H\mathrm{,}{w}_H\right)=\left(f\left(u_H\right)\mathrm{,}{w}_H\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{w}_H\in W_H\mathrm{.}$ (3.6)

Step 2: On the fine grid ${\mathcal{T}}_h$, compute $\left({\tilde{p}}^h,{\tilde{u}}^h\right)\in V_h\times W_h$ to satisfy the following linear system:

$\left({\tilde{p}}^h\mathrm{,}{v}_h\right)-\left({\tilde{u}}^h\mathrm{,}\text{div}{v}_h\right)=\mathrm{0,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{v}_h\in V_h^0\mathrm{,}$ (3.7)

$\left(\text{div}{\tilde{p}}^h\mathrm{,}{w}_h\right)-\left(f^{\prime }\left(u_H\right){\tilde{u}}^h\mathrm{,}{w}_h\right)=\left(f\left(u_H\right)\mathrm{,}{w}_h\right)-\left(f^{\prime }\left(u_H\right)u_H\mathrm{,}{w}_h\right)\mathrm{,}\forall \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}{w}_h\in W_h\mathrm{.}$ (3.8)

For semilinear system (3.1) and (3.2) and (3.5) and (3.6), we use the Newton iteration to compute $\left(p_H\mathrm{,}{u}_H\right)$ in the implementation.

In averaging techniques, the error estimator is based on a smoother approximation in $S^1{\left({\mathcal{T}}_h\right)}^2$, the continuous ${\mathcal{T}}_h$ -piecewise linears approximation to the discrete solution $p^h$ (or ${\tilde{p}}^h$ ), for instance,

${\eta }_z=\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert p^h-q^h\Vert \mathrm{,}$

which can be served as a computable estimator.

The triangle inequality shows that ${\eta }_z$ is efficient up to higher order terms of exact solution $p$, indeed,

$\begin{array}{c}{\eta }_z=\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert p^h-q^h\Vert \\ =\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert p^h-p+p-q^h\Vert \\ \le \Vert p-p^h\Vert +\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert p-q^h\Vert \mathrm{.}\end{array}$

The last term converges as $O\left(h^2\right)$ is of higher order than the error $\Vert p-p^h\Vert =O\left(h\right)$. So we have

${\eta }_z-\text{h}\text{.o}\text{.t}\text{.}\le \Vert p-p^h\Vert \mathrm{.}$ (3.9)

where “h.o.t.” denotes the higher-order term.

In the following, by using the solution ( $p^h\mathrm{,}{u}^h$ ) (or ( ${\tilde{p}}^h\mathrm{,}{\tilde{u}}^h$ )), Helmholtz decomposition and interpolation operator $\mathcal{J}$, we analyze the upper bound of $\Vert p-p^h\Vert$ and $\Vert p-{\tilde{p}}^h\Vert$ for the two algorithms.

3.1. A Upper Bound for the Error of Two-Grid Algorithm 1

By Helmholtz decomposition, we get the following lemma.

Lemma 3.1 ( [14] ). There exist $\alpha \mathrm{,}\beta \in H^1\left(\Omega \right)$ that satisfy boundary condition ${\alpha |}_{{\Gamma }_D}=0$ and ${\beta |}_{{\Gamma }_N}$ is constant

$p-p^h=\nabla \alpha +\text{curl}\beta \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}{\Vert p-p^h\Vert }^2={\Vert \nabla \alpha \Vert }^2+{\Vert \nabla \beta \Vert }^2\mathrm{,}$ (3.10)

and then

${\Vert p-p^h\Vert }^2={\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \alpha \text{d}x+{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\beta \text{d}x\mathrm{.}$ (3.11)

In order to estimate the right hand side of (3.11), and using the theoretical analysis in [17], we have

Lemma 3.2. Suppose the u and $u_H$ are the solutions of (2.2) and (2.3) and (3.1) and (3.2), there exists a constant C independent of H such that

$\Vert u-u_H\Vert \le CH{\Vert u\Vert }_2\mathrm{.}$ (3.12)

We can also get ${\Vert u-u_H\Vert }_{L^4\left(\Omega \right)}\le CH{\Vert u\Vert }_3$.

Then, by using the Green’s fromula and Lemma 3.2, we can bound the first contribution of (3.11).

Lemma 3.3. Let $p$ and $p^h$ are the solutions of (2.2) and (2.3) and (3.3) and (3.4), and $p\cdot n\mathrm{,}{p}^h\cdot n\in L^2\left({\Gamma }_N\right)$ . Then we have

$\begin{array}{c}{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \alpha \text{d}x\le CHh\Vert f^{\prime }\left(\hat{u}\right)\Vert \Vert f\left(u\right)\Vert \Vert \nabla \alpha \Vert +C{h}^2\Vert f\left(u_H\right)\Vert \Vert \nabla \alpha \Vert \\ +C{\Vert h_{\epsilon }^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \mathrm{.}\end{array}$ (3.13)

Proof. Employ the Green’s formula and $L^2$ -projection $P_h$, we can get the first contribution on the right-hand side of (3.11), that

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \alpha \text{d}x={\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \left(\alpha -\hat{\alpha }\right)\text{d}x\\ =-{\displaystyle {\int }_{\Omega }}\text{div}\left(p-p^h\right)\left(\alpha -\hat{\alpha }\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_N}}\left(p-p^h\right)\cdot n\left(\alpha -\hat{\alpha }\right)\text{d}s\\ ={\displaystyle {\int }_{\Omega }}\left(-f\left(u\right)+P_{h}f\left(u_H\right)\right)\left(\alpha -\hat{\alpha }\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_N}}\left(p-p^h\right)\cdot n\left(\alpha -\hat{\alpha }\right)\text{d}s\\ ={\displaystyle {\int }_{\Omega }}\left(-f\left(u\right)+f\left(u_H\right)\right)\left(\alpha -\hat{\alpha }\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(P_{h}f\left(u_H\right)-f\left(u_H\right)\right)\left(\alpha -\hat{\alpha }\right)\text{d}x\\ +{\displaystyle {\int }_{{\Gamma }_N}}\left(p-p^h\right)\cdot n\left(\alpha -\hat{\alpha }\right)\text{d}s\\ =I_1+I_2+I_3\mathrm{.}\end{array}$ (3.14)

Here mean value $\hat{\alpha }\in {\mathcal{L}}^0\left({\mathcal{T}}_h\right)$ of $\alpha \in H_D^1\left(\Omega \right)$. Now we estimate the right-hand side terms. For $I_1$, using Lemma 3.2, we conclude that

$\begin{array}{c}{I}_1={\displaystyle {\int }_{\Omega }}-f^{\prime }\left(\hat{u}\right)\left(u-u_H\right)\left(\alpha -{\alpha }_T\right)\text{d}x\\ \le Ch\Vert f^{\prime }\left(\hat{u}\right)\Vert \Vert u-u_H\Vert \Vert \nabla \alpha \Vert \\ \le CHh\Vert f^{\prime }\left(\hat{u}\right)\Vert {\Vert u\Vert }_2\Vert \nabla \alpha \Vert \mathrm{,}\end{array}$

here $\hat{u}$ is a value between u and $u_H$. For $I_2$, employ the inequality (2.9), and $I_3$ with Lemma 4.2 of [14], we conclude that

$I_2\le C{h}^2\Vert f\left(u_H\right)\Vert \Vert \nabla \alpha \Vert \mathrm{,}$

$I_3\le C{\Vert h_{\epsilon }^{1/2}\left(g_N-p^h\cdot n\right)\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \le C{\Vert h_{\epsilon }^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \mathrm{.}$ (3.15)

which completes proof.

Lemma 3.4. Suppose $p\mathrm{,}q\in V$, $p\cdot n\mathrm{,}{p}^h\cdot n\in L^2\left({\Gamma }_N\right)$, and $p^h\in V_h$ is the two-grid solution satisfying (3.3) and (3.4), for $w\in H_D^1\left(\Omega \right)$ and $\Vert \nabla w\Vert =1$, there holds

$\begin{array}{l}\underset{w\in H_0^1\left(\Omega \right)}{\sup }{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla w\text{d}x\\ \le C\Vert p^h-q\Vert +C{\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{div}\left(p-q\right)-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}+CHh\Vert f^{\prime }\left(\hat{u}\right)\Vert {\Vert u\Vert }_2+C{h}^2\Vert f\left(u_H\right)\Vert \mathrm{.}\end{array}$ (3.16)

Proof. Employ the Green’s formula, then,

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla w\text{d}x\\ ={\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \left(w-\mathcal{J}w\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \mathcal{J}w\text{d}x\\ ={\displaystyle {\int }_{\Omega }}\left(p-q\right)\cdot \nabla \left(w-\mathcal{J}w\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(q-p^h\right)\cdot \nabla \left(w-\mathcal{J}w\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \mathcal{J}w\text{d}x\\ =-{\displaystyle {\int }_{\Omega }}\text{div}\left(p-q\right)\left(w-\mathcal{J}w\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_N}}\left(w-\mathcal{J}w\right)\left(p-q\right)\cdot n\text{d}s\\ +{\displaystyle {\int }_{\Omega }}\left(q-p^h\right)\cdot \nabla \left(w-\mathcal{J}w\right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \mathcal{J}w\text{d}x\mathrm{.}\end{array}$ (3.17)

Employ the inequalities (2.10), (2.11) and (2.13) of Lemma 2.1, the first two terms of the right hand side of (3.17) can be written

$-{\displaystyle {\int }_{\Omega }}\text{div}\left(p-q\right)\left(w-\mathcal{J}w\right)\text{d}x\le C{\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{div}\left(p-q\right)-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\mathrm{,}$ (3.18)

${\displaystyle {\int }_{{\Gamma }_N}}\left(w-\mathcal{J}w\right)\left(p-q\right)\cdot n\text{d}s\le C{\Vert h_{\epsilon }^{1/2}\left(p-q\right)\cdot n\Vert }_{L^2\left({\Gamma }_N\right)}\le C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\mathrm{,}$ (3.19)

${\displaystyle {\int }_{\Omega }}\left(q-p^h\right)\cdot \nabla \left(w-\mathcal{J}w\right)\text{d}x\le C\Vert p^h-q\Vert \mathrm{.}$ (3.20)

And here we have $S^1\left({\mathcal{T}}_h\right)\subset H^1\left(\Omega \right)$, so we estimate the last term ${\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \mathcal{J}w\text{d}x$ just like estimate ${\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \alpha \text{d}x$, and by using triangle inequality. Therefore, we complete the proof.

For the second term on the right hand side (3.11), denote $\mathcal{Q}\left(a_1\mathrm{,}{a}_2\right)=\left(-a_2\mathrm{,}{a}_1\right)$ for vectors, then, $\text{curl}\left(q-p\right)=\text{div}\mathcal{Q}\left(q-p\right)$ and use that $\text{curl}p=0$, let $\bar{p}=\mathcal{Q}p$ and ${\bar{p}}^h=\mathcal{Q}{p}^h$, use ${\displaystyle {\int }_{\Omega }}\left(\bar{p}-{\bar{p}}^h\right)\cdot \nabla w_h\text{d}x=-{\displaystyle {\int }_{\Omega }}\text{curl}\left(p-p^h\right)w_h\text{d}x$ for all $w_h\in {\tilde{S}}_D^1\left({\mathcal{T}}_h\right)=S_N^1\left({\mathcal{T}}_h\right)$, where $S_N^1\left({\mathcal{T}}_h\right)=\left\{v_h\in S^1\left({\mathcal{T}}_h\right)\mathrm{<mo>\; :\; </mo>}{v_h|}_{{\Gamma }_N}=0\right\}$. Therefore, we can get the following lemma.

Lemma 3.5. Let $p$ and $p^h$ are the solutions of (2.2) and (2.3) and (3.3) and (3.4), then,

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\beta \text{d}x\\ \le C\Vert \nabla \beta \Vert \Vert q-p^h\Vert +C\Vert \nabla \beta \Vert {\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}q-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +C\Vert \nabla \beta \Vert {\Vert h_{\epsilon }^{1/2}\left(q\cdot t\right)\Vert }_{L^2\left({\Gamma }_D\right)}+CHh\Vert \nabla \beta \Vert \Vert f^{\prime }\left(\hat{u}\right)\Vert {\Vert u\Vert }_2\\ +C{h}^2\Vert \nabla \beta \Vert \Vert f\left(u_H\right)\Vert +C\Vert \nabla \beta \Vert {\Vert h_{\epsilon }^{1/2}\left(p^h\cdot t\right)\Vert }_{L^2\left({\Gamma }_D\right)}\mathrm{.}\end{array}$ (3.21)

Proof. Using Lemma 3.2 and Lemma 3.4 and (2.8), let $\beta =0$ on ${\Gamma }_N$. So we have

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\beta \text{d}x\\ ={\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\left(\beta -\mathcal{J}\beta \right)\text{d}x+{\displaystyle {\int }_{\Omega }}\mathrm{<mo\; stretchy="false">\; (\; </mo>}p-p^h\mathrm{<mo\; stretchy="false">\; )\; </mo>}\cdot \text{curl}\mathcal{J}\beta \text{d}x\\ ={\displaystyle {\int }_{\Omega }}\left(p-q\right)\cdot \text{curl}\left(\beta -\mathcal{J}\beta \right)\text{d}x+{\displaystyle {\int }_{\Omega }}\left(q-p^h\right)\cdot \text{curl}\left(\beta -\mathcal{J}\beta \right)\text{d}x\\ +{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\mathcal{J}\beta \text{d}x\\ \le {\displaystyle {\int }_{\Omega }}\text{curl}\left(p-q\right)\cdot \left(\beta -\mathcal{J}\beta \right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_D}}\left(\left(p-q\right)-\left(\left(p-q\right)\cdot n\right)n\right)\cdot \left(\beta -\mathcal{J}\beta \right)\text{d}s\\ +C\Vert \nabla \beta \Vert \Vert q-p^h\Vert +{\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\mathcal{J}\beta \text{d}x\mathrm{,}\end{array}$ (3.22)

By using Lemma 2.1, we get

${\displaystyle {\int }_{\Omega }}\text{curl}\left(p-q\right)\cdot \left(\beta -\mathcal{J}\beta \right)\text{d}x\le C\Vert \nabla \beta \Vert {\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}\left(p-q\right)-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\mathrm{,}$ (3.23)

$\begin{array}{l}{\displaystyle {\int }_{{\Gamma }_D}}\left(\left(p-q\right)-\left(\left(p-q\right)\cdot n\right)n\right)\cdot \left(\beta -\mathcal{J}\beta \right)\text{d}s\\ \le C\Vert \nabla \beta \Vert {\Vert h_{\epsilon }^{1/2}\left(\left(p-q\right)-\left(\left(p-q\right)\cdot n\right)n\right)\Vert }_{L^2\left({\Gamma }_D\right)}\\ \le C\Vert \nabla \beta \Vert {\Vert h_{\epsilon }^{1/2}\left(q\cdot t\right)\Vert }_{L^2\left({\Gamma }_D\right)}\mathrm{,}\end{array}$ (3.24)

where $t=\mathcal{Q}n$ denotes unit tangential vector, using $\text{curl}p=0$, (3.23) can be written as

${\displaystyle {\int }_{\Omega }}\text{curl}\left(p-q\right)\cdot \left(\beta -\mathcal{J}\beta \right)\text{d}x\le C\Vert \nabla \beta \Vert {\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}q-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\mathrm{,}$ (3.25)

the last term ${\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \text{curl}\mathcal{J}\beta \text{d}x$ similar to ${\displaystyle {\int }_{\Omega }}\left(p-p^h\right)\cdot \nabla \mathcal{J}w\text{d}x$ of Lemma 3.4, which completes the proof.

Combine with Lemma 3.3 and Lemma 3.5, and using Lemma 5.2 of [14], we can know ${\Vert h_{\epsilon }^{1/2}\left(p^h\cdot n\right)\Vert }_{L^2\left({\Gamma }_N\right)}$ and ${\Vert h_{\epsilon }^{1/2}\left(q\cdot n\right)\Vert }_{L^2\left({\Gamma }_N\right)}$ vanished, we have the following conclusion

Lemma 3.6. Let $p$ and $p^h$ are the solutions of (2.2) and (2.3) and (3.3) and (3.4), $\text{div}{p}^h\in {\mathcal{L}}^0\left({\mathcal{T}}_h\right)$, then, we have,

$\begin{array}{c}\Vert p-p^h\Vert \le \underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }C\Vert p^h-q^h\Vert +\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }C{\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}{q}^h-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +CHh\Vert f^{\prime }\left(\hat{u}\right)\Vert \Vert f\left(u\right)\Vert +C{h}^2\Vert f\left(u_H\right)\Vert +C{\Vert h_{\epsilon }^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\mathrm{.}\end{array}$ (3.26)

The last three terms are high order terms compared with error $\Vert p^h-q^h\Vert$. By using the inverse inequality, we have the following upper bound result

Lemma 3.7. Suppose that the discrete $p^h$ satisfies $\text{curl}{p}^h=0$, $\text{div}{p}^h\in {\mathcal{L}}^0\left({\mathcal{T}}_h\right)$, then

$\Vert p-p^h\Vert \le C\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert p^h-q^h\Vert +\text{h}\text{.o}\text{.t}\mathrm{..}$ (3.27)

Combined with the efficiency (3.9) of the estimator, we have

Theorem 3.1. Let $p$ be the exact solution of (2.2) and (2.3) and $p^h\in V_h$ be the solution of (3.3) and (3.4), then, we have

${\eta }_z-\text{h}\text{.o}\text{.t}\text{.}\le \Vert p-p^h\Vert \le C{\eta }_z+\text{h}\text{.o}\text{.t}\mathrm{..}$ (3.28)

3.2. A Upper Bound for the Error of Two-Grid Algorithm 2

In this subsection, we will get the upper bounds for Algorithm 2. In order to make a theoretical analysis of Algorithm 2, we need to assume that the first derivative of $f\left(u\right)$ satisfies $f^{\prime }\left(u\right)<0$ in this subsection. First, we need the following priori error estimate for approximate solution ${\tilde{u}}^h$ from Algorithm 2.

Lemma 3.8. Let u be the exact solution of (2.2) and (2.3) and ${\tilde{u}}^h$ be the solution of (3.7) and (3.8), then, we have,

$\Vert u-{\tilde{u}}^h\Vert \le C{H}^2{\Vert u\Vert }_3^2\mathrm{.}$ (3.29)

Proof. Using (2.2) and (2.3) and (3.7) and (3.8), we have

$\left(p-{\tilde{p}}^h\mathrm{,}{v}_h\right)-\left(u-{\tilde{u}}^h\mathrm{,}\text{div}{v}_h\right)=\mathrm{0,}$

$\left(\text{div}\left(p-{\tilde{p}}^h\right)\mathrm{,}{w}_h\right)+\left(f^{\prime }\left(u_H\right){\tilde{u}}^h\mathrm{,}{w}_h\right)=\left(f\left(u\right)-f\left(u_H\right)\mathrm{,}{w}_h\right)+\left(f^{\prime }\left(u_H\right)u_H\mathrm{,}{w}_h\right)\mathrm{.}$

Using Taylor expansion for $f\left(u\right)$ at $u_H$, let $w_h=u-{\tilde{u}}^h$, $v_h=p-{\tilde{p}}^h$ and add the last two equations together to get

${\Vert p-{\tilde{p}}^h\Vert }^2-\left(f^{\prime }\left(u_H\right)\left(u-{\tilde{u}}^h\right)\mathrm{,}u-{\tilde{u}}^h\right)=\left(\frac{1}{2}{f}^{″}\left(\bar{u}\right){\left(u-u_H\right)}^2\mathrm{,}u-{\tilde{u}}^h\right)\mathrm{,}$ (3.30)

where $\bar{u}$ is some value between u and $u_H$. By using the assumption of $f\left(u\right)$ as well as the Cauchy inequality, we have the following estimation

${\Vert u-{\tilde{u}}^h\Vert }^2\le C{\Vert u-u_H\Vert }_{L^4\left(\Omega \right)}^2\Vert u-{\tilde{u}}^h\Vert \le C{H}^2{\Vert u\Vert }_3^2\Vert u-{\tilde{u}}^h\Vert \mathrm{,}$

which completes the proof.

Lemma 3.9. Let $p$ and ${\tilde{p}}^h$ are the solutions of (2.2) and (2.3) and (3.7) and (3.8). Then

$\begin{array}{c}{\displaystyle {\int }_{\Omega }}\left(p-{\tilde{p}}^h\right)\cdot \nabla \alpha \text{d}x\le C{H}^{2}h\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_3\Vert \nabla \alpha \Vert +C{H}^{2}h\Vert f^{″}\left(\bar{u}\right)\Vert {\Vert u\Vert }_3^2\Vert \nabla \alpha \Vert \\ +C{h}^2\Vert f\left(u_H\right)\Vert \Vert \nabla \alpha \Vert +CH{h}^2\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_2\Vert \nabla \alpha \Vert \\ +C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \mathrm{.}\end{array}$ (3.31)

Proof. Employ the Cauchy’s inequality and Poincaré inequality, we have

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-{\tilde{p}}^h\right)\cdot \nabla \alpha \text{d}x\\ ={\displaystyle {\int }_{\Omega }}\left(p-{\tilde{p}}^h\right)\cdot \nabla \left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}x\\ ={\displaystyle {\int }_{\Omega }}-\text{div}\left(p-{\tilde{p}}^h\right)\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}x+{\displaystyle {\int }_{{\Gamma }_N}}\left(p-{\tilde{p}}^h\right)\cdot n\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}s\\ ={\displaystyle {\int }_{\Omega }}\left(-f\left(u\right)+P_h\left(f\left(u_H\right)+f^{\prime }\left(u_H\right)\left({\tilde{u}}^h-u_H\right)\right)\right)\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}x\\ +{\displaystyle {\int }_{{\Gamma }_N}}\left(p-{\tilde{p}}^h\right)\cdot n\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}s\end{array}$

$\begin{array}{l}={\displaystyle {\int }_{\Omega }}\left(-f\left(u\right)+f\left(u_H\right)+f^{\prime }\left(u_H\right)\left({\tilde{u}}^h-u_H\right)\right)\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}x\\ +{\displaystyle {\int }_{\Omega }}(P_h\left(f\left(u_H\right)+f^{\prime }\left(u_H\right)\left({\tilde{u}}^h-u_H\right)\right)\\ -\left(f\left(u_H\right)+f^{\prime }\left(u_H\right)\left({\tilde{u}}^h-u_H\right)\right))\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}x\\ +{\displaystyle {\int }_{{\Gamma }_N}}\left(p-{\tilde{p}}^h\right)\cdot n\left(\alpha -{\alpha }_{\mathcal{T}}\right)\text{d}s\\ =Q_1+Q_2+Q_3\mathrm{,}\end{array}$ (3.32)

here ${\alpha }_{\mathcal{T}}\in {\mathcal{L}}^0\left({\mathcal{T}}_h\right)$ is the average of $\alpha$. Similar to $I_1-I_3$, we can estimate $Q_1-Q_3$ as

$\begin{array}{c}{Q}_1\le Ch\Vert f^{\prime }\left(u_H\right)\left(u-{\tilde{u}}^h\right)\Vert \Vert \nabla \alpha \Vert +Ch{\Vert \frac{1}{2}{f}^{″}\left(\bar{u}\right){\left(u-u_H\right)}^2\Vert }_{L^4\left(\Omega \right)}\Vert \nabla \alpha \Vert \\ \le C{H}^{2}h\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_3\Vert \nabla \alpha \Vert +C{H}^{2}h\Vert f^{″}\left(\bar{u}\right)\Vert {\Vert u\Vert }_3^2\Vert \nabla \alpha \Vert \mathrm{,}\end{array}$

$\begin{array}{c}{Q}_2\le C{h}^2\Vert \left(f\left(u_H\right)+f^{\prime }\left(u_H\right)\left(u-u_H\right)\right)\Vert \Vert \nabla \alpha \Vert \\ \le C{h}^{2}f\left(u_H\right)\Vert \nabla \alpha \Vert +CH{h}^2{f}^{\prime }\left(u_H\right){\Vert u\Vert }_2\Vert \nabla \alpha \Vert \mathrm{,}\end{array}$

$Q_3\le C{\Vert h^{1/2}\left(g_N-{\tilde{p}}^h\cdot n\right)\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \le C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\Vert \nabla \alpha \Vert \mathrm{.}$

Here $\bar{u}$ is some value between u and $u_H$, which completes the proof.

Similar to Lemma 3.4, we have the following result.

Lemma 3.10. Let $p$ and ${\tilde{p}}^h$ are solutions of (2.2) and (2.3) and (3.7) and (3.8), then, we have

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-{\tilde{p}}^h\right)\cdot \nabla w\text{d}x\\ \le C\Vert {\tilde{p}}^h-q\Vert +C{\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{div}\left(p-q\right)-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}+C{H}^{2}h\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_3+C{H}^{2}h\Vert f^{″}\left(\bar{u}\right)\Vert {\Vert u\Vert }_3^2\\ +C{h}^2\Vert f\left(u_H\right)\Vert +CH{h}^2\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_2\mathrm{.}\end{array}$ (3.33)

Lemma 3.11. Let $p$ and ${\tilde{p}}^h$ are the solutions of (2.2) and (2.3) and (3.7) and (3.8), then, we have

$\begin{array}{l}{\displaystyle {\int }_{\Omega }}\left(p-{\tilde{p}}^h\right)\cdot \text{curl}\beta \text{d}x\\ \le C\Vert \nabla \beta \Vert \Vert q-{\tilde{p}}^h\Vert +C\Vert \nabla \beta \Vert {\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}q-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +\Vert \nabla \beta \Vert \{C{\Vert h_{\epsilon }^{1/2}\left(q\cdot t\right)\Vert }_{L^2\left({\Gamma }_D\right)}+C{H}^{2}h\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_3+C{H}^{2}h\Vert f^{″}\left(\bar{u}\right)\Vert {\Vert u\Vert }_3^2\\ +C{h}^2\Vert f\left(u_H\right)\Vert +CH{h}^2\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_2\}\mathrm{.}\end{array}$ (3.34)

Combine with Lemma 3.9 and Lemma 3.11, and similar to Lemma 3.6, we finally get the following result.

Lemma 3.12. Let $p$ and ${\tilde{p}}^h$ are the exact and numerical solutions satisfying (2.2) and (2.3) and (3.7) and (3.8) respectively, then,

$\begin{array}{c}\Vert p-{\tilde{p}}^h\Vert \le \underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }C\Vert {\tilde{p}}^h-q^h\Vert +\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }C{\left({\displaystyle \underset{z\in \mathcal{K}}{\sum }}{h}_z^2\underset{f_z\in R}{\min }{\Vert \text{curl}{q}^h-f_z\Vert }_{L^2\left({\Omega }_z\right)}^2\right)}^{1/2}\\ +C{H}^{2}h\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_3+C{H}^{2}h\Vert f^{″}\left(\bar{u}\right)\Vert {\Vert u\Vert }_3^2+C{h}^2\Vert f\left(u_H\right)\Vert \\ +CH{h}^2\Vert f^{\prime }\left(u_H\right)\Vert {\Vert u\Vert }_2+C{\Vert h^{3/2}\partial g_N/\partial s\Vert }_{L^2\left({\Gamma }_N\right)}\mathrm{.}\end{array}$ (3.35)

By Lemma 3.12 and similar efficiency result (3.9), we have.

Theorem 3.2. Suppose ${\tilde{p}}^h$ satisfies $\text{curl}{\tilde{p}}^h=0$ and $\text{div}{\tilde{p}}^h\in {\mathcal{L}}^0\left({\mathcal{T}}_h\right)$, then,

$C\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert {\tilde{p}}^h-q^h\Vert +\text{h}\text{.o}\text{.t}\text{.}\le \Vert p-{\tilde{p}}^h\Vert \le C\underset{q^h\in S^1{\left({\mathcal{T}}_h\right)}^2}{\min }\Vert {\tilde{p}}^h-q^h\Vert +\text{h}\text{.o}\text{.t}\mathrm{..}$ (3.36)

4. Posteriori Error Estimator with Averaging Technique

In this section, we use the averaging technique to construct an averaging operator $\mathcal{A}\mathrm{<mo>\; :\; </mo>}{V}_h\to S^1{\left({\mathcal{T}}_h\right)}^2$ of posteriori error estimator, and prove $\Vert p^h-\mathcal{A}{p}^h\Vert$ (or $\Vert {\tilde{p}}^h-\mathcal{A}{\tilde{p}}^h\Vert$ ) is very approximation to $\Vert p-p^h\Vert$ (or $\Vert p-{\tilde{p}}^h\Vert$ ). In practise, we use the averaging operator $\mathcal{A}$ to compute the upper bound of the $\Vert p-p^h\Vert$ (or $\Vert p-{\tilde{p}}^h\Vert$ ). Take Algorithm 1 for instance (all the following conclusions are also hold for Algorithm 2). Set

${\eta }_{\mathcal{A}}=\Vert p^h-\mathcal{A}{p}^h\Vert \mathrm{,}$

then, the minimum ${\eta }_z$ is frequently replaced by an upper bound ${\eta }_{\mathcal{A}}$,

${\eta }_z\le {\eta }_{\mathcal{A}}\mathrm{.}$

From [18], one of a popular averaging operator $\mathcal{A}$ is defined, for each node z, by

$\mathcal{A}{p}^h\left(z\right)={\mathcal{A}}_z\left({p^h|}_{w_z}\right)\in R^2\mathrm{,}$

where ${\mathcal{A}}_z={\pi }_z\circ M_z$ with ${\pi }_z\mathrm{<mo>\; :\; </mo>}{R}^2\to R^2$ being an orthogonal projection and linear and continuous averaging $M_z$ being defined as

$M_z\left(p^h\right)={\displaystyle {\int }_{w_z}}{p}^h\text{d}x/\left|w_z\right|\mathrm{,}$

where $\left|w_z\right|$ denotes the area of patch $w_z$ related to node z.

For the efficiency of estimator ${\eta }_{\mathcal{A}}$, we have the following result.

Lemma 4.1 ( [18] ). There exists a mesh-size independent positive constant C with

${\eta }_{\mathcal{A}}\le C{\eta }_z\mathrm{.}$

According to the relationship between $\Vert p^h-q^h\Vert$ and $\Vert p^h-\mathcal{A}{p}^h\Vert$, we conclude the following result.

Theorem 4.1. Let $p$ and $p^h$ are the exact solution and numerical solutions from Algorithm 1 respectively, then, we have

$C{\eta }_{\mathcal{A}}-\text{h}\text{.o}\text{.t}\text{.}\le \Vert p-p^h\Vert \le C{\eta }_{\mathcal{A}}+\text{h}\text{.o}\text{.t}\mathrm{..}$ (4.1)

The result also holds for $\Vert p-{\tilde{p}}^h\Vert$, i.e.,

$C{\eta }_{\mathcal{A}}-\text{h}\text{.o}\text{.t}\text{.}\le \Vert p-{\tilde{p}}^h\Vert \le C{\eta }_{\mathcal{A}}+\text{h}\text{.o}\text{.t}\mathrm{..}$ (4.2)

5. Numerical Experiments

In this section, we will validate the a posteriori error estimates of averaging technique of two-grid mixed finite element solutions for semilinear elliptic equations by some numerical examples. Our focus is to deserve the ability of the error estimates to imitate the convergence behaviour of the error in the $L^2$ -norm.

To simplify our problem, we consider the following semilinear elliptic equations with entire boundary $\partial \Omega :={\Gamma }_D\left({\Gamma }_N=\varnothing \right)$ :

$\{\begin{array}{ll}\text{div}p=-u^3+g\left(x\right)\mathrm{,}\hfill & x\in \Omega \mathrm{,}\hfill \\ p=-\nabla u\mathrm{,}\hfill & x\in \Omega \mathrm{,}\hfill \\ u\left(x\right)=\mathrm{0,}\hfill & x\in \partial \Omega \mathrm{.}\hfill \end{array}$ (5.1)

And use the lowest order Raviart-Thomas mixed finite elements in the implementation. We consider the mesh $H\mathrm{,}h$ satisfy $h=H^2$ in the following numerical experiments.

Example 1: We choose $g\left(x\right)$ in a way such that the exact solution:

$u\left(x\right)=x_1\left(1-x_1\right)x_2\left(1-x_2\right),$

so we get the explicit expression of $g\left(x\right)$ as:

$g\left(x\right)=2\left(x_1\left(1-x_1\right)+x_2\left(1-x_2\right)\right)+{\left(x_1\left(1-x_1\right)x_2\left(1-x_2\right)\right)}^3\mathrm{,}$

with the domain $\Omega ={\left[\mathrm{0,1}\right]}^2$.

Example 2: We choose $g\left(x\right)$ in a way such that the exact solution:

$u\left(x\right)=\sin \left(\pi x_1\right)\sin \left(\pi x_2\right)\mathrm{,}$

so we get the explicit expression of $g\left(x\right)$ as:

$g\left(x\right)=2{\pi }^2\sin \left(\pi x_1\right)\sin \left(\pi x_2\right)+{\left(\sin \left(\pi x_1\right)\sin \left(\pi x_2\right)\right)}^3\mathrm{,}$

with the L-shape domain $\Omega =\left(-\mathrm{1,1}\right)\times \left(-\mathrm{1,1}\right)\mathrm{<mo>\; /\; </mo>}\left(\mathrm{0,1}\right)\times \left(\mathrm{0,1}\right)$.

From the numerical results presented in Table 1 and Table 2 for Example 1 and Table 3 and Table 4 for Example 2, we conclude that error estimators by averaging technique is reliable and efficient for both Algorithm 1 and Algorithm 2.

Table 1. A posteriori error of the Algorithm 1 for the example 1.

Table 2. A posteriori error of the Algorithm 2 for the example 1.

Table 3. A posteriori error of the Algorithm 1 for the example 2.

Table 4. A posteriori error of the Algorithm 2 for the example 2.

6. Conclusion

In this paper, we present a posteriori error estimate of two-grid mixed finite element method for semilinear elliptic equations by using averaging techniques. Theoretical analysis as well as the numerical experiments is provided to prove the efficiency and reliability of error estimators. In our following work, we will construct the adaptive two-grid mixed finite element method for the semilinear elliptic equations using the posteriori error estimators studied in this paper.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Durango, F. and Novo, J. (2019) A Posteriori Error Estimations for Mixed Finite Element Approximations to the Navier-Stokes Equations Based on Newton-Type Linearization. Journal of Computational and Applied Mathematics, 367, Article ID: 112429. https://doi.org/10.1016/j.cam.2019.112429
[2]	Larson, M.G. and Målqvist, A. (2008) A Posteriori Error Estimates for Mixed Finite Element Approximations of Elliptic Problems. Numerische Mathematik, 108, 487-500. https://doi.org/10.1007/s00211-007-0121-y
[3]	Carstensen, C. (1997) A Posteriori Error Estimate for the Mixed Finite Element Method. Mathematics of Computation, 66, 465-476. https://doi.org/10.1090/S0025-5718-97-00837-5
[4]	Chen, Y. and Liu, W. (2008) A Posteriori Error Estimates for Mixed Finite Element Solutions of Convex Optimal Control Problems. Journal of Computational and Applied Mathematics, 211, 76-89. https://doi.org/10.1016/j.cam.2006.11.015
[5]	Xu, J. (1994) A Novel Two-Grid Method for Semilinear Equations. SIAM Journal on Scientific Computing, 15, 231-237. https://doi.org/10.1137/0915016
[6]	Xu, J. (1996) Two-Grid Discretization Techniques for Linear and Non-Linear PDEs. SIAM Journal on Numerical Analysis, 33, 1759-1777. https://doi.org/10.1137/S0036142992232949
[7]	Chen, L. and Chen, Y. (2011) Two-Grid Method for Nonlinear Reaction Diffusion Equations by Mixed Finite Element Methods. Journal of Scientific Computing, 49, 383-401. https://doi.org/10.1007/s10915-011-9469-3
[8]	Chen, C. and Liu, W. (2012) A Two-Grid Method for Finite Element Solutions of Nonlinear Parabolic Equations. Abstract and Applied Analysis, 2012, Article ID: 391918. https://doi.org/10.1155/2012/391918
[9]	Sun, Y., Shi, F., Zheng, H., Li, H. and Wang, F. (2021) Two-Grid Domain Decomposition Methods for the Coupled Stokes-Darcy System. Computer Methods in Applied Mechanics and Engineering, 385, Article ID: 114041. https://doi.org/10.1016/j.cma.2021.114041
[10]	Chen, Y., Li, Q., Wang, Y. and Huang, Y. (2020) Two-Grid Methods of Finite Element Solutions for Semi-Linear Elliptic Interface Problems. Numerical Algorithms, 84, 307-330. https://doi.org/10.1007/s11075-019-00756-0
[11]	Bi, C., Wang, C. and Lin, Y. (2018) A Posteriori Error Estimates of Two-Grid Finite Element Methods for Nonlinear Elliptic Problems. Journal of Scientific Computing, 74, 23-48. https://doi.org/10.1007/s10915-017-0422-y
[12]	Chen, C., Chen, Y. and Zhao, X. (2018) A Posteriori Error Estimates of Two-Grid Finite Volume Element Methods for Nonlinear Elliptic Problems. Computers & Mathematics with Applications, 75, 1756-1766. https://doi.org/10.1016/j.camwa.2017.11.035
[13]	Li, Y. and Zhang, Y. (2021) Analysis of Adaptive Two-Grid Finite Element Algorithms for Linear and Nonlinear Problems. SIAM Journal on Scientific Computing, 43, 908-928. https://doi.org/10.1137/19M1285615
[14]	Carstensen, C. and Bartels, S. (2002) Each Averaging Technique Yields Reliable a Posteriori Error Control in FEM on Unstructured Grids. Part1: Low Order Conforming, Nonconforming, and Mixed FEM. Mathematics of Computation, 71, 945-969. https://doi.org/10.1090/S0025-5718-02-01402-3
[15]	Verfürth, R. (2013) A Posteriori Estimation Techniques for Finite Element Methods. Oxford University Press, Oxford. https://doi.org/10.1093/acprof:oso/9780199679423.001.0001
[16]	Hou, T., Chen, L. and Yang, Y. (2018) Two-Grid Methods for Expanded Mixed Finite Element Approximations of Semi-Linear Parabolic Integro-Differential Equations. Applied Numerical Mathematics, 132, 163-181. https://doi.org/10.1016/j.apnum.2018.06.001
[17]	Milner, F.A. (1985) Mixed Finite Element Methods for Quasilinear Second-Order Elliptic Problems. Mathematics of Computation, 44, 303-320. https://doi.org/10.1090/S0025-5718-1985-0777266-1
[18]	Carstensen, C. (2004) Some Remarks on the History and Future of Averaging Techniques in a Posteriori Finite Element Error Analysis. ZAMM—Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 84, 3-21. https://doi.org/10.1002/zamm.200410101