Uniform error analysis of a rectangular Morley finite element method on a Shishkin mesh for a 4th-order singularly perturbed boundary value problem

TL;DR

This paper analyzes the error of a rectangular Morley finite element method on a Shishkin mesh for a 4th-order singularly perturbed boundary value problem, showing optimal convergence rates in challenging regimes.

Contribution

It provides a uniform error analysis and sharp convergence estimates for the Morley finite element method on Shishkin meshes for a complex singularly perturbed problem.

Findings

Achieves an $O(N^{-3/2})$ convergence rate in the most difficult regime.

Demonstrates the method's superiority over Adini finite elements in similar settings.

Validates theoretical results with numerical experiments.

Abstract

The singularly perturbed reaction-diffusion problem $\varepsilon^2\Delta^2 u - \mathrm{div}\left(c\nabla u\right) = f$ is considered on the unit square $\Omega$ in $\mathbb{R}^2$ with homogenous Dirichlet boundary conditions. Its solution typically contains boundary layers on all sides of~$\Omega$. It is discretised by a finite element method that uses rectangular Morley elements on a Shishkin mesh. In an associated energy-type norm that is natural for this problem, we prove an $O(\varepsilon^{1/2}N^{-1}+\varepsilon N^{-1}\ln N + N^{-3/2})$ rate of convergence for the error in the computed solution, where $N$~is the number of mesh intervals in each coordinate direction. Thus in the most troublesome regime when $\varepsilon \approx N^{-1}$, our method is proved to attain an $O(N^{-3/2})$ rate of convergence, which is shown to be sharp by our numerical experiments and is superior to the $O(N^{-1/2})$ rate that is proved in Meng & Stynes, Adv. Comput. Math. 2019 when Adini finite elements are used to solve the same problem on the same mesh.

Tables

Table 1. Table 1: Numerical results for Example  5.1

$\epsilon$	$N$ = 16	$N$ = 32	$N$ = 64	$N$ = 128
1.0e-00	2.08e-03	1.02e-03	5.07e-04	2.52e-04
	1.02	1.01	1.00
1.0e-01	2.46e-02	1.12e-02	5.33e-03	2.59e-03
	1.13	1.07	1.03
1.0e-02	3.79e-02	1.19e-02	4.46e-03	1.90e-03
	1.66	1.42	1.23
1.0e-03	3.52e-02	1.01e-02	3.23e-03	1.14e-03
	1.79	1.65	1.50
1.0e-04	3.14e-02	8.44e-03	2.34e-03	7.05e-04
	1.89	1.84	1.73
1.0e-05	2.99e-02	7.77e-03	2.02e-03	5.43e-04
	1.94	1.94	1.89
1.0e-06	2.94e-02	7.54e-03	1.91e-03	4.92e-04
	1.96	1.97	1.96
1.0e-07	2.92e-02	7.47e-03	1.88e-03	4.75e-04
	1.96	1.98	1.98
1.0e-08	2.92e-02	7.45e-03	1.87e-03	4.70e-04
	1.97	1.99	1.99

Table 2. Table 2: Numerical results for Example  5.1 with ε = N − 1 \varepsilon=N^{-1}

$N$ = 16	$N$ = 32	$N$ = 64	$N$ = 128	$N$ = 256
2.66e-02	1.18e-02	4.76e-03	1.79e-03	6.50e-04
1.16	1.31	1.41	1.46

Table 3. Table 3: Numerical results for Example  5.2

$\epsilon$	$N$ = 16	$N$ = 32	$N$ = 64	$N$ = 128
1.0e-00	1.45e-02	7.15e-03	3.54e-03	1.76e-03
	1.02	1.01	1.00
1.0e-01	6.17e-02	3.13e-02	1.55e-02	7.67e-03
	0.97	1.01	1.01
1.0e-02	3.94e-02	1.71e-02	8.66e-03	4.75e-03
	1.20	0.98	0.86
1.0e-03	3.05e-02	1.01e-02	3.99e-03	1.81e-03
	1.58	1.35	1.13
1.0e-04	2.69e-02	7.61e-03	2.36e-03	8.53e-04
	1.82	1.68	1.47
1.0e-05	2.57e-02	6.78e-03	1.84e-03	5.45e-04
	1.92	1.87	1.76
1.0e-06	2.53e-02	6.52e-03	1.68e-03	4.47e-04
	1.95	1.95	1.91
1.0e-07	2.52e-02	6.44e-03	1.63e-03	4.16e-04
	1.96	1.98	1.96
1.0e-08	2.51e-02	6.41e-03	1.61e-03	4.07e-04
	1.97	1.98	1.98

Table 4. Table 4: Numerical results for Example  5.3

$\epsilon$	$N$ = 16	$N$ = 32	$N$ = 64	$N$ = 128
1.0e-00	1.22e-02	6.01e-03	2.98e-03	1.48e-03
	1.02	1.01	1.00
1.0e-01	1.00e-01	4.65e-02	2.21e-02	1.07e-02
	1.11	1.07	1.03
1.0e-02	1.59e-01	4.48e-02	1.69e-02	7.95e-03
	1.82	1.40	1.08
1.0e-03	1.82e-01	5.08e-02	1.50e-02	4.96e-03
	1.84	1.76	1.59
1.0e-04	1.79-01	4.77e-02	1.27e-02	3.59e-03
	1.90	1.90	1.82
1.0e-05	1.77e-01	4.65e-02	1.20e-02	3.13e-03
	1.93	1.95	1.93
1.0e-06	1.76e-01	4.61e-02	1.17e-02	2.99e-03
	1.93	1.97	1.97
1.0e-07	1.76e-01	4.60e-02	1.16e-02	2.94e-03
	1.93	1.97	1.98
1.0e-08	1.76e-01	4.60e-02	1.16e-02	2.93e-03
	1.94	1.97	1.99

Equations

$${\epsilon }^2{\Delta }^{2}u(x,y)-\mathrm{div}\left(c(x,y)\nabla u(x,y)\right)$$

$$u(x,y)=\frac{\partial u(x,y)}{\partial n}$$

$$u(x,y)$$

$$:=S(x,y)+u^E(x,y)$$

$$\left|\frac{{\partial }^{i+j}S(x,y)}{\partial x^i\partial y^j}\right|$$

$$\left|\frac{{\partial }^{i+j}{E}_1(x,y)}{\partial x^i\partial y^j}\right|$$

$$\left|\frac{{\partial }^{i+j}{E}_{12}(x,y)}{\partial x^i\partial y^j}\right|$$

$$p_1(x,y)=\frac{1}{8}\left[2\left(1-\frac{x-x_c}{h_x}\right)\left(1-\frac{y-y_c}{h_y}\right)+{\left(\frac{x-x_c}{h_x}\right)}^3-\frac{x-x_c}{h_x}+{\left(\frac{y-y_c}{h_y}\right)}^3-\frac{y-y_c}{h_y}\right],$$

$$p_2(x,y)=\frac{1}{8}\left[2\left(1-\frac{x-x_c}{h_x}\right)\left(1+\frac{y-y_c}{h_y}\right)+{\left(\frac{x-x_c}{h_x}\right)}^3-\frac{x-x_c}{h_x}-{\left(\frac{y-y_c}{h_y}\right)}^3+\frac{y-y_c}{h_y}\right],$$

$$p_3(x,y)=\frac{1}{8}\left[2\left(1+\frac{x-x_c}{h_x}\right)\left(1+\frac{y-y_c}{h_y}\right)-{\left(\frac{x-x_c}{h_x}\right)}^3+\frac{x-x_c}{h_x}-{\left(\frac{y-y_c}{h_y}\right)}^3+\frac{y-y_c}{h_y}\right],$$

$$p_4(x,y)=\frac{1}{8}\left[2\left(1+\frac{x-x_c}{h_x}\right)\left(1-\frac{y-y_c}{h_y}\right)-{\left(\frac{x-x_c}{h_x}\right)}^3+\frac{x-x_c}{h_x}+{\left(\frac{y-y_c}{h_y}\right)}^3-\frac{y-y_c}{h_y}\right],$$

$$q_1(x,y)=-\frac{h_y}{4}\left(\frac{y-y_c}{h_y}+1\right){\left(\frac{y-y_c}{h_y}-1\right)}^2,q_2(x,y)=-\frac{h_x}{4}\left(\frac{x-x_c}{h_x}+1\right){\left(\frac{x-x_c}{h_x}-1\right)}^2,$$

$$q_3(x,y)=\frac{h_y}{4}{\left(\frac{y-y_c}{h_y}+1\right)}^2\left(\frac{y-y_c}{h_y}-1\right),q_4(x,y)=\frac{h_x}{4}{\left(\frac{x-x_c}{h_x}+1\right)}^2\left(\frac{x-x_c}{h_x}-1\right).$$

