Nitsche's method for a Robin boundary value problem in a smooth domain

TL;DR

This paper establishes optimal error estimates for Nitsche's finite-element method applied to Robin boundary value problems in smooth domains, analyzing the influence of penalty parameters and confirming results with numerical examples.

Contribution

It provides new optimal error estimates for Nitsche's method applied to Robin problems, including the optimal relation between mesh size and penalty parameter, without requiring strong regularity assumptions.

Findings

Optimal error estimates are derived for the method.

The dependence of error on the penalty parameter is clarified.

Numerical examples confirm theoretical results.

Abstract

We prove several optimal-order error estimates for a finite-element method applied to an inhomogeneous Robin boundary value problem (BVP) for the Poisson equation defined in a smooth bounded domain in $\mathbb{R}^n$, $n=2,3$. The boundary condition is weakly imposed using Nitsche's method. The Robin BVP is interpreted as the classical penalty method with the penalty parameter $\varepsilon$. The optimal choice of the mesh size $h$ relative to $\varepsilon$ is a non-trivial issue. This paper carefully examines the dependence of $\varepsilon$ on error estimates. Our error estimates require no unessential regularity assumptions on the solution. Numerical examples are also reported to confirm our results.

Equations

$$-\Delta u$$

$$\frac{\partial u}{\partial \nu }+\frac{1}{\epsilon }u$$

$$u=u_0\text{mbox}on\Gamma$$

$$\Vert \tilde{u}-u_h^{\epsilon }{\Vert }_{H^1({\Omega }_h)}\le C[(h+{\epsilon }^{-\frac{1}{2}}{h}^2+{\epsilon }^{\frac{1}{2}})\Vert u{\Vert }_{H^4(\Omega )}+{\epsilon }^{-\frac{1}{2}}{h}^2\Vert {\tilde{u}}_0{\Vert }_{H^3({\Omega }_h)}].$$

$${\Vert \tilde{u}-u_h\Vert }_{h,\star }\le Ch\left[\left(1+\frac{h}{\epsilon }\right){\Vert u\Vert }_{2,\Omega }+{\Vert f\Vert }_{\Omega }+\frac{h}{\epsilon }{\Vert u_0\Vert }_{2,\Omega }+{\Vert g\Vert }_{1,\Omega }\right],$$

$${\left|v\right|}_{i,\mathcal{O}}={\left(\sum _{|\alpha |=i}{\Vert \frac{{\partial }^{\alpha }v}{\partial x^{\alpha }}\Vert }_{\mathcal{O}}^2\right)}^{1/2},{\Vert v\Vert }_{j,\mathcal{O}}={\left(\sum _{i=0}^j{\left|v\right|}_{i,\mathcal{O}}^2\right)}^{1/2}.$$

$$0\le \epsilon \le {\epsilon }_0$$

$${\Omega }_h=\mathrm{Int}\left(\bigcup _{K\in {\mathcal{T}}_h}K\right);$$

$$\mathcal{E}_{h}\mathrel{\mathop{:}}=\{E\subset\partial\Omega_{h}\colon\mbox{$E$ is an $(n-1)$-face of some $K\in\mathcal{T}_{h}$}\},$$

$$E\in {\mathcal{E}}_h$$

$$V_h:=\{\chi \in C(\overline{\Omega }):\chi {|}_K\in {\mathcal{P}}^1(K)(\forall K\in {\mathcal{T}}_h)\}.$$

$$\overline{\Omega \cup {\Omega }_h}\subset \tilde{\Omega }.$$

$$(Pv){|}_{\Omega }=v\text{mbox}in\Omega ,\Vert Pv{\Vert }_{W^{k,r}(\tilde{\Omega })}\le C_{l,r,\Omega }\Vert v{\Vert }_{W^{k,r}(\Omega )},$$

$$\tilde{f}:=Pf,{\tilde{u}}_0:=P{u}_0,\tilde{g}:=Pg.$$

$$a_h(w,v)$$

$$b_h(w,v)$$

$$+\frac{1}{\epsilon +\gamma h_E}⟨w,v{⟩}_E-\frac{\epsilon \gamma h_E}{\epsilon +\gamma h_E}⟨\frac{\partial w}{\partial {\nu }_h},\frac{\partial v}{\partial {\nu }_h}{⟩}_E],$$