$$\{\begin{array}{ll}{p}_i(a_j)={\delta }_{ij}\text{ and }\frac{1}{|e_j|}{\int }_{e_j}\frac{\partial p_i}{\partial n_{e_j}}\mathrm{d}s=0& \text{ for }i,j=1,2,3,4,\\ q_i(a_j)=0\text{ and }\frac{1}{|e_j|}{\int }_{e_j}\frac{\partial q_i}{\partial n_{e_j}}\mathrm{d}s={\delta }_{ij}& \text{for }i,j=1,2,3,4,\end{array}$$

$${\epsilon }^2(\Delta u,\Delta v)+\left(c\nabla u,\nabla v\right)=(f,v)\text{ for all }v\in H_0^2(\Omega ),$$

$${\epsilon }^2{a}_N(u_N,v_N)+b_N(u_N,v_N)=(f,v_N{)}_{L^2(\Omega )}\forall v_N\in V_N.$$

$$(v,w_N{)}_K=(P_K^{0}v,w_N{)}_K\forall v\in L^2(K),\forall w_N\in P_0(K).$$

$${\Vert v-P_K^{0}v\Vert }_{L^2(K)}^2\le C\left(h_x^2{\Vert \frac{\partial v}{\partial x}\Vert }_{L^2(K)}^2+h_y^2{\Vert \frac{\partial v}{\partial y}\Vert }_{L^2(K)}^2\right),$$

$${\Vert v-I_{N}v\Vert }_{L^2(K)}^2\le C\left(h_x^2{\Vert \frac{\partial v}{\partial x}\Vert }_{L^2(K)}^2+h_y^2{\Vert \frac{\partial v}{\partial y}\Vert }_{L^2(K)}^2\right)$$

$${\Vert I_{N}v\Vert }_{H^1(K)}\le C{\Vert v\Vert }_{H^1(K)}.$$

$$h_x{\Vert \frac{{\partial }^2{\xi }_N}{\partial x^2}\Vert }_{L^2(K)}+h_y{\Vert \frac{{\partial }^2{\xi }_N}{\partial x\partial y}\Vert }_{L^2(K)}$$

$$h_y{\Vert \frac{{\partial }^2{\xi }_N}{\partial y^2}\Vert }_{L^2(K)}+h_x{\Vert \frac{{\partial }^2{\xi }_N}{\partial x\partial y}\Vert }_{L^2(K)}$$

$${\epsilon }^2{a}_N({\xi }_N,{\xi }_N)$$

$$={\epsilon }^2{a}_N({\eta }_N,{\xi }_N)+{\epsilon }^2{a}_N(u_N-u,{\xi }_N)+b_N({\eta }_N,{\xi }_N)+b_N(u_N-u,{\xi }_N)$$

$$={\epsilon }^2{a}_N({\eta }_N,{\xi }_N)+b_N({\eta }_N,{\xi }_N)+({\epsilon }^2{\Delta }^{2}u-\mathrm{div}(c\nabla u),I_N{\xi }_N)$$

$$+(f,{\xi }_N-I_N{\xi }_N)-{\epsilon }^2{a}_N(u,{\xi }_N)-b_N(u,{\xi }_N).$$

$$({\Delta }^{2}u,I_N{\xi }_N)=-(\nabla \Delta u,\nabla (I_N{\xi }_N))\text{ and }-(\mathrm{div}(c\nabla u),I_N{\xi }_N)-b_N(u,{\xi }_N)=b_N(u,I_N{\xi }_N-{\xi }_N).$$

$$a_N(u,{\xi }_N)=\sum _{K\in {\mathcal{T}}_N}(\Delta u,\Delta {\xi }_N{)}_{L^2(K)}$$

$$=-\sum _{K\in {\mathcal{T}}_N}(\nabla \Delta u,\nabla {\xi }_N{)}_{L^2(K)}+\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial n^2}\frac{\partial {\xi }_N}{\partial n}\mathrm{d}s+\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial s\partial n}\frac{\partial {\xi }_N}{\partial s}\mathrm{d}s,$$

$${\epsilon }^2{a}_N({\xi }_N,{\xi }_N)$$

$$={\epsilon }^2{a}_N({\eta }_N,{\xi }_N)+b_N({\eta }_N,{\xi }_N)-b_N(u,{\xi }_N-I_N{\xi }_N)+(f,{\xi }_N-I_N{\xi }_N)$$

$$-{\epsilon }^2\left[\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial n^2}\frac{\partial {\xi }_N}{\partial n}\mathrm{d}s+\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}u}{\partial s\partial n}\frac{\partial {\xi }_N}{\partial s}\mathrm{d}s\right]$$

$$+{\epsilon }^2\sum _{K\in {\mathcal{T}}_N}{\int }_K(\nabla \Delta u)\nabla ({\xi }_N-I_N{\xi }_N)\mathrm{d}x\mathrm{d}y.$$

$${\Pi }_{e_i}^{0}w=\frac{1}{|e_i|}{\int }_{e_i}w\mathrm{d}s\text{ and }{\mathcal{R}}_{e_i}^{0}w=w-{\Pi }_{e_i}^{0}w.$$

$$\left|\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}w}{\partial n^2}\frac{\partial {\xi }_N}{\partial n}\mathrm{d}s\right|$$

$$\le C\left[{\left(\sum _{K\in {\mathcal{T}}_N}{h}_x^4{\Vert \frac{{\partial }^{4}w}{\partial x\partial y^3}\Vert }_{L^2(K)}^2\right)}^{1/2}+{\left(\sum _{K\in {\mathcal{T}}_N}{h}_y^4{\Vert \frac{{\partial }^{4}w}{\partial x^3\partial y}\Vert }_{L^2(K)}^2\right)}^{1/2}\right]|{\xi }_N{|}_{2,N}$$

$$\left|\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}w}{\partial n^2}\frac{\partial {\xi }_N}{\partial n}\mathrm{d}s\right|$$

$$\le C\left[\sum _{K\in {\mathcal{T}}_N}{\left(h_x^2{\Vert \frac{{\partial }^{4}w}{\partial x\partial y^3}\Vert }_{L^2(K)}^2\right)}^{1/2}+\sum _{K\in {\mathcal{T}}_N}{\left(h_y^2{\Vert \frac{{\partial }^{4}w}{\partial x^3\partial y}\Vert }_{L^2(K)}^2\right)}^{1/2}\right]|{\xi }_N{|}_{1,N}.$$

$$\sum _{K\in {\mathcal{T}}_N}{\int }_{\partial K}\frac{{\partial }^{2}w}{\partial n^2}\frac{\partial {\xi }_N}{\partial n}\mathrm{d}s$$

$$=\sum _{K\in {\mathcal{T}}_N}\sum _{i=1}^4{\int }_{e_i}\frac{{\partial }^{2}w}{\partial n^2}\left({\mathcal{R}}_{e_i}^0\frac{\partial {\xi }_N}{\partial n}\right)\mathrm{d}s=:\sum _{K\in {\mathcal{T}}_N}\sum _{i=1}^4{J}_i,\text{say.}$$

$$J_2+J_4$$

Full text

1 Introduction

Set $\Omega:=(0,1)^{2}\subset\mathbb{R}^{2}$, with closure $\bar{\Omega}$ and boundary $\partial\Omega$. The bending of a simply supported plate is modelled [10, 15, 26] by the following 4th-order singularly perturbed boundary value problem in which $u$ represents the deflection of the plate under the transverse loading $f$:

[TABLE]

where the parameter $\varepsilon$ (which is the ratio of bending rigidity to tensile stiffness of the plate) satisfies $0<\varepsilon\ll 1$, the function $c$ lies in $C(\bar{\Omega})$ with $c\geq c_{0}>0$ for some constant $c_{0}$ and $f$ lies in $C^{1}(\bar{\Omega})$. The vector $n=(n_{x},n_{y})^{T}$ is the outward-pointing unit normal vector to the boundary $\partial\Omega$ of $\Omega$, so $\partial u(x,y)/\partial n$ denotes the normal derivative on $\partial\Omega$.

During the past two decades, various finite element methods (FEMs) have ben proposed for the numerical solution of 4th-order singularly perturbed problems. Error analyses on quasi-uniform meshes using nonconforming FEMs are discussed in [5, 15, 19, 20, 24, 31, 33], and a $C^{0}$ interior penalty method in [2]. Because it is troublesome to impose global $C^{1}$ continuity on the computed solution in the 2D domain $\Omega$ for 4th-order problems, the $C^{0}$ interior penalty method relaxes the global $C^{1}$ continuity by adding penalty terms [2, eq.(2.3)] to a globally $C^{0}$ finite element space. Similarly, the discontinuous Galerkin method adds stabilisation terms [9, (3.4) and (3.13)] to the finite element space without imposing global $C^{0}$ continuity. Compared with these two approaches, nonconforming FEMs discard the global $C^{1}$ continuity condition but do not usually add any stabilisation terms.

The well-known nonconforming triangular Morley element for 4th-order problems was constructed in [23] and appears in many textbooks [3, 6, 13, 29]. This element does not however yield convergence for 2nd-order boundary value problems and is therefore unsuitable for our singularly perturbed problem (1.1); this issue is discussed in [24, Section 1]. Thus to solve (1.1), various modifications of the triangular Morley element have been proposed [24, 33].

An alternative departure from the triangular Morley element is the rectangular Morley element which was introduced in [32] and [29, Section 2.7.1]. It uses the same degrees of freedom as the triangular Morley element, viz., the point value at each node and the integral mean of the outer normal derivative on each edge. But unlike its triangular counterpart, the rectangular Morley element is convergent for both 4th-order and 2nd-order problems, so it is reasonable to consider it for solving (1.1), as has already been pointed out in [31].

Note also that the rectangular Morley element has only eight degrees of freedom in each element, which compares well with the twelve degrees of freedom of the nonconforming Adini element and the sixteen d.o.f. of the conforming Bogner-Fox-Schmit element (these two finite elements are also suitable for solving the problem (1.1)). We shall use the rectangular Morley element method to solve (1.1).

To solve (1.1) effectively, one must also use an appropriate mesh. Since the exact solution of (1.1) exhibits boundary layers (see Assumption 1.1 below), the mesh should be refined near $\partial\Omega$. Suitable meshes were used with an Adini nonconforming FEM in [21], a $C^{0}$ interior penalty method in [12], a mixed FEM in [11] and a mixed $hp$-FEM in [7]. One sees easily that the meshes most suitable for (1.1) are severely anisotropic, so one cannot use results from the literature that assume quasi-uniformity of the mesh, but mesh quasi-uniformity is assumed in many of the earlier papers that we mentioned earlier. Consequently, for example, while Lemmas 3.2 and 3.4 below were proved previously assuming mesh quasi-uniformity, we have to prove them on meshes without this assumption.

We shall use a tensor-product Shishkin mesh when solving (1.1). In an energy-type norm that is natural for this problem, we prove $O(\varepsilon^{1/2}N^{-1}+\varepsilon N^{-1}\ln N+N^{-3/2})$ convergence, where $N$ is the number of mesh elements in each coordinate direction and this convergence is uniform in the small parameter $\varepsilon$. This result is much better than the $O(N^{-1/2})$ rate that is attained on quasi-uniform meshes [31, Corollary 4.1].

Furthermore, comparing our Morley result with the nonconforming Adini element studied in [21], in the most troublesome regime when $\varepsilon\approx N^{-1}$, our method achieves optimal $O(N^{-3/2})$ convergence, which is better than the $O(N^{-1/2})$ attained in [21]; see Remark 4.16.

The paper is structured as follows. Section 1.1 gives a decomposition of the solution $u$ of (1.1) that is needed for the entire analysis. The rectangular Morley element is defined in Section 2. In Section 3 we define the discrete analogue of (1.1) that is solved by our FEM, and we derive various preliminary results for this discretisation. The analysis in this section is valid on arbitrary tensor-product meshes. The Shishkin mesh is defined in Section 4, then we continue our earlier analysis of the numerical method, finally proving its convergence in Theorem 4.14. The paper is completed by numerical experiments in Section 5.

Notation. Throughout the paper, $C$ denotes a generic constant that is independent of $\varepsilon$ and of the mesh. Standard Lebesgue spaces $L^{p}(\Omega^{\prime})$ and Sobolev spaces $H^{k}(\Omega^{\prime})$ are used, where $\Omega^{\prime}$ is any measurable subset of $\Omega$. Their associated norms are $\|\cdot\|_{L^{2}(\Omega^{\prime})}$, $\|\cdot\|_{H^{k}(\Omega^{\prime})}$ and the associated semi-norm is $\left|\cdot\right|_{H^{k}(\Omega^{\prime})}$. The inner product in $L^{2}(\Omega^{\prime})$ is denoted by $(\cdot,\cdot)_{\Omega^{\prime}}$, but we omit the subscript when $\Omega^{\prime}=\Omega$.

1.1 A decomposition of the solution $u$

If $c\in C(\bar{\Omega})\subset L^{2}(\Omega)$ and $f\in C^{1}(\bar{\Omega})\subset L^{2}(\Omega)$, then (cf. [13, Theorem 5.1]) one can use the Lax-Milgram lemma to show that (1.1) has a unique solution in $H^{2}(\Omega)$, but for our analysis we need more precise information about the behaviour of $u$, especially near the boundary $\partial\Omega$.

Number the sides of $\bar{\Omega}$ in clockwise order as $1,2,3,4$, where $1$ is associated with the interval $[0,1]$ on the $x$-axis. Then denote the corner where the sides $i$ and $j$ meet by $(i;j)$ with $i<j$. As in the papers [11, 21], we make the following assumption.

Assumption 1.1**.**

The solution $u$ of the boundary value problem (1.1) can be decomposed as

[TABLE]

for $(x,y)\in\bar{\Omega}$, where $S(x,y)$ is a smooth function, each $E_{i}(x,y)$ is a boundary layer component associated with the side $i$ of $\bar{\Omega}$, and each $E_{ij}(x,y)$ is a corner layer component associated with the corner $(i;j)$, i.e., there exists a constant $C$ such that for all $(x,y)\in\overline{\Omega}$ and $0\leq i+j\leq 4$ one has

[TABLE]

with analogous bounds for the remaining components of the decomposition.

2 The rectangular Morley element on anisotropic meshes

For the moment we use an arbitrary tensor product mesh on $\Omega$. In Section 4 we shall specialise our analysis to a suitable Shishkin mesh.

Let $N$ be a positive integer. Let $0=x_{0}<x_{1}\dots<x_{N}=1$ and $0=y_{0}<y_{1}\dots<y_{N}=1$ be arbitrary meshes in the interval $[0,1]$ on the $x$ and $y$ axes respectively. Draw vertical and horizontal lines through these points to partition $\bar{\Omega}$ into mesh rectangular elements with vertices $\{(x_{i},y_{j}):0\leq i,j\leq N\}$. We write $\mathcal{T}_{N}$ for the union of these mesh rectangles, so $\mathcal{T}_{N}$ is a partition of $\bar{\Omega}$.

Let $K$ be a typical mesh rectangle as shown in Figure 1. Its vertices are $a_{1},a_{2},a_{3},a_{4}$ and its edges are $e_{1},e_{2},e_{3},e_{4}$ with unit normal vectors $n_{e_{j}}$ for $j=1,2,3,4$.

Let $(x_{c},y_{c})$ be the center of $K$ and suppose $K$ has dimensions $2h_{x}\times 2h_{y}$. Then the coordinates of the vertices $a_{j}$ are $a_{1}=(x_{c}-h_{x},y_{c}-h_{y})$, $a_{2}=(x_{c}-h_{x},y_{c}+h_{y})$, $a_{3}=(x_{c}+h_{x},y_{c}+h_{y})$, $a_{4}=(x_{c}+h_{x},y_{c}-h_{y})$. For $j=1,2,3,4$, let $|e_{j}|$ denote the length of edge $e_{j}$ and write $n_{e_{j}}$ for the outward-pointing unit normal to $e_{j}$.

We now describe the rectangular Morley element of  [28, 29, 32]. For each mesh rectangle $K$, let $P(K)$ denote the span of $\{1,x,y,x^{2},xy,y^{2},x^{3},y^{3}\}$.

Define a set of 8 functions in $P(K)$ as follows:

[TABLE]

One can verify (see [22, eq.(2.5)]) that

[TABLE]

where $\delta_{ij}$ is the Kronecker symbol. Hence the $p_{i}$ and $q_{j}$ are linearly independent and form a basis for $P(K)$.

The rectangular Morley element space $V_{N}$ is now defined to be the set of functions $v\in L^{2}(\Omega)$ such that $v|_{K}\in P(K)$ for each mesh rectangle $K$, $v$ is continuous at all mesh vertices in $\Omega$ as one moves from one mesh rectangle to another, $v$ vanishes at all vertices in $\partial\Omega$, and for all $K$ one has $\int_{e}\left[\frac{\partial v}{\partial n_{e}}\right]\,\mathrm{d}s=0$ for all edges $e$ of $K$ where $n_{e}$ is the unit normal to $e$ that points out of $K$; here $\left[\frac{\partial v}{\partial n_{e}}\right]$ denotes the jump in $\frac{\partial v}{\partial n_{e}}$ as one moves from one side of $e$ to the other side (i.e., from one mesh rectangle to another) if $e\not\subset\partial\Omega$, while $\left[\frac{\partial v}{\partial n_{e}}\right]=\frac{\partial v}{\partial n}$ if $e\subset\partial\Omega$.

3 The discrete problem and its error analysis

In this section we discretise (1.1) using a finite element method based on the rectangular Morley element and begin the error analysis of this method.

3.1 Discretisation of (1.1)

Any function $v_{N}\in V_{N}$ may be only piecewise continuous on $\Omega$ so one cannot apply differential operators (such as $\nabla$ and $\Delta$) to $v_{N}$ globally. Thus, we apply these operators only piecewise on the interiors of the elements $K\in\mathcal{T}_{N}$ and indicate this modification by writing $\nabla_{N}$ and $\Delta_{N}$ instead of $\nabla$ and $\Delta$, e.g., $(\nabla_{N}v_{N})|_{K}=\nabla(v_{N}|_{K})$. For $v_{N},w_{N}\in V_{N}$, define the bilinear forms $b_{N}(v_{N},w_{N}):=(c\nabla_{N}v_{N},\nabla_{N}w_{N})$ and $a_{N}(v_{N},w_{N}):=(\Delta_{N}v_{N},\Delta_{N}w_{N})$. We shall also use expressions like $b_{N}(u,w_{N})$ and $a_{N}(u,w_{N})$, which are defined in the same way.