$$\begin{gathered} l_h(v):=(\tilde{f},v{)}_{{\Omega }_h}+\sum _{E\in {\mathcal{E}}_h}[\{\frac{1}{\epsilon +\gamma h_E}⟨{\tilde{u}}_0,v{⟩}_E-\frac{\gamma h_E}{\epsilon +\gamma h_E}⟨{\tilde{u}}_0,\frac{\partial v}{\partial {\nu }_h}{⟩}_E \\ +\frac{\epsilon }{\epsilon +\gamma h_E}⟨\tilde{g},v{⟩}_E-\frac{\epsilon \gamma h_E}{\epsilon +\gamma h_E}⟨\tilde{g},\frac{\partial v}{\partial {\nu }_h}{⟩}_E], \end{gathered}$$

$$u_h\in V_h,a_h(u_h,\chi )=l_h(\chi )(\forall \chi \in V_h).$$

$${\Vert v\Vert }_h^2$$

$${\Vert v\Vert }_{h,\star }^2$$

$$\Vert v{\Vert }_h\le \Vert v{\Vert }_{h,\star }\le C\Vert v{\Vert }_h(v\in V_h).$$

$$\sum _{E\in {\mathcal{E}}_h}{h}_E{\Vert \frac{\partial v}{\partial {\nu }_h}\Vert }_E^2\le C\Vert \nabla v{\Vert }_K^2(v\in V_h,K\in {\mathcal{T}}_h,E\subset K).$$

$$\Vert v{\Vert }_E^2\le C{h}_E^{-1}\left(\Vert v{\Vert }_K^2+h_K^2\Vert \nabla v{\Vert }_K^2\right)(E\in {\mathcal{E}}_h,E\subset \partial K,v\in H^1(K)).$$

$$a_h(w,v)\le C{\Vert w\Vert }_{h,\star }{\Vert v\Vert }_{h,\star }(\forall w,v\in H^2({\Omega }_h)+V_h).$$

$$a_h(\chi ,\chi )\ge C{\Vert \chi \Vert }_h^2(\forall \chi \in V_h).$$

$$0<\gamma \le {\gamma }_0.$$

$$h\le C\underset{E\in {\mathcal{E}}_h}{\min }{h}_E.$$

$${\Vert \tilde{u}-u_h\Vert }_{h,\star }\le Ch\left[\left(1+\frac{h}{\epsilon }\right){\Vert u\Vert }_{2,\Omega }+{\Vert f\Vert }_{\Omega }+\frac{h}{\epsilon }{\Vert u_0\Vert }_{2,\Omega }+{\Vert g\Vert }_{1,\Omega }\right],$$

$${\Vert \tilde{u}-u_h\Vert }_{h,\star }\le Ch({\Vert u\Vert }_{2,\Omega }+{\Vert f\Vert }_{\Omega }+{\Vert u_0\Vert }_{2,\Omega }).$$

$${\Vert \tilde{u}-u_h\Vert }_{{\Omega }_h}$$

$$\displaystyle\hskip 50.00008pt+\left(\frac{h}{\varepsilon}+\frac{h^{2}}{\varepsilon^{2}}\right)\left\lVert u_{0}\right\rVert_{2,\Omega}+\left\lVert g\right\rVert_{3,\Omega}+\frac{h}{\varepsilon}\left\lVert g\right\rVert_{1,\Omega}\Biggr{]},$$

$${\Vert \tilde{u}-u_h\Vert }_{\Omega }\le C{h}^2({\Vert u\Vert }_{3,\Omega }+{\Vert u_0\Vert }_{2,\Omega }).$$

$$\mathrm{dist}(\Omega ,{\Omega }_h)\le C{h}^2$$

Full text

1. Introduction

Nitsche’s method [16] is well-known as a powerful method for imposing the Dirichlet boundary condition (DBC) in the finite element method (FEM). DBC is usually imposed by specifying the function values themselves at boundary nodal points. In contrast, Nitsche’s method is based on the method of “weak imposition” of DBC using penalty parameter. Actually, this strategy is useful for resolving the issue of spurious oscillations for non-stationary Navier–Stokes and convection–diffusion equations as was mentioned in Bazilevs et al. [6, 7].