One can write (1.1) in the equivalent weak form

[TABLE]

where as usual H_{0}^{2}(\Omega)=\big{\{}g\in H^{2}(\Omega):g|_{\partial\Omega}=\frac{\partial g}{\partial n}|_{\partial\Omega}=0\big{\}}.

Then our discretisation of (1.1) is: Find $u_{N}\in V_{N}$ such that

[TABLE]

Existence and uniqueness of the solution $u_{N}$ of (3.1) follows readily from the Lax-Milgram lemma.

3.2 Preliminaries

For each $K\in{\mathcal{T}}_{N}$, let $P_{0}(K)$ be the space of constant functions defined on $K$ and let $P_{K}^{0}:L^{2}(K)\rightarrow P_{0}(K)$ be the $L^{2}$ projector, viz.,

[TABLE]

From [14, Lemma 7.5], for all $v\in H^{1}(K)$ one has the anisotropic error bound

[TABLE]

with $C$ independent of $v$ and $K$.

Let $Q_{1}(K)$ denote the usual space of polynomials of the form $c_{0}+c_{1}x+c_{2}y+c_{3}xy$ defined on $K$. Set Q_{N}=\big{\{}v_{N}\in H^{1}(\Omega):v_{N}|_{K}\in Q_{1}(K)\,\forall K\in\mathcal{T}_{N}\big{\}}. Let $I_{N}:C(K)\rightarrow Q_{N}$ denote the standard piecewise bilinear nodal interpolation operator, viz., given $v\in C(K)$, for each vertex $P$ of $K$ we set $I_{N}v(P)=v(P)$. Then for each $K$ and all $v\in H^{1}(K)\cap C(K)$, by [1, Theorem 2.7] one has

[TABLE]

and

[TABLE]

By applying a standard scaling argument to the $x$ and $y$ variables separately, one can modify the proof of [3, Lemma 4.5.3] to obtain the following anisotropic inverse inequalities: for each $K\in{\mathcal{T}}_{N}$ and any $\xi_{N}\in P(K)$, one has

[TABLE]

where the constants $C$ are independent of $K$ and $\xi_{N}$.

3.3 The error equation

Let $v_{N}\in V_{N}$ be arbitrary. Set $\xi_{N}=u_{N}-v_{N}$ and $\eta_{N}=u-v_{N}$. Then (3.1) and (1.1) yield

[TABLE]

Integrating by parts and using $I_{N}\xi_{N}=0$ on $\partial\Omega$, as in [18, (3.4)] we get

[TABLE]

Integration by parts on each $K\in\mathcal{T}_{N}$ gives, like [18, (3.5)],

[TABLE]

where $\partial/\partial s$ and $\partial/\partial n$ are, respectively, the tangential and outward-pointing normal derivatives along the boundary $\partial K$. From these calculcations we get the error equation

[TABLE]

Remark 3.1**.**

The error equation (3.6) is closely related to [31, (3.7)]; see also [18, (3.6)]. The analysis of (3.6) that we now develop has some similarity to the analysis in [31], but our work is greatly complicated (compared with [31]) by the extreme anisotropy of the Shishkin mesh and by our aim of proving error bounds that show convergence even when $\varepsilon$ is very small; in [31] the mesh is shape-regular and the bounds established in Theorems 4.1 and 4.2 blow up if one takes the limit as $\varepsilon\to 0$, while the bound of [31, Theorem 4.3] (which remains valid as $\varepsilon\to 0$) is only $O(N^{-1/2})$, which is inferior to the $O(N^{-3/2})$ bound that we shall prove in Theorem 4.14.

3.4 Bounds on error equation integrals along $\partial K$

In this subsection we bound the terms $\varepsilon^{2}\sum_{K\in\mathcal{T}_{N}}\int_{\partial K}\frac{\partial^{2}u}{\partial n^{2}}\frac{\partial\xi_{N}}{\partial n}\;\mathrm{d}s$ and $\varepsilon^{2}\sum_{K\in\mathcal{T}_{N}}\int_{\partial K}\frac{\partial^{2}u}{\partial s\partial n}\frac{\partial\xi_{N}}{\partial s}\;\mathrm{d}s$ from the error equation (3.6).

For all $v_{N}\in V_{N}$, define $|v_{N}|_{1,N}^{2}:=(\nabla_{N}v_{N},\nabla_{N}v_{N}),\,|v_{N}|_{2,N}^{2}:=(\Delta_{N}v_{N},\Delta_{N}v_{N})$ and the norm $\|v_{N}\|_{\varepsilon,N}^{2}:=\varepsilon^{2}|v_{N}|_{2,N}^{2}+|v_{N}|_{1,N}^{2}$. One can define $|u|_{1,N},|u|_{2,N}$ and $\|u\|_{\varepsilon,N}$ in the same way.

Let the edges of each $K\in\mathcal{T}_{N}$ be labelled $e_{i}$ for $i=1,\dots,4$ as in Figure 1. For each $w\in L^{2}(e_{i})$ define

[TABLE]

The argument used in the next lemma imitates in part the proof of [18, Lemma 3.2], where only a uniform mesh was considered.

Lemma 3.2**.**

Assume that $w\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$. Then there exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

For any internal edge $e_{i}$, the edge-jump condition in the definition of $V_{N}$ implies that $\Pi^{0}_{e_{i}}\frac{\partial\xi_{N}}{\partial n}$ has the same magnitude but different signs when calculated on the two adjacent mesh rectangles $K$ that share $e_{i}$. Consequently $\sum_{K\in\mathcal{T}_{N}}\sum_{i=1}^{4}\int_{e_{i}}\frac{\partial^{2}w}{\partial n^{2}}\left(\Pi^{0}_{e_{i}}\frac{\partial\xi_{N}}{\partial n}\right)\;\mathrm{d}s=0$, as in these integrals $\frac{\partial^{2}w}{\partial n^{2}}$ depends on $e_{i}$ on each internal edge but not on the associated $K$, while if $e_{i}\subset\partial\Omega$ then $\frac{\partial\xi_{N}}{\partial n}=0$. Hence

[TABLE]

We now imitate the proof of [27, Theorem 4.1] by considering $|J_{2}+J_{4}|$ and also borrow some techniques from [18, Lemma 3.2]. For any $K\in\mathcal{T}_{N}$, from Figure 1 one sees that $\frac{\partial\xi_{N}}{\partial n}=-\frac{\partial\xi_{N}}{\partial x}$ on $e_{2}$, $\frac{\partial\xi_{N}}{\partial n}=\frac{\partial\xi_{N}}{\partial x}$ on $e_{4}$ and $\frac{\partial^{2}w}{\partial n^{2}}=\frac{\partial^{2}w}{\partial x^{2}}$ on $e_{2}$ and $e_{4}$. Thus

[TABLE]

As $\xi_{N}\in V_{N}=\text{span}\{1,x,y,x^{2},xy,y^{2},x^{3},y^{3}\}$, it is easy to see that

[TABLE]

and the definition of $\mathcal{R}^{0}_{e_{2}}$ implies that $\int_{y_{c}-h_{y}}^{y_{c}+h_{y}}\left(\mathcal{R}^{0}_{e_{2}}\frac{\partial\xi_{N}}{\partial x}\right)(x_{c}-h_{x},y)\,\mathrm{d}y=0.$ Hence

[TABLE]

A Cauchy-Schwarz inequality now yields

[TABLE]

In the second integral here, one has

[TABLE]

so another Cauchy-Schwarz inequality gives

[TABLE]

where we used the property that $\frac{\partial^{2}\xi_{N}}{\partial x\partial s}$ is independent of $x$, which follows quickly from the definition of the Morley shape function space $P(K)$ in Section 2. Hence

[TABLE]

For the first integrand in (3.11), one has

[TABLE]

Again appealing to Cauchy-Schwarz, we get

[TABLE]

hence

[TABLE]

Substituting (3.12) and (3.13) into (3.11) then taking a square root yields

[TABLE]

One can show similarly that

[TABLE]

Combining (3.9), (3.14) and (3.15), a Cauchy-Schwarz inequality gives

[TABLE]

Thus, we have proved (3.7).

Applying the anisotropic inverse inequalities (3.5a) to (3.14) and (3.5b) to (3.15) gives

[TABLE]

Then, like the above derivation of (3.7), one obtains (3.8):

[TABLE]

∎

Remark 3.3**.**

At first sight (3.7), which bounds $\sum_{K\in\mathcal{T}_{N}}\int_{\partial K}\frac{\partial^{2}w}{\partial n^{2}}\frac{\partial\xi_{N}}{\partial n}\;\mathrm{d}s$ with factors $h_{x}^{2}$ and $h_{y}^{2}$, looks more attractive than (3.8), which yields factors $h_{x}$ and $h_{y}$. But to continue the analysis with (3.7) means working with $|\xi_{N}|_{2,N}$ and using the inequality $|\xi_{N}|_{2,N}\leq\varepsilon^{-1}\|\xi_{N}\|_{\varepsilon,N}$, which introduces an undesirable factor $\varepsilon^{-1}$ into our error bounds. Consequently we shall use (3.8) in our subsequent analysis, for then we have available $|\xi_{N}|_{1,N}\leq\|\xi_{N}\|_{\varepsilon,N}$, which causes no difficulties.

The next result was proved in [18, Lemma 3.3] while assuming that the mesh was quasi-uniform, but here we remove this restriction.

Lemma 3.4**.**

Assume that $w\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$. Then there exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

The definition of $V_{N}$ implies that on each edge $e_{i}$ that is shared by two mesh rectangles $K_{1}$ and $K_{2}$ (say), $\xi_{N}|_{e_{i}}$ is the same in both $K_{1}$ and $K_{2}$; but on $e_{i}$ the derivative $\frac{\partial^{2}w}{\partial s\partial n}$ will have opposite signs in these two rectangles, so

[TABLE]

Similarly to the derivation of (3.10), in the notation of Figure 1 we have

[TABLE]

As $\xi_{N}\in V_{N}=\text{span}\{1,x,y,x^{2},xy,y^{2},x^{3},y^{3}\}$, it is easy to see that

[TABLE]

and the definition of $\mathcal{R}^{0}_{e_{2}}$ implies that $\int_{y_{c}-h_{y}}^{y_{c}+h_{y}}\left(\mathcal{R}^{0}_{e_{2}}\frac{\partial\xi_{N}}{\partial y}\right)(x_{c}-h_{x},y)\;\mathrm{d}y=0$. Hence

[TABLE]

By a Cauchy-Schwarz inequality we get

[TABLE]

Replacing $\partial^{2}w/\partial x^{2}$ by $\partial^{2}w/\partial y\partial x$ in the derivation of (3.13) yields

[TABLE]

while similarly to the derivation of (3.12), one has

[TABLE]

Hence (3.19) implies that

[TABLE]

One can show similarly that

[TABLE]

From (3.18), (3.20) and (3.21) and some Cauchy-Schwarz inequalities, we can derive

[TABLE]

which proves (3.16).

Invoking the anisotropic inverse inequalities (3.5a)–(3.5b), from (3.20) and (3.21) one has

[TABLE]

Then one obtains (3.17) by imitating the above derivation of (3.16). ∎

Lemmas 3.2 and 3.4 are anisotropic generalisations on nonuniform meshes of the consistency error bounds proved in [17, 18] for discretisations of a classical biharmonic problem ($\varepsilon=1$ and $c\equiv 0$ in (1.1)) on uniform meshes.

4 Error analysis on a Shishkin mesh

All our analysis so far is valid on the general tensor-product mesh defined at the start of Section 2, but to obtain good accuracy in our numerical solution of (1.1) one should use a mesh that is designed to handle boundary layers in singularly perturbed problems. We now restrict our attention to a well-known tensor-product mesh of this type: the Shishkin mesh.

4.1 The Shishkin mesh

Recall that $N$ is the number of mesh intervals in each coordinate direction. Assume from now on that $N$ is divisible by $4$. Define the transition parameter

[TABLE]

We assume that $\lambda=\varepsilon\ln N$ as otherwise $N$ would be exponentially large compared with $1/\varepsilon$, and then the problem can be solved accurately using a uniform mesh. In practice $\varepsilon$ is small, so it is extremely unlikely that $\varepsilon\ln N>1/4$.

Partition the interval $[0,1]$ by a piecewise uniform mesh that is obtained by dividing $[0,\lambda]$ into $N/4$ equal subintervals, $[\lambda,1-\lambda]$ into $N/2$ equal subintervals and $[1-\lambda,1]$ into $N/4$ equal subintervals. The tensor product of two such 1D meshes is our 2D Shishkin mesh for the problem (1.1); it is shown in Figure 2. (See [25, 30] for further discussion of Shishkin meshes.) We continue to write $\mathcal{T}_{N}$ for the 2D mesh. For any mesh sub-domain $\tilde{\Omega}$ of $\Omega$, denote the mesh on $\tilde{\Omega}$ by $\mathcal{T}_{N}(\tilde{\Omega})$.

In our error analysis, we use the Shishkin mesh subdomains $\Omega_{1}$,…, $\Omega_{9}$ that are displayed in Figure 3. For convenience, define $\widehat{\Omega}_{1}:=\Omega_{1}\cup\Omega_{4}\cup\Omega_{7}$, $\widehat{\Omega}_{2}:=\Omega_{2}\cup\Omega_{5}\cup\Omega_{8}$ and $\widehat{\Omega}_{3}:=\Omega_{3}\cup\Omega_{6}\cup\Omega_{9}$.

We shall now derive an error bound for the numerical solution $u_{N}$ of (3.1) by specialising our earlier results to the Shishkin mesh. To do this we analyse the consistency error of the discretisation of each term in the PDE (1.1a) in several lemmas, then combine these bounds in the proof of Theorem 4.14 to give the final convergence result.

4.2 Consistency error of $a_{N}(\cdot,\cdot)$

Lemma 4.1**.**

Suppose that $u\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$ satisfies Assumption 1.1. Then there exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

We apply the bound (3.8) of Lemma 3.2 to each of the components of $u$ that are provided by Assumption 1.1.

For the smooth component $S$, the bound (1.2a) and $h_{x}\leq 2N^{-1}$ yield

[TABLE]

For the layer component $E_{1}$, using (1.2b) on $\widehat{\Omega}_{2}$ we get

[TABLE]

and on $\widehat{\Omega}_{1}\cup\widehat{\Omega}_{3}$

[TABLE]

Combining this pair of estimates and recalling that $\varepsilon\ln N\leq 1/4$, we obtain

[TABLE]

The estimates for $E_{2}$, $E_{3}$ and $E_{4}$ are similar.

For the layer component $E_{12}$, using (1.2c) on $\widehat{\Omega}_{2}$ gives

[TABLE]

and on $\widehat{\Omega}_{1}\cup\widehat{\Omega}_{3}$

[TABLE]

Thus,

[TABLE]

The estimates for $E_{23}$, $E_{34}$ and $E_{41}$ are similar.

Adding (4.1), (4.2) and (4.3), we get

[TABLE]

Similarly,

[TABLE]

Now the bound (3.7) of Lemma 3.2 and $|\xi_{N}|_{1,N}\leq\|\xi_{N}\|_{\varepsilon,N}$ yield

[TABLE]

∎

Lemma 4.2**.**

Suppose that $u\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$ satisfies Assumption 1.1. Then there exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

Imitate the proof of Lemma 4.1 but invoke Lemma 3.4 instead of Lemma 3.2. To be specific, Lemma 3.4 gives bounds on the smooth and corner layer components that are the same as (4.1) and (4.3), while for the boundary layer components, the bound $C\varepsilon^{1/2}N^{-1}$ of (4.2) is replaced by the smaller bound $C\varepsilon^{3/2}N^{-1}$. ∎

Remark 4.3**.**

[Comparison with Adini element] For the highest-order term $\varepsilon^{2}\Delta^{2}u(x,y)$ of (1.1a), the consistency error of the rectangular Morley element on the Shishkin mesh is better than the corresponding bound for the rectangular Adini element, which was discussed in [21]. For in [21, eq.(4.27)], the consistency error of the Adini element is bounded by

[TABLE]

or

[TABLE]

after using an inverse inequality. In the above bounds, terms such as $h_{x}^{2}\left\|\frac{\partial^{4}u}{\partial y^{4}}\right\|_{L^{2}(K)}$ or $h_{x}\left\|\frac{\partial^{4}u}{\partial y^{4}}\right\|_{L^{2}(K)}$ are troublesome because of the mismatch where derivatives in one variable are multiplied by mesh widths in the other variable—this is problematic for example on mesh rectangles $K\subset\Omega_{2}$ where $|\partial^{4}u/\partial y^{4}|=O(\varepsilon^{-3})$ by (1.2b) but $h_{x}=O(N^{-1})$, leading to bounds involving negative powers of $\varepsilon$. Remarkably, the rectangular Morley element does not have this drawback; the estimates of Lemmas 3.2 and 3.4 do not have similar mismatches and consequently we are able to prove the error bounds of Lemmas 4.1 and 4.2 in which no negative powers of $\varepsilon$ appear.

Recall the piecewise bilinear interpolation operator $I_{N}$ of Section 3.1 and the piecewise constant projector $P^{0}_{K}$ of Section 3.2.

Lemma 4.4**.**

Suppose that $u\in H_{0}^{2}(\Omega)\cap H^{4}(\Omega)$ satisfies Assumption 1.1. Then there exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

Let $\xi_{N}\in V_{N}$ and $K\in\mathcal{T}_{N}$. Similarly to the proof of [17, Lemma 3.1], one sees that

[TABLE]

by verifying this property for each of the functions $1,x,y,x^{2},xy,y^{2},x^{3},y^{3}$. Hence

[TABLE]

where we used (3.4) and Cauchy-Schwarz inequalities. But the anisotropic projection error bound (3.2) yields

[TABLE]