In recent years, demand for computing complex boundary conditions has been increasing. Boundary conditions involving the Laplace–Beltrami operator, such as a dynamic boundary condition and a generalized Robin boundary condition play important roles in application to the reduced fluid–structure interaction model and Cahn–Hilliard equation (see, e.g., [10], [17] and [8]). Nitsche’s method may be an effective approach to address these boundary conditions, and therefore, is worthy of a thorough investigation.

When numerically solving PDEs in a smooth domain, we often utilize polyhedral approximations of the domain. Generally, a facile approximation of the problem may result in a wrong numerical solution; the so-called Babuška’s paradox in [2, §5] is a remarkable example. Therefore, investigating not only the error caused by discretizations but also that caused by domain approximations is important. For the standard FEM, approximating domains is a common problem, and analysis of the energy norm is well-developed thus far. Only recently, the optimal order $W^{1,\infty}$ and $L^{\infty}$ stability and error estimates were established; refer to [13] for detail.

Consequently, we evaluate Nitsche’s method for PDEs in a smooth bounded domain. We study a finite-element method (FEM) applied to an inhomogeneous Robin boundary value problem (BVP) for the Poisson equation defined in a smooth bounded domain $\Omega$ in $\mathbb{R}^{n}$, $n=2,3$. The boundary condition on the boundary $\Gamma=\partial\Omega$ is weakly imposed using Nitsche’s method. We then derive several optimal-order error estimates under reasonable regularity assumptions on the solution. Specifically, we consider

[TABLE]

instead of (1.1b).

The case of a polyhedral domain with $\varepsilon>0$ has already been addressed in Juntunen and Stenberg [11]. We are motivated by [11] and this paper is a generalization of [11] to a smooth domain. We study simultaneously the case $\varepsilon=0$, that is the case of DBC.

If we are concerned with the Dirichlet BVP (1.1a) and (1.1c), the Robin BVP (1.1a) and (1.1b) with $g=0$ implies the classical penalty method with the penalty parameter $\varepsilon\to 0$. (The $\varepsilon$ is interpreted as the penalty parameter in the classical penalty method. On the other hand, Nitsche’s method is introduced using the penalty parameter, which we will write as $\gamma$. The readers have to care not to confuse.) FEM for this method is well studied so far; we refer to [3, 4, 5, 15] for example. In particular, Barret and Elliott [4] presented the error estimate in a smooth domain as

[TABLE]

Therein, $h$ denotes the granularity parameter of the triangulation and $\Omega_{h}$ denotes a polyhedral approximation of $\Omega$ satisfying (2.2). (See also Remark 2.3.) The continuous $\mathcal{P}^{1}$ finite element solution is represented by $u_{h}^{\varepsilon}$. The definition of function spaces and their norms are described in the end of this section. Moreover, $\tilde{u}$ and $\tilde{u}_{0}$ are suitable extensions of $u$ and $u_{0}$, respectively. The precise definition will be mentioned in the next section. The estimate (1.2) gives the optimal-order estimate for the $H^{1}$ norm by setting $\varepsilon=Ch^{2}$. However, we need a surplus regularity $u\in H^{4}(\Omega)$. Barret and Elliott [5] later studied the iso-parametric FEM for a similar problem and obtained similar results as ours. However, regularity assumptions slightly vary from ours.

This paper carefully examines the dependence of $\varepsilon$ on error estimates. As a matter of fact, how to choose $h$ relative to $\varepsilon$ is a non-trivial issue for smooth domain cases. A suitable regularity assumption is another non-trivial issue as we recalled for the standard FEM above. This point was not discussed in Juntunen and Stenberg [11], because these issues do not appear in polyhedral cases. In fact, we succeed in deriving (see Theorem 1)

[TABLE]

where $\left\lVert\cdot\right\rVert_{h,\star}$ denotes the DG $H^{1}$ norm defined as (2.10). Consequently, we deduce the optimal-order estimate for the DG $H^{1}$ norm by setting $h=C\varepsilon$ under no further assumptions the smoothness of the solution and data. It makes this possible by applying some estimates reported in [13]. On the other hand, we assume a surplus regularity $u\in W^{3,q}(\Omega)$, $q>n$, for deducing the optimal-order estimate for the $L^{2}$ norm (see Theorem 2). In our opinion, this is an essential requirement; see Remark 2.2 and Section 6.