Here the worst-behaved terms are $h_{x}^{2}\left\|\frac{\partial^{4}u}{\partial x\partial y^{3}}\right\|_{L^{2}(K)}^{2}$ and $h_{y}^{2}\left\|\frac{\partial^{4}u}{\partial x^{3}\partial y}\right\|_{L^{2}(K)}^{2}$. These terms can be bounded like (4.4), obtaining

[TABLE]

The lemma now follows since $\left(\sum_{K\in\mathcal{T}_{N}}\|\nabla\xi_{N}\|_{L^{2}(K)}^{2}\right)^{1/2}\leq\|\xi_{N}\|_{\varepsilon,N}$. ∎

4.3 Consistency error of $b_{N}(\cdot,\cdot)$

In this subsection we consider the consistency error of the lower-order term in the PDE (1.1a). The components $u^{E}$ and $S$ of $u$ that were defined in Assumption 1.1 are handled separately.

Lemma 4.5**.**

There exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

Recall that $u^{E}(x,y)=\sum_{i=1}^{4}E_{i}(x,y)+E_{12}(x,y)+E_{23}(x,y)+E_{34}(x,y)+E_{14}(x,y)$. Similarly to the proof of Lemma 4.4, one has

[TABLE]

as $\left|\nabla u^{E}\right|\leq C$ on $\Omega$ by Assumption 1.1 and $\left(\sum_{K\in\mathcal{T}_{N}}\|\nabla\xi_{N}\|_{L^{2}(K)}^{2}\right)^{1/2}\leq\|\xi_{N}\|_{\varepsilon,N}$. These two terms have a similar structure, but $c$ is well behaved so we give a detailed analysis only for the term containing $u^{E}$.

The anisotropic projection error estimate (3.2) yields

[TABLE]

The boundary layer and corner layer components of $u^{E}$ are handled as in the proof of Lemma 4.1. The worst-behaved terms are $h_{x}^{2}\left\|\frac{\partial^{2}u^{E}}{\partial x\partial y}\right\|_{L^{2}(K)}^{2}$ and $h_{y}^{2}\left\|\frac{\partial^{2}u^{E}}{\partial x\partial y}\right\|_{L^{2}(K)}^{2}$, whereas in Lemma 4.1 one has $h_{x}^{2}\left\|\frac{\partial^{4}u^{E}}{\partial x\partial y^{3}}\right\|_{L^{2}(K)}^{2}$; that is, here there are two fewer orders of derivative but this is balanced by the additional factor $\varepsilon^{2}$ that is present in Lemma 4.1. Thus, like (4.2), we get finally

[TABLE]

and like (4.3) we get

[TABLE]

These estimates yield

[TABLE]

Combining these bounds, we are done. ∎

It now remains only to bound the smooth component term $\left|b_{N}(S,\xi_{N}-I_{N}\xi_{N})\right|$, but obtaining a sharp estimate for this term turns out to be troublesome.

For the rectangular Shishkin mesh $\mathcal{T}_{N}$, let $e$ be an edge of some mesh element $K$. If $e\not\subset\partial\Omega$, then there are two elements $K_{L}$ and $K_{R}$ in $\mathcal{T}_{N}$ that share $e$ as a common edge. Let $\omega_{e}:=K_{L}\cup K_{R}$ denote the patch associated with $e$. Let $|K_{L}|$ and $|K_{R}|$ denote the areas of $K_{L}$ and $K_{R}$ respectively. If $|K_{L}|=|K_{R}|$, then we call the patch uniform.

Let $\mathcal{E}_{N}=\mathcal{E}_{N}^{1}\bigcup\mathcal{E}_{N}^{2}\bigcup\mathcal{E}_{N}^{3}$ be the set of all edges of elements in $\mathcal{T}_{N}$, where $\mathcal{E}_{N}^{1}$ contains the internal edges with uniform patches, $\mathcal{E}_{N}^{2}$ contains the internal edges with nonuniform patches, and $\mathcal{E}_{N}^{3}$ contains the edges in $\partial\Omega$.

Recall the properties (2.1) of the Morley basis. Let $\mathcal{N}_{N}:=\{(x_{i},y_{j}):i,j=0,\dots,N\}$ denote the set of nodes of $\mathcal{T}_{N}$. For each $v_{N}\in V_{N}$, we shall treat separately the components of $v_{N}$ that correspond to nodes or edges. Thus, set $V_{N}^{\mathcal{X}}:=\left\{v_{N}\in V_{N}:\int_{e}\frac{\partial v_{N}}{\partial\mathbf{n}_{e}}ds=0\ \forall\,e\in\mathcal{E}_{N}\right\}$ and $V_{N}^{\mathcal{E}}:=\left\{v_{N}\in V_{N}:v_{N}(x)=0\ \forall\,x\in\mathcal{N}_{N}\right\}$.

We decompose $V_{N}^{\mathcal{E}}$ further. For each $e\in\mathcal{E}_{N}$, define

[TABLE]

Set $P^{\mathcal{X}}(K)={\rm span}\{p_{i},1\leqslant i\leqslant 4\}$, where the $p_{i}$ were defined in Section 2. For each mesh element $K$ and each edge $e$ of $K$, set

[TABLE]

If $e$ is an internal edge with patch $\omega_{e}=K_{L}\cup K_{R}$, set

[TABLE]

The next two lemmas establish special properties of the rectangular Morley element.

Lemma 4.6**.**

Let $K\in{\mathcal{T}}_{N}$ be a mesh element. Let $\phi\in P^{\mathcal{X}}(K)$ and $q_{K}^{1}\in Q_{1}(K)$. Then

[TABLE]

Proof.

Recall that $(x_{c},y_{c})$ is the center of $K$ and $K$ has dimensions $2h_{x}\times 2h_{y}$. To simplify the presentation, set $\xi_{x}=(x-x_{c})/h_{x}$ and $\xi_{y}=(y-y_{c})/h_{y}$. One can deduce from the definitions of the $p_{i}$ in Section 2 that for any $\phi\in P^{\mathcal{X}}(K)$, one has

[TABLE]

Hence

[TABLE]

and similarly for $\partial(\phi-I_{N}\phi)/\partial y$. But

[TABLE]

so $\int_{K}q_{K}^{1}\frac{\partial(\phi-I_{N}\phi)}{\partial x}\,\mathrm{d}x\,\mathrm{d}y=0$ follows. The other identity in (4.6) is proved similarly. ∎

Let $Q_{1}(\omega_{e})$ denote the space of polynomials of the form $c_{0}+c_{1}x+c_{2}y+c_{3}xy$ defined on $\omega_{e}$ and let $P_{0}(\omega_{e})$ denote the space of constant functions defined on $\omega_{e}$.

Lemma 4.7**.**

Let $e\in\mathcal{E}_{N}$ be an edge.

1. (i)

If $e\in\mathcal{E}_{N}^{1}\bigcup\mathcal{E}_{N}^{2}$, its associated patch $\omega_{e}$ is uniform, $\psi\in P^{e}(\omega_{e})$ and $q_{\omega_{e}}^{1}\in Q_{1}(\omega_{e})$, then

[TABLE] 2. (ii)

If $e\in\mathcal{E}_{N}^{1}\bigcup\mathcal{E}_{N}^{2}$, its associated patch $\omega_{e}$ is nonuniform, $\psi\in P^{e}(\omega_{e})$ and $p_{\omega_{e}}^{0}\in P_{0}(\omega_{e})$, then

[TABLE] 3. (iii)

If $e\in\mathcal{E}_{N}^{3}$ with $e\subset\partial K$, $\psi\in P^{e}(K)$ and $p_{K}^{0}\in P_{0}(K)$, then

[TABLE]

Proof.

Assume that the edge $e$ is parallel to the $y$-coordinate axis. (The case where $e$ is parallel to the $x$-coordinate axis is similar.) Let the center of $e$ be the point $(e_{x,c},e_{y,c})$ and let the half-lengths of $K_{L}$ and $K_{R}$ in the $x$ direction be $h_{x}^{L}$ and $h_{x}^{R}$ respectively. Then $e=\{(e_{x,c},y):y=e_{y,c}+\xi_{y}h_{y},\;-1\leq\xi_{y}\leq 1\}$ and

[TABLE]

On $\omega_{e}$, take increasing $x$ as the positive direction. Then (see Section 2) the basis function $\psi^{\omega_{e}}$ associated with $\omega_{e}$ is like $q_{4}$ in $K_{L}$ and $-q_{2}$ in $K_{R}$:

[TABLE]

Clearly $\psi^{\omega_{e}}=0$ on the right and left edges of $K_{L}$ and of $K_{R}$. Hence $I_{N}\psi^{\omega_{e}}\equiv 0$ and one gets

[TABLE]

and

[TABLE]

from which (4.8) and (4.9) follow.

One also has $\partial\psi^{\omega_{e}}/\partial y\equiv 0$, so $\int_{\omega_{e}}q_{\omega_{e}}^{1}\frac{\partial\left(\psi^{\omega_{e}}-I_{N}\psi^{\omega_{e}}\right)}{\partial y}\,\mathrm{d}x\,\mathrm{d}y=0$. Moreover, if $h_{x}^{L}=h_{x}^{R}$ then from

[TABLE]