This paper comprises six sections. In Section 2, the continuous $\mathcal{P}^{1}$ finite element approximation using Nitsche’s method and the main error estimates (Theorems 1 and 2) are described. After having presented some preliminary results in Section 3, we prove Theorems 1 and 2 in Sections 4 and 5, respectively. Finally, numerical examples are also reported to confirm our results in Section 6.

Notation

We list the notations used in this paper. We follow the standard notation of, for example, [1] for function spaces and their norms. In particular, for $1\leq p\leq\infty$ and a positive integer $j$, we use the standard Lebesgue space $L^{p}(\mathcal{O})$ and Sobolev space $W^{j,p}(\mathcal{O})$. Hereinafter, $\mathcal{O}$ denotes a bounded domain in $\mathbb{R}^{n}$. The inner product and norm of $L^{2}(\mathcal{O})$ are denoted, respectively, by $(\cdot,\cdot)_{\mathcal{O}}$ and $\|\cdot\|_{\mathcal{O}}$. The norm of $W^{j,p}(\mathcal{O})$ is denoted by $\|\cdot\|_{W^{j,p}(\mathcal{O})}$. As usual, we set $H^{j}(\mathcal{O})=W^{j,2}(\mathcal{O})$, and the semi-norm and norm of $H^{j}(\mathcal{O})$ are denoted by, respectively,

[TABLE]

For $S\subset\partial\mathcal{O}$, we define $L^{p}(S)$ using a surface measure $dS$ in a common approach. The inner product and norm of $L^{2}(S)$ is denoted by, respectively, $\langle\cdot,\cdot\rangle_{S}$ and $\|\cdot\|_{S}$. Moreover, $\mathcal{P}^{r}(\mathcal{O})$ denotes the set of all polynomials in $\mathcal{O}$ of degree $\leq r$.

2. Nitsche’s method and the main results

We recall that $\Omega$ is a bounded domain in $\mathbb{R}^{n}$, $n=2,3$. Throughout this paper, we assume that the boundary $\Gamma$ of $\Omega$ is a $C^{l}$ boundary, where $l$ is an integer $\geq 2$.

From the general theory of elliptic PDEs, we know that the unique solution $u$ of (1.1) belongs to $H^{2}(\Omega)$ and satisfies $\|u\|_{2,\Omega}\leq C_{\varepsilon}(\|f\|_{\Omega}+\|u_{0}\|_{2,\Omega}+\|g\|_{1,\Omega})$, where $C_{\varepsilon}$ denotes a positive constant depending only on $\Omega$ and $\varepsilon$. The pure Neumann problem is out of our interest. Therefore, we assume

[TABLE]

for a suitably large $\varepsilon_{0}>0$.

Let $\{\mathcal{T}_{h}\}_{h}$ be a regular family of triangulations $\mathcal{T}_{h}$ of a polyhedral domain $\Omega_{h}\subset\mathbb{R}^{n}$ in the sense of [9]. That is,

1. (1)

$\mathcal{T}_{h}$ is a set of closed $n$-simplices (elements) $K$, and

[TABLE] 2. (2)

The granularity parameter $h$ is defined as $h=\displaystyle\max_{K\in\mathcal{T}_{h}}h_{K}$, where $h_{K}$ denotes the diameter of $K$; 3. (3)

Any two elements of $\mathcal{T}_{h}$ meet only in entire common faces or sides or in vertices; 4. (4)

There exists a positive constant $c_{\textup{reg}}$ satisfying $h_{K}\leq c_{\textup{reg}}\rho_{K}$ for all $K\in\mathcal{T}_{h}\in\{\mathcal{T}_{h}\}_{h}$, where $\rho_{K}$ denotes the diameter of the inscribed ball of $K$.

We then introduce the boundary mesh $\mathcal{E}_{h}$ inherited from $\mathcal{T}_{h}$ by

[TABLE]

and the boundary $\Gamma_{h}\mathrel{\mathop{:}}=\partial\Omega_{h}$ is expressed as $\Gamma_{h}=\bigcup_{E\in\mathcal{E}_{h}}E$. We assume that $\Gamma_{h}$ is an approximate surface/polygon of $\Gamma$ in the sense that

[TABLE]

We use the continuous $\mathcal{P}^{1}$ finite element space

[TABLE]

Below we fix a sufficiently smooth domain $\tilde{\Omega}\subset\mathbb{R}^{n}$ such that

[TABLE]

Since $\Gamma$ is of class $C^{l}$, $l\geq 2$, the domain $\Omega$ admits a strong $l$-extension operator $P$. That is, $P$ is a linear operator of $W^{l,r}(\Omega)\to W^{l,r}(\tilde{\Omega})$ for any $0\leq k\leq l$ and $1\leq r<\infty$, and it satisfies

[TABLE]

where $C_{l,r,\Omega}$ denotes a positive constant depending only on $l$, $r$ and $\Omega$; see [1, Theorem 5.22] for example. Using this, we write

[TABLE]

Following [11], we set

[TABLE]

where $\gamma>0$ is a penalty parameter, $h_{E}$ the diameter of $E$, and $\nu_{h}$ the outer unit normal vector to $\Gamma_{h}$.

Nitsche’s method for (1.1) is stated as follows:

[TABLE]

We use the following norms that depend on $\varepsilon$ and $h_{E}$:

[TABLE]

We recall that $\left\lVert\cdot\right\rVert_{h}$ and $\left\lVert\cdot\right\rVert_{h,\star}$ are equivalent on $V_{h}$ uniformly in $h$. That is, there exists a positive constant $C_{0}$ independent of $h$ such that

[TABLE]

Here and hereinafter, $C$ denotes a generic positive constant which is independent of $h$ and $\varepsilon$. The value of $C$ may be different at each occurrence. The inequalities (2.11) follow from the well-known inequality

[TABLE]

In fact, (2.12) is a readily obtainable consequence of the standard trace inequality,

[TABLE]

Nitsche’s method (2.8) admits a unique solution in view of the following basic result; see [11, Theorem 3.2].

Lemma 2.1**.**

We have

[TABLE]

Actually, $\gamma_{0}$ can be taken as any positive number strictly smaller than $1/C$, where $C$ denotes the constant appearing in (2.12); see [11, Theorem 3.2]. Below we always assume that

[TABLE]

To deduce convergence results, we need an inverse assumption as

[TABLE]

We are now in a position to state our main result. We recall that $C$ denotes a positive constant which is independent of $\varepsilon$ and $h$.

Theorem 1** ($H^{1}$ estimates).**

Suppose that $\Gamma$ is a $C^{2}$ boundary. Let $u\in H^{2}(\Omega)$ and $u_{h}\in V_{h}$ represent the solutions of (1.1) and (2.8), respectively. Assume that (2.2), (2.15) and (2.16) are satisfied. Then, if $\varepsilon>0$, we have

[TABLE]

where $\tilde{u}=Pu$. If $\varepsilon=0$, we have

[TABLE]

Theorem 2** ($L^{2}$ estimates).**

Suppose that $\Gamma$ is a $C^{3}$ boundary. Let $u\in H^{2}(\Omega)$ and $u_{h}\in V_{h}$ represent the solutions of (1.1) and (2.8), respectively. Assume that (2.2), (2.15) and (2.16) are satisfied. Then, if $\varepsilon>0$, $u\in W^{3,q}(\Omega)$ for some $q>n$, $u_{0}\in H^{3}(\Omega)$ and $g\in H^{3}(\Omega)$, we have

[TABLE]

where $\tilde{u}=Pu$. On the other hand, if $\varepsilon=0$ and $u\in H^{3}(\Omega)$, we have

[TABLE]

Remark 2.2*.*

Theorem 1 reports that the optimal rate of convergence for the $H^{1}$ error is achieved under a reasonable (minimal) regularity assumption on $u$. On the other hand, we pose a somewhat surplus regularity $u\in W^{3,q}(\Omega)$, $q>n$, for deducing the optimal rate of convergence for the $L^{2}$ error. In our opinion, this is an essential requirement. Actually, a numerical example reported in Section 6 shows the second-order convergence may not take place if $u\in W^{3,q}(\Omega)$, $q>n$.

Remark 2.3*.*

We are assuming (2.2) for $\Omega_{h}$. This can be replaced by

[TABLE]

with some obvious modification of proofs.

3. Boundary-skin estimates

We collect some auxiliary results that will be used in the proof of the main results.