we get $\int_{\omega_{e}}q_{\omega_{e}}^{1}\frac{\partial\left(\psi^{\omega_{e}}-I_{N}\psi^{\omega_{e}}\right)}{\partial x}\,\mathrm{d}x\,\mathrm{d}y=0$, which yields (4.7). ∎

Let $I_{\omega_{e}}:C(\omega_{e})\rightarrow Q_{1}(\omega_{e})$ denote the bilinear nodal interpolation operator on the patch $\omega_{e}$ that, given $v\in C(\omega_{e})$, satisfies $I_{\omega_{e}}v(P)=v(P)$ for each of the four vertices $P$ of $\omega_{e}$ that do not lie on $e$. Let $P_{\omega_{e}}^{0}:L^{2}(\omega_{e})\rightarrow P_{0}(\omega_{e})$ be the $L^{2}$ projector on patch $\omega_{e}$, viz.,

[TABLE]

Lemma 4.8**.**

There exists a constant $C$ such that for all $\xi_{N}\in V_{N}$, one has

[TABLE]

Proof.

For each $\xi_{N}\in V_{N}$, by [22, Lemma 3.6] there exist $\xi_{N}^{\mathcal{X}}\in V_{N}^{\mathcal{X}}$ and $\xi_{N}^{e}\in V_{N}^{e}\,\forall e\in\mathcal{E}_{N}$ such that $\xi_{N}=\xi_{N}^{\mathcal{X}}+\sum_{e\in\mathcal{E}_{N}}\xi_{N}^{e}$ and

[TABLE]

Consequently

[TABLE]

We will deal with the right-hand side terms one by one. First, using (4.6), Cauchy-Schwarz inequalities and (1.2a) one has

[TABLE]

For the terms in $\sum_{e\in\mathcal{E}_{N}}\dots$ in (4.11), consider the three cases $e\in\mathcal{E}_{N}^{1}$, $e\in\mathcal{E}_{N}^{2}$ and $e\in\mathcal{E}_{N}^{3}$. The definition of $V_{N}$ implies that $\xi_{N}^{e}=0$ if $e\in\mathcal{E}_{N}^{3}$, so we need only consider $e\in\mathcal{E}_{N}^{1}$ and $e\in\mathcal{E}_{N}^{2}$.

Using (4.7), Cauchy-Schwarz inequalities and (1.2a) we get

[TABLE]

Using (4.8), (1.2a), (3.2) and (3.3), we obtain

[TABLE]

as $\text{meas}(\omega_{e})=O(N^{-2})$ and there are $O(N)$ edges in $\mathcal{E}_{N}^{2}$.

Combining (4.11)–(4.14) then recalling (4.10) gives $\left|b_{N}(S,\xi_{N}-I_{N}\xi_{N})\right|\leq CN^{-3/2}\|\xi_{N}\|_{\varepsilon,N}$. ∎

Lemma 4.9**.**

There exists a constant $C$ such that for all $\xi_{N}\in V_{N}$ one has

[TABLE]

Proof.

Using the decomposition of $\xi_{N}\in V_{N}$ in (4.10), we have

[TABLE]

We will deal with these terms one by one.

Integrating by parts in (4.6), one sees that on each element $K\in{\mathcal{T}}_{N}$, for every $\phi\in P^{\mathcal{X}}(K)$ and every $p_{K}^{0}\in P_{0}(K)$ one has $\int_{K}p_{K}^{0}(\phi-I_{N}\phi)\,\mathrm{d}x\,\mathrm{d}y=0$. Hence we have

[TABLE]

where we used (3.2) and (3.3).

Similarly, integrating by parts in (4.7), we find that if $e\in\mathcal{E}_{N}^{1}\bigcup\mathcal{E}_{N}^{2}$, its associated patch $\omega_{e}$ is uniform, $\psi\in P^{e}(\omega_{e})$ and $p_{\omega_{e}}^{0}\in P_{0}(\omega_{e})$, then $\int_{\omega_{e}}p_{\omega_{e}}^{0}(\psi-I_{N}\psi)\,\mathrm{d}x\,\mathrm{d}y=0$. Hence

[TABLE]

using a variant of (3.2) to bound $\left\|f-P_{\omega_{e}}^{0}f\right\|_{L^{2}(\omega_{e})}$ and the definition of $V_{N}^{\mathcal{E}}$ to deduce that $I_{N}\xi_{N}^{e}=0$.

Using Cauchy-Schwarz inequalities and (3.3) we have

[TABLE]

as $\text{meas}(\omega_{e})=O(N^{-2})$ and there are $O(N)$ edges in $\mathcal{E}_{N}^{2}$.

Combining (4.16), (4.17) and (4.18) with (4.15) and using (4.10) we could get

[TABLE]

∎

Lemma 4.10**.**

Let $v_{N},\xi_{N}\in V_{N}$ be arbitrary and set $\eta_{N}=u-v_{N}$. Then there exists a constant $C$, which is independent of $v_{N}$ and $\xi_{N}$, such that

[TABLE]

Proof.

Two Cauchy-Schwarz inequalities and the definition of $|\cdot|_{2,N}$ give

[TABLE]

The result now follows from the definition of $\|\xi_{N}\|_{\varepsilon,N}$. ∎

Lemma 4.11**.**

Let $v_{N},\xi_{N}\in V_{N}$ be arbitrary and set $\eta_{N}=u-v_{N}$. Then there exists a constant $C$, which is independent of $v_{N}$ and $\xi_{N}$, such that

[TABLE]

Proof.

Two Cauchy-Schwarz inequalities and the definition of $|\xi_{N}|_{\varepsilon,N}$ give

[TABLE]

∎

4.4 Reduced Morley interpolation

To finish the convergence analysis, we define an interpolation operator $\Pi_{K}^{-}$ that uses reduced forms of the basis functions $p_{i},q_{i}\,(i=1,2,3,4)$ of Section 2. One does not use the full basis functions because anisotropic equivalents of Lemmas 4.12 and 4.13 are not known for them.

As in [28], on each mesh rectangle $K$ as in Figure 1 and each $v\in C^{1}(K)$, define the reduced Morley interpolant $\Pi_{K}^{-}v(x,y)$ for $(x,y)\in K$ by

[TABLE]

where we set

[TABLE]

for $(x,y)\in K$. This set of functions is obtained by removing all the cubic polynomials from the basis functions of the standard rectangular Morley basis of Section 2. For each mesh rectangle $K$, let $P^{-}(K)$ denote the space spanned by $\{p_{i}^{-},q_{i}^{-}:i=1,2,3,4\}$.

Lemma 4.12**.**

For each $K\in{\mathcal{T}}_{N}$ and any $v\in H^{2}(K)$, there exist constants $C$ (which are independent of $K$ and $v$) such that

[TABLE]

and

[TABLE]

For each $K\in{\mathcal{T}}_{N}$ and any $v\in H^{3}(K)$, there exist constants $C$ (which are independent of $K$ and $v$) such that

[TABLE]

and

[TABLE]

Proof.

We combine [4, Theorem 2.3] with a standard scaling argument to prove (4.19) and (4.20), imitating the analysis of the Wilson element in [4, Section 3]. Let $\hat{K}:=[-1,1]^{2}$ be the standard reference element. Define the nodes $\hat{a}_{i}$ and edges $\hat{e}_{j}$ of $\hat{K}$ analogously to Figure 1. Given any $\hat{v}\in C^{1}(\hat{K})$, its reduced rectangular Morley interpolant is defined by

[TABLE]

where we set

[TABLE]

We check that $\Pi_{\hat{K}}^{-}$ satisfies the hypotheses of [4, Theorem 2.3] for $|\alpha|=1$, $p=q=2$, $l=m=0$ and $r=3$. Let $D^{\alpha}$ denote $\partial/\partial x$ or $\partial/\partial y$. It is easy to see that $P_{0}(\hat{K})\subset D^{\alpha}P^{-}(\hat{K})=P_{1}(\hat{K})$ and $\Pi_{\hat{K}}^{-}$ is a continuous linear mapping from $H^{2}(\hat{K})$ to $H^{1}(\hat{K})$. We must verify [4, eqs.(2.14) and (2.15)]. We shall do this for the case $\alpha=(1,0)$ (that is, $D^{\alpha}=\partial/\partial x$); the other case $\alpha=(0,1)$ is similar. A computation yields

[TABLE]

where $\beta_{i}(\hat{v}):=F_{i}(\frac{\partial\hat{v}}{\partial\hat{x}})$ for $i=1,2,3$ with $F_{1}(\hat{w}):=\frac{1}{4}\int_{\hat{K}}\frac{\partial\hat{w}}{\partial\hat{y}}\,d\hat{x}\,d\hat{y},\ F_{2}(\hat{w}):=\frac{1}{4}\int_{\hat{K}}\frac{\partial\hat{w}}{\partial\hat{x}}\,d\hat{x}\,d\hat{y}$,

[TABLE]

and $\hat{q}_{1}(\hat{x},\hat{y}):=\hat{y},\,\hat{q}_{2}(\hat{x},\hat{y}):=\hat{x},\,\hat{q}_{3}(\hat{x},\hat{y}):=1$. Clearly $\{\hat{q}_{1},\hat{q}_{2},\hat{q}_{3}\}$ is a basis for $\hat{D}^{\alpha}P^{-}(\hat{K})=P_{1}(\hat{K})$, so we have now verified [4, eq.(2.14)]. To verify [4, eq.(2.15)] is straightforward: Cauchy-Schwarz inequalities yield $\left|F_{i}(\hat{w})\right|\leq\hat{C}\|\hat{w}\|_{H^{1}(\hat{K})}$ for $i=1,2,3$.

One can now invoke [4, Theorem 2.3] for $|\alpha|=1$, $p=q=2$, $l=m=0$ and $r=3$. Then a standard scaling argument gives (4.19) and (4.20).

The interpolation error estimates (4.21) and (4.22) can be proved similarly, except that now $l=1$ (we had $l=0$ for (4.19) and (4.20)), while the other parameters are unchanged: $|\alpha|=1$, $p=q=2$, $m=0$ and $r=3$. It is easy to see that $P_{1}(\hat{K})\subset D^{\alpha}P^{-}(\hat{K})=P_{1}(\hat{K})$ and $\Pi_{\hat{K}}^{-}$ is a continuous linear mapping from $H^{3}(\hat{K})$ to $H^{1}(\hat{K})$. The verification of [4, eq.(2.14)] is same as before and the verification of [4, eq.(2.15)] is similar since $\left|F_{i}(\hat{w})\right|\leq\hat{C}\|\hat{w}\|_{H^{1}(\hat{K})}\leq\hat{C}\|\hat{w}\|_{H^{2}(\hat{K})}$ for $i=1,2,3$. Then we invoke [4, Theorem 2.3], followed by a standard scaling argument, thereby obtaining (4.21) and (4.22). ∎

The next result is the analogue of Lemma 4.12 for 2nd-order derivatives.

Lemma 4.13**.**

For each $K\in{\mathcal{T}}_{N}$ and any $v\in H^{3}(K)$, there exist constants $C$ (independent of $K$ and $v$) such that

[TABLE]

[TABLE]

and

[TABLE]

Proof.

Apply the first inequality of [28, eq.(12)] on the reference element $\hat{K}=[-1,1]^{2}$, then use a standard scaling argument. ∎

4.5 Error bound on the Shishkin mesh

We can now complete the error analysis of our numerical method.

Theorem 4.14**.**

Let $u$ be the solution of (1.1) and $u_{N}$ the solution of (3.1). Let Assumption 1.1 be satisfied. Then there exists a constant $C$ (which is independent of $\varepsilon$ and $N$) such that

[TABLE]

Proof.

We imitate the argument of [28, Theorem 1]. Define $v_{N}\in V_{N}$ by $v_{N}=\Pi_{K}^{-}u$ on each $K\in\mathcal{T}_{N}$. Set $\xi_{N}=u_{N}-v_{N}$. and $\eta_{N}=u-v_{N}$. Recalling the error equation (3.6), we estimate the terms on its right-hand side using Lemmas 4.1, 4.2 and 4.4, 4.5, 4.8 and the definitions of $a_{N}$ and $b_{N}$: this gives

[TABLE]

for some constants $C$. But $u-u_{N}=u-v_{N}+v_{N}-u_{N}=\eta_{N}-\xi_{N}$, so

[TABLE]

from the definitions of $\left\|\cdot\right\|_{\varepsilon,N}$ and $\eta_{N}$.

The bounds of Lemma 4.13 yield

[TABLE]

To estimate each of these terms, one can proceed similarly to the proof of Lemma 4.1 above. Specifically, the worst-behaved terms are $h_{x}^{2}\left\|\frac{\partial^{3}u}{\partial x\partial y^{2}}\right\|_{L^{2}(K)}^{2}$ and $h_{y}^{2}\left\|\frac{\partial^{3}u}{\partial x^{2}\partial y}\right\|_{L^{2}(K)}^{2}$. Comparing $h_{x}^{2}\left\|\frac{\partial^{3}u}{\partial x\partial y^{2}}\right\|_{L^{2}(K)}^{2}$ with the term $h_{x}^{2}\left\|\frac{\partial^{4}u}{\partial x\partial y^{3}}\right\|_{L^{2}(K)}^{2}$ in Lemma 4.1, here the order of the derivative is one less, which effectively removes a factor $\varepsilon^{-1}$ from the calculation leading to (4.4), so instead of (4.4) one obtains

[TABLE]

Next, for the boundary and corner layer components, the bounds of Lemma 4.12 give

[TABLE]

by (4.3). For the smooth component, by (4.21) and (4.22) we have

[TABLE]

Substituting (4.27), (4.28) and (4.29) into (4.26), we are done. ∎

It is well known that in the analysis of singularly perturbed problems on layer-adapted meshes such as the Shishkin mesh, the most troublesome regime is when $\varepsilon\approx N^{-1}$. For this regime, Theorem 4.14 gives the following result.

Corollary 4.15**.**

Let $u$ be the solution of (1.1) and $u_{N}$ the solution of (3.1). Let Assumption 1.1 be satisfied. If $\varepsilon\approx N^{-1}$, then there exists a constant $C$ (which is independent of $\varepsilon$ and $N$) such that

[TABLE]

Remark 4.16** (Adini versus Morley).**

The Adini element is a well-known example of a nonconforming finite element that is suitable for our singularly perturbed problem (1.1). For the Adini element on an appropriate Shishkin mesh, one has the error bound [21, Corollary 4.1]

[TABLE]

Thus in the most challenging regime when $\varepsilon\approx N^{-1}$, the Adini element obtains $O(N^{-1/2})$ convergence, but Corollary 4.15 shows that the rectangular Morley element attains $O(N^{-3/2})$ convergence.

Remark 4.17**.**

An inspection of the entire analysis leading to the result of Theorem 4.14 shows that we are close to proving the sharper error bound

[TABLE]

which would agree with our numerical results in Section 5. The weaker term $N^{-3/2}$ in Theorem 4.14 appears only because of the two inequalities (4.14) and (4.18)—all other steps of the analysis yield terms in (4.30). Unfortunately we are unable at present to improve (4.14) and (4.18).

5 Numerical experiments

To test the accuracy of our finite element method we present three numerical examples, one with a known solution that displays typical layer behaviour and two others with unknown solutions.

Example 5.1**.**

In (1.1) choose $c(x,y)\equiv 1$ and choose $f(x,y)$ such that the exact solution is $u(x,y)=g(x)h(y)$, where

[TABLE]

and

[TABLE]

with $l=1-e^{-1/\varepsilon}$, $q=2-l$ and $d=1/(q-2\varepsilon l)$.

The derivatives of $u$ in this example match the bounds of Assumption 1.1. The solution computed by our finite element method (3.1) is $u_{N}$. Table 1 displays the errors $\|u-u_{N}\|_{\varepsilon,N}$ for various values of $\varepsilon$ and $N$, and the corresponding convergence rates

[TABLE]

where the convergence is assumed to be $O(N^{-r})$ for each fixed value of $\varepsilon$, i.e., the convergence rate is measured along each row of the table.

To test the case $\varepsilon\approx N^{-1}$ that was discussed in Remark 4.16, we choose $\varepsilon=N^{-1}$ in Example 5.1. Table 2 displays the errors $\|u-u_{N}\|_{N^{-1},N}$ for various values of $N$, and the corresponding convergence rates

[TABLE]

where the convergence is assumed to be $O(N^{-r})$. The table shows that our method attains $O(N^{-3/2})$ convergence, in agreement with Corollary 4.15.

Example 5.2**.**

In (1.1) we impose the transverse loading $f(x,y)=100(1-x-y+2xy)(x+y-2xy)$, which is the same as the numerical example of [11, equation (15)], and choose $c(x,y)=3+(1+x)y^{2}+(2-y)e^{x}$. The exact solution of this problem is unknown.

Example 5.3**.**

In (1.1) we impose the transverse loading $f(x,y)=2\pi^{2}[1-\cos(2\pi x)\cos(2\pi y)]$ and choose $c(x,y)=1$; this numerical example was tested in [31, Example 5.3] and [16, Example 5.3]. Its exact solution is unknown.

Numerical results for Example 5.2 and 5.3 are shown in Table 3 and 4. As the exact solution is unknown, we use the double mesh principle [8, Chapter 8] to estimate the errors and rates of convergence. To be specific, in the table $\|u-u_{N}\|_{\varepsilon,N}$ is replaced by $\|\tilde{u}_{2N}-u_{N}\|_{\varepsilon,N}$, where $\tilde{u}_{2N}$ is the solution computed by the same method on a Shishkin-type mesh that is constructed by bisecting each element $K\in\mathcal{T}_{N}$ in the $x$ and $y$ directions, i.e., dividing each $K$ into 4 rectangles of equal size. Similarly, in the above formula for estimating the numerical rate $r$ of convergence, we replace $u-u_{M}$ by $\tilde{u}_{2M}-u_{M}$ for $M=N,2N$.

Tables 1, 3 and 4 seem to indicate that the errors in the computed solutions obey the bound

[TABLE]

of equation (4.30). Thus Theorem 4.14 above appears not to be sharp; this weakness in the analysis was discussed in Remark 4.17.