Since $\Omega$ is a bounded $C^{l}$ domain, there exists a local coordinate system $\{U_{r},y_{r},\phi_{r}\}_{r=1}^{M}$ to ensure the following:

$\{U_{r}\}_{r=1}^{M}$ is an open covering of $\Gamma$. 2. 2)

For any $1\leq r\leq M$, there exists a congruent transformation $A_{r}$ such that $y_{r}=(y_{r1},y^{\prime}_{r})=A_{r}(x)$, where $x$ is the original coordinate. 3. 3)

For any $1\leq r\leq M$, $\phi_{r}$ is a $C^{l}$ function in $\Delta_{r}\mathrel{\mathop{:}}=\{y_{r1}\in\mathbb{R}^{n-1}\colon\left\lvert y_{r1}\right\rvert\leq\alpha\}$ and $\Gamma\cap U_{r}$ is a graph of $\phi_{r}$ with respect to the coordinate $y_{r}$.

Assuming that $h$ is sufficiently small if necessary, our possible assumptions are as follows:

For any $1\leq r\leq M$, there exists a function $\phi_{rh}$ such that $\Gamma_{h}\cap U_{r}$ is a graph of $\phi_{rh}$ with respect to the coordinate $y_{r}$.

In addition, we assume that $h_{0}$ is sufficiently small to ensure that for any $x\in\Gamma$ and $r=1,\ldots,M$, the open ball $B(x,h_{0})$ with center $x$ and radius $h_{0}$ is contained in a neighborhood $U_{r}$. Let $d(x)$ be the signed distance function defined by

[TABLE]

We define $\Gamma(\delta)\mathrel{\mathop{:}}=\{x\in\mathbb{R}^{n}\colon\left\lvert d(x)\right\rvert<\delta\}$, which we call the boundary-skin region. Then, for a sufficiently small $\delta$, the orthogonal projection $\pi$ onto $\Gamma$ exists such that

[TABLE]

Because $h$ is sufficiently small, $\pi$ is defined on $\Gamma_{h}\subset\Gamma(\delta)$ and for each $E\in\mathcal{E}_{h}$, and $\pi(E)$ comprises some local neighborhood $U_{r}$. In this case, $\pi|_{\Gamma_{h}}$ has the inverse operator $\pi^{*}(x)=x+t^{*}(x)\nu(x)$, and $\left\lVert t^{*}\right\rVert_{L^{\infty}(\Gamma)}\leq Ch^{2}$. Moreover, $\pi(\mathcal{E}_{h})\mathrel{\mathop{:}}=\{\pi(E)\colon E\in\mathcal{E}_{h}\}$ is a partition of $\Gamma$.

We assume that all these properties hold for any $h\leq h_{0}$ by assuming that $h_{0}$ is sufficiently small if necessary.

Now we can state the boundary-skin estimates. For the proof, refer to [14, Theorems 8.1, 8.2, and 8.3 and Lemma 9.1] and [13, Lemma A.1].

Lemma 3.1** (Boundary-skin estimates).**

Let $\delta\in[c^{-1}h^{2},~{}ch^{2}]$ with a positive constant $c$. We have

[TABLE]

Lemma 3.2**.**

[TABLE]

• Proof.

It is a direct consequence of (3.2d) in view of (2.1) and (2.15). ∎

Lemma 3.3**.**

[TABLE]

• Proof.

First, consider the case $\varepsilon=0$. Let $\delta$ be a constant satisfying $c^{-1}h^{2}\leq\delta\leq ch^{2}$ as in Lemma 3.1. Since $u\circ\pi=u\circ u_{0}$ on $\Gamma_{h}$, we have by (3.2b)

[TABLE]

which implies (3.4a). Therein, we have used $\|\tilde{u}\|_{1,\Gamma(\delta)}\leq Ch\|\tilde{u}\|_{2,\Omega}$ and $\|\tilde{u}_{0}\|_{1,\Gamma(\delta)}\leq Ch\|\tilde{u}_{0}\|_{2,\Omega}$, which are direct consequences of (3.2c) and the trace theorem.

We proceed to the case $\varepsilon>0$. Since $(\nabla\tilde{u}\circ\pi)\cdot(\nu\circ\pi)+\varepsilon^{-1}\tilde{u}\circ\pi=\tilde{g}\circ\pi+\varepsilon^{-1}\tilde{u}_{0}\circ\pi$ on $\Gamma_{h}$, we have

[TABLE]

As above, we estimate as

[TABLE]

On the other hand, by (3.2b) and (3.2a)

[TABLE]

Summing up, we deduce (3.4b). ∎

4. Proof of Theorem 1

We start with a version of the Strang lemma. To state it, we set

[TABLE]

for $v\in H^{2}(\tilde{\Omega})+V_{h}$ and $\chi\in V_{h}$.

Lemma 4.1**.**

Under the same assumption of Theorem 1, we have

[TABLE]

• Proof.

Letting $\xi\in V_{h}$ and $\chi=u_{h}-\xi$, we have

[TABLE]

where (2.14b), (2.14a), and (2.11) are applied. This, together with the triangular inequality, implies (4.2) ∎

The standard Lagrange interpolation operator of $C(\overline{\Omega})\to V_{h}$ is denoted by $\Pi_{h}$. It is well-known that

[TABLE]

As a direct application of (4.3) and (2.13), we derive

[TABLE]

We can state the following proof.

• Proof of Theorem 1.

First consider the case $\varepsilon>0$. In view of (4.2) and (4.4), it suffices to prove

[TABLE]

Let $\chi\in V_{h}$ be arbitrary. Applying the integrate by parts, we have

[TABLE]

Since $-\Delta\tilde{u}=\tilde{f}$ in $\Omega$, we have by (3.3)

[TABLE]

Using (3.4b),

[TABLE]

Summing up, we deduce (4.5).

We proceed to the case $\varepsilon=0$. It suffices to prove

[TABLE]

In this case, we have $r_{h}(\tilde{u},\chi)=\textup{I}_{1}+\textup{I}_{3}$, where

[TABLE]

We apply (2.16) and (3.4a) to obtain

[TABLE]

where $h_{\min}=\min_{E}h_{E}$. Therefore, (4.8) is proved. ∎

5. Proof of Theorem 2

We use the Green operator $G:L^{2}(\Omega_{h})\to H^{2}(\Omega)$ defined as follows. For $\varepsilon>0$, $z=G\eta\in H^{2}(\Omega)$ is the unique solution of

[TABLE]

For any $\eta\in L^{2}(\Omega_{h})$, $G\eta$ admits the a priori estimate

[TABLE]

In particular, $C$ does not depend on $\varepsilon$. We omit the proof since it is outside the scope of this paper. As a matter of fact, this can be verified by a standard method of difference quotient. For example, if tracing the proof of [12, Theorem 3.3] carefully, we find that the estimate (5.2) holds true. Moreover, if $\Gamma$ is a $C^{4}$ boundary, the same proof of [18, Lemma 4.1] is also applicable.

Lemma 5.1**.**

Under the same assumption of Theorem 2, we have

[TABLE]

where $\tilde{z}=Pz\in H^{2}(\tilde{\Omega})$, $z=G{\eta}\in H^{2}(\Omega)$ and

[TABLE]

• Proof.

First, suppose that $\varepsilon>0$. In the similar way as the derivation of (4.6), we deduce

[TABLE]

for $\chi\in H^{2}(\Omega_{h})+V_{h}$ and $v\in H^{2}(\Omega_{h})$.

Let $\eta\in L^{2}(\Omega_{h})$ and $\tilde{z}=Pz=P(G\eta)$. We use the same symbol $\eta$ to express the zero extension of $\eta$ into $\Omega\backslash\Omega_{h}$. Substituting (5.4) for $\chi=\tilde{u}-u_{h}$ and $v=\tilde{z}$, we obtain

[TABLE]

Thanks to (5.2), an estimation for $\textup{J}_{3}$ is readily;

[TABLE]

We apply Lemma 3.3 with $u_{0}=0$ and $g=0$ and derive

[TABLE]

Summing up, we deduce

[TABLE]

Finally, using (2.14a) and (5.2), we deduce (5.3) for $\varepsilon>0$.

If $\varepsilon=0$, $\textup{J}_{4}$ is replaced by

[TABLE]

We apply Lemma 3.3 with $u_{0}=0$ and get

[TABLE]

Therefore, (5.3) holds even for $\varepsilon=0$. ∎

We finally state the following proof.

• Proof of Theorem 2.

We define following bilinear and linear forms:

[TABLE]

Then, we obtain

[TABLE]

First consider the case $\varepsilon>0$. Since $\tilde{u}$ and $\tilde{z}$ are continuous in $\Omega_{h}$, we have

[TABLE]

Using (3.2c), we obtain

[TABLE]

Given that

[TABLE]

we have

[TABLE]

Using (3.2a), we obtain

[TABLE]

Using (5.1b) and (3.2b), we have

[TABLE]

Hence, we obtain

[TABLE]

We rearrange $\textup{A}_{2}$ as

[TABLE]

Using boundary skin estimates, we get

[TABLE]

Therein, we used the Sobolev inequality and (2.5) to derive

[TABLE]

(To apply (2.5) we are assuming that $\Gamma$ is a $C^{3}$ boundary.)

Therefore, we deduce

[TABLE]

Similarly, we have

[TABLE]

$\textup{A}_{3}$ is rewritten as

[TABLE]

Using boundary skin estimates and (5.1b), we can perform estimations as

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

So, we get

[TABLE]

Summing up, we deduce

[TABLE]

We apply (3.4b), we obtain

[TABLE]

Consequently,

[TABLE]

Using the same way as the proof of (4.5), we derive

[TABLE]

for $w\in H^{2}(\tilde{\Omega})+V_{h}$. By substituting this for $w=\tilde{z}-\Pi_{h}\tilde{z}$, we have

[TABLE]

Finally, we have

[TABLE]

and we obtain the estimate (2.19).

We proceed to the case $\varepsilon=0$. We replace $\textup{I}_{2}$ by

[TABLE]

and $\textup{I}_{3}$ by

[TABLE]

Using (3.2a) and (3.2b), we have

[TABLE]

Using (3.4a) and (3.2b), we obtain

[TABLE]

Since $\textup{A}_{2}=\textup{A}_{21}$, we get

[TABLE]

Using (3.4a), we have

[TABLE]

Hence, we obtain

[TABLE]

Therefore, (2.20) is proved. ∎

6. Numerical examples

In this section, we present some numerical results to verify the validity of our error estimates. We consider the Poisson problem (1.1) in a disk $\Omega=\{x\in\mathbb{R}^{2}\colon\left\lvert x\right\rvert<1\}$.

First, we confirm the validity of the estimates in Theorems 1 and 2. We set $f$, $g$ and $u_{0}$ so that a function $u(x_{1},x_{2})=\sin(x_{1})\sin(x_{2})$ solves (1.1). Let $u_{h}\in V_{h}$ be the solution of (2.8). Figure 1 shows the the $H^{1}$ error $\left\lVert u-u_{h}\right\rVert_{h,\star}$ and the $L^{2}$ error $\left\lVert u-u_{h}\right\rVert_{\Omega_{h}}$ for $\varepsilon=0,1$. We observe that the convergence rates are almost $O(h)$ for the $H^{1}$ error and $O(h^{2})$ for the $L^{2}$ error. Thus, the optimal convergence rates actually take place and the estimates of Theorems 1 and 2 are confirmed.

Subsequently, we consider the exact solution $u(x_{1},x_{2})=((x_{1}+1)^{2}+x_{2}^{2})^{\frac{2}{3}}$ and the corresponding $f$, $g$, and $u_{0}$. Let $\varepsilon=1$. In this case, $u\in H^{2}(\Omega)$ and $u\not\in W^{3,q}(\Omega)$ for any $q>n$. That is, the assumption of Theorem 2 does not satisfied. Figure 2 reports the $H^{1}$ error and the $L^{2}$ error for $\varepsilon=1$. We see from Figure 2 that the convergence rate for the $H^{1}$ error is $O(h)$. However, the second-order convergence does not achieve for the $L^{2}$ error. Actually, we observe that the convergence rate for the $L^{2}$ error is $O(h^{2-\sigma})$ with some small $\sigma>0$. This result is consistent with Theorem 2. Therefore, we can conclude that the regularity condition $u\in W^{3,q}(\Omega)$ with $q>n$ is an essential requirement for deducing the optimal order convergence.

Acknowledgements

This work was supported by CREST (JPMJCR15D1) of JST, Japan, by Grant-in-Aid for Scientific Research B (15H03635) of JSPS, Japan, and by Grant-in-Aid for Scientific Research A (21H04431) of JSPS, Japan.