---
abstract: |
    A numerical algorithm is presented to solve a benchmark problem proposed by Hemker[@hemker].  The algorithm incorporates asymptotic information into the design of appropriate piecewise-uniform Shishkin meshes.  Moreover, different co-ordinate systems are utilized due to the different geometries and associated layer structures that are involved in this problem.  Numerical results are presented to demonstrate the effectiveness of the proposed numerical algorithm. 

    [**Keywords:**]{} Singularly perturbed, Shishkin mesh, Hemker problem. 

    [**AMS subject classifications:**]{} 65N12, 65N15, 65N06. 
author:
- 'A. F. Hegarty [^1][^2]'
- 'E. O’Riordan [^3]'
title: '**A numerical algorithm to computationally solve the Hemker problem using Shishkin meshes**'
---

*Dedicated to the memory of P. W. Hemker, who inspired this research work*

Introduction
============

In [@hemker] Hemker proposed a model test problem in two space dimensions, defined on the unbounded domain $ \mathbb{R}^2 \setminus \{ x^2+y^2 \leq 1\} $, which is exterior to the unit circle.  The problem involves the simple constant coefficient linear singularly perturbed convection-diffusion equation $$\label{Piet's-problem}
-\ve \triangle  u + u_x=0, \quad \hbox{for} \quad x^2+y^2 >1;$$ with the boundary conditions $u(x,y) =1$, if $x^2+y^2 = 1$ and $ u(x,y) \rightarrow 0$ as $x^2+y^2 \rightarrow \infty$. An exponential boundary layer and characteristic interior layers appear in the solution of this problem.  In neighbourhoods of the two points $(0,\pm1)$, where the characteristic lines $y=C$ of the reduced problem ($ v_x=0$) are tangential to the circle, there are transition regions, where the steep gradients in the solution migrate from the exponential boundary layers (located on the left side of the unit disk) to the characteristic internal layers which are emanating from the characteristic points $(0,\pm1)$ (see Figure \[Fig.0\]).  To design a numerical method, which produces stable and accurate approximations over the entire domain, for arbitrary small values of the perturbation parameter $\ve$, is seen in [@john1; @bosco; @havik; @hemker; @volker1] as a reasonable challenge to the numerical analysis community.  Here we concentrate on the layers that are appear in the vicinity of the circle by considering the problem restricted to a finite domain.  Hence, we do not address the potential merging of the two characteristic layers, which can occur at a sufficentlly large distance of $O(1/{\ve})$ downwind (i.e., $x >> 1$) of the circle. 

In several publications (e.g.[@john1; @gabriel]) the Hemker problem is used to test the stability of numerical algorithms designed for a wide class of convection-dominated convection-diffusion problems, as classical finite element methods produce spurious oscillations for this type of singularly perturbed problems.  In [@tailor1], a stable numerical method was constructed on a quasi-uniform mesh for problem (\[Piet’s-problem\]), but there was a limited discussion of the accuracy of the numerical approximations.  In this paper, we guarantee parameter-uniform stability of the discrete operators by using simple upwinded finite difference operators.  Our main focus is on the design of an appropriate layer-adapted mesh, so that we guarantee that a significant proportion of the mesh point lie within the layers.  In Theorem \[main\], several pointwise bounds on the solution of the continuous problem are established, from which the location and widths of any layers are identified.  Asymptotic analysis [@eckhaus; @hemker; @ilin; @lagerstrom; @waechter] has also been used to determine the location and scale of all the layers that can appear in the solution of problem (\[Piet’s-problem\]). 

The numerical algorithm constructed in this paper is composed of several different Shishkin meshes [@gis1] defined across different co-ordinate systems aligned to three overlapping subdomains.  As we lack sufficient theoretical information about the localized character of the partial derivatives of the continuous solution, we have no meaningful pointwise bounds for the approximation error associated with any computational algorithm applied to the Hemker problem.  Here, we test for convergence of the numerical approximations using the double-mesh principle [@fhmos] and, more importantly, to identify when any numerical method fails to be convergent.  We emphasize that we shall estimate the global pointwise convergence of the interpolated computed approximations across the entire domain.  Numerical results with a preliminary version of this algorithm were reported in [@glasgow-BAIL]. 

The Shishkin mesh [@fhmos; @mos] is a central component of the algorithm.  The simplicity of this mesh is one of the key attributes of this particular layer-adapted mesh, which allows easy extensions to more complicated problems.  Shishkin meshes have the additional property that, if one has established parameter-uniform nodal accuracy [@fhmos] in a subdomain, then (for problems with regular exponential boundary layers [@hemkera] and characteristic boundary layers [@hemkerb]) this nodal accuracy extends to global accuracy across the subdomain using basic bi-linear interpolation.  This interpolation feature of the mesh permits us employ overlapping subdomains, with different co-ordinate systems aligned to the local geometry of the layers, and to subsequently test computationally for global accuracy across the entire domain. 

In §2, we identify bounds on the solution of the Hemker problem restricted to a bounded domain.  In §3, we discuss how we computationally estimate the order of parameter-uniform convergence for any numerical method.  In §4, we construct and describe the numerical algorithm, which involves four distinct stages.  The first two stages generate an initial approximation which has defects only near the characteristic points.  The third and fourth stage correct this initial approximation.  In §5, we present numerical results to illustrate the performance of the final algorithm.  The numerical results indicate that this new algorithm is generating numerical approximations which are converging, over an extensive range of the singular perturbation parameter, to the solution of a bounded-domain version of the original Hemker problem. 

[**Notation**]{}: We will employ three distinct co-ordinate systems in this paper.  A Cartesian co-ordinate system $(x,y)$, a polar co-ordinate system $(r,\theta)$ and a particular parabolic co-ordinate system $(s,t)$. In each co-ordinate system, we adopt the following notation for functions: $$f(x,y)=\hat f(s,t) =\tilde f(r,\theta).$$ We use these co-ordinate systems to solve various sub-problems on an annulus $ A$, a rectangle $S$ and planar regions $ Q^+,  Q^-$. The numerical solution is determined using four sub-components $U_A, U_B, U_C$ and $U_D$; where $U_A$ is defined over the annulus $A$, $U_C$ is defined over the planar region $ Q^+ \cup Q^-$ and $U_B, U_D$ are defined over the rectangular region $S$. In addition, the algorithm produces an initial global approximation $\bar U^N_1$, which is shown to lose accuracy in some of the layers.  This initial approximation is subsequently corrected to produce a globally pointwise accurate approximation $\bar U^N_2$. Throughout the paper, $\Vert \cdot \Vert _D$ denotes the pointwise norm measured over the domain $D$.

The continuous problem
======================

In this paper, we confine our discussion to the Hemker problem (\[Piet’s-problem\]) posed on a bounded domain of dimension $O(R^2)$. Consider the singularly perturbed elliptic problem: Find $ u (x,y)$ such that

\[cont-prob\] $$\begin{aligned}
 L u:=- \ve \triangle u +u_x =0, \quad (x,y) \in   \Omega ; \\
u=0, \  (x,y) \in \partial \Omega _O; \quad   u_x =0, \  (x,y) \in \partial \Omega _N;\quad
 u=1, \  (x,y) \in \partial \Omega _I; \end{aligned}$$ where the bounded domain $\Omega$ and the boundaries $\partial \Omega _O, \partial \Omega _N, \partial \Omega _I$ are defined to be $$\begin{aligned}
\Omega := \{(x,y) \vert 1 < x^2 +y^2 < R, x \leq 0\} \cup \{(x,y) \vert -R \leq y \leq R, 0 \leq x \leq R, 1 < x^2 +y^2 \},
\\
\partial \Omega _N:= \{ (R,y) |-R < y < R \}, \quad
  \partial \Omega _I:= \{(x,y) \vert x^2 +y^2 =1 \}; \\
\partial \Omega _O:= \{ (x,y) \vert x^2 +y^2 = R, x \leq 0\} \cup \{(x,y) \vert y =\pm R, 0 \leq x \leq R \}.\end{aligned}$$

A sample computed solution (using the algorithm (\[corrected-approx\])) is displayed in Figure \[Fig.0\], which illustrates the location of the layers that can appear in the solution.  In all of our numerical experiments, we have simply taken $R=4$.

We have a minimum principle associated with this problem. 

[@prot pg.  61] For any $ w \in C^0(\bar D) \cap C^2(D), D \subset  \Omega$, where $ L w (x,y) \geq 0, (x,y) \in D$, then $\min _{\bar D}  w \geq \min _{\partial D} w$.

This follows from observing that if $ v := e^{-\frac{x}{2\ve} } w$, then $ L  w = e^{\frac{x}{2\ve} }(-\ve \triangle v +\frac{1}{4\ve}  v)$ which will lead to a contradiction if $ w$ had an internal minimum within $D$.

Thus, we have a comparison principle. 

If $ w ,v \in C^0(\bar D) \cap C^2(D), D \subset  \Omega$ are such that $ L  w (x,y)\geq  L v(x,y), \ \forall (x,y) \in  D$ and $ w \geq  v$ on the boundary $\partial D$, then $w (x,y) \geq v(x,y),\ \forall (x,y) \in \bar D$.

For any open connected subdomain $D \subset \Omega$ we define the boundaries $$\partial D _N:=  \{ \vec p \in \partial D \vert  (1,0)\cdot \vec n_p >0\} \quad
\hbox{and} \quad 
\partial D _O:= \partial D \setminus \partial D _N;$$ where $\vec n_p$ is the outward normal to $\partial D$ at $\vec p$. As in [@prot pg.65], we can establish that

\[compar2\] If $ w , v \in C^0(\bar D ) \cap C^2(D ), D \subset  \Omega$ are such that $ L w (x,y)\geq  Lv(x,y), \ \forall (x,y) \in  D$; $ w \geq  v$ on the boundary $\partial D _O  $ and $ w_x \geq  v_x$ on the boundary $\partial D _N  $, then $ w (x,y) \geq  v(x,y),\ \forall (x,y) \in \bar D$.

Using these results, we can establish the following bounds on the solution:

\[main\] Assuming $\ve$ is sufficiently small, then the solution $u$ of problem (\[cont-prob\]) satisfies the following bounds

\[solutions-bounds\] $$\begin{aligned}
0 \leq  u(x,y) \leq 1, && (x,y) \in \bar \Omega ; \label{basic}\\
\tilde u(r,\theta) \leq Ce^{\frac{ \kappa \cos (\theta ) (r-1)}{\ve}}, && \cos \theta \leq -\ve ^{1/3},  \  \kappa <1 ;
\label{left}\\
 u(x,y) \leq Ce^{-\frac{(0.5x^2+\vert y \vert -1)}{3\ve ^{2/3}}}, &&   x \in \ve ^{1/3}[-1,1]; \label{center} \\
 u (x,y)\leq Ce^{-\frac{\mu (y-1)}{\sqrt{\ve}}}, \ y \geq 1; &&   u(x,y) \leq Ce^{-\frac{\mu (y+1)}{\sqrt{\ve}}}, \ y \leq -1,\ x\geq 0, \quad \mu >0.   \label{right}\end{aligned}$$

The first bound follows easily from the minimum principle.  See details in the appendix for all of the remaining bounds. 

Computationally testing for convergence
=======================================

As identified in [@hemker], we do not have a useable closed form representation of the exact solution to the Hemker problem to allow us evaluate the accuracy of any computed approximation.  The infinite series representation [@hemker] for the exact solution has difficulties for moderately small values of the singular perturbation parameter $\ve$. Hence, to test for convergence we rely on the double-mesh method of estimating the order of convergence [@fhmos Chapter 8].  We elaborate on this experimental approach in this section. 

For every particular value of $\ve$ and $N$, let $U_\ve^N$ be the computed solutions on certain meshes $\Omega _\ve ^N$, where $N$ denotes the number of mesh elements used in each co-ordinate direction.  Define the maximum local two-mesh global differences $ D^N_\ve$ and the maximum parameter-uniform two-mesh global differences $D^N$ by $$D^N_\ve:= \Vert \bar U_\ve^N-\bar U_\ve ^{2N}\Vert _{\Omega _\ve^N \cup \Omega _\ve^{2N}} \quad \hbox{and}  \quad D^N:= \sup _{0 < \ve \leq 1} D^N_\ve,$$ where $\bar U_\ve ^N$ denotes the bilinear interpolation of the discrete solution $U_\ve ^N$ on the mesh $\Omega _\ve ^{N}$. Then, for any particular value of $\ve $ and $N$, the local orders of global convergence are denoted by $\bar p ^{N}_\ve$ and, for any particular value of $N$ and [*all values of $\ve $*]{}, the [**parameter-uniform**]{} global orders of convergence $\bar p ^N$ are defined, respectively, by $$\bar p^N_\ve:=  \log_2\left (\frac{D^N_\ve}{D^{2N}_\ve} \right) \quad \hbox{and} \quad  \bar p^N :=  \log_2\left (\frac{D^N}{D^{2N}} \right).$$ If, for a certain class ${\cal{C}}$ of singularly perturbed problems, there exists a theoretical error bound of the form: There exists a constant $C$ independent of $\ve$ and $N$ such that for all $\ve >0$ $$\label{theory}
\Vert \bar U^N -u \Vert _{\Omega}\leq CN^{-p},\quad p >0;$$ then it follows that $$D^N \leq C(1+2^{-p}) N^{-p}.$$ Hence, for any particular sample problem from this class ${\cal{C}}$, we expect to observe this theoretical convergence rate $p$ in the computed rates of convergence $\bar p^N $.  That is, we expect that $\bar p^N \approx p$.

A useful attribute of the two-mesh method is that we can use it to identify when a numerical method is not parameter-uniform.  Observe that $$\Vert \bar U_\ve ^N - \bar U_\ve ^{2N} \Vert _{\Omega} \leq \Vert \bar U_\ve^N - u \Vert_{\Omega} +\Vert u - \bar U_\ve ^{2N} \Vert_{\Omega}.$$ Hence, if the parameter-uniform two mesh differences $D^N$ fail to converge to zero, then the numerical method is also not a parameter-uniform numerical method.  In our quest for a parameter-uniform numerical method, we used this feature to identify necessary components to construct a parameter-uniform numerical method.  On the other hand, without the existence of a theoretical error bound (\[theory\]) (as is the case with the Hemker problem), if the global two mesh differences $D^N$ are seen to converge then we can only conclude that the numerical method may be a parameter-uniform numerical method.  We would require a theoretical parameter-uniform error bound on the numerical approximations, before one can assert that the numerical method is indeed parameter-uniform. 

For any numerical method applied to a class of singularly perturbed problems, our primary interest is in determining the parameter-uniform orders of global convergence $\bar p^N$. However, we can also examine the local orders of convergence $\bar p_\ve^N$ to see how the numerical method performs for each possible value of $\ve$ over the range $0<\ve \leq 1$. In general, we note that $\bar p^N \neq \min _{\ve}  \bar p^N_\ve$. In the case of piecewise-uniform meshes, certain anomalies can sometimes be observed in the local orders of convergence (i.e., $\bar p_\ve ^N \not\approx \bar p^N$) for certain values of $\ve$. We illustrate this effect with the following theoretical example.  Based on the nature of the typical errors on a piecewise-uniform Shishkin mesh in a one dimensional convection-diffusion problem, suppose that the two-mesh differences $D_\ve^N$ were of the form $$D_\ve ^N := \Bigl\{ \begin{array}{ll} \frac{\rho}{1+\rho} \qquad   \hbox{if} \quad k\ve \ln N \geq 1 \\
N^{-1} \qquad \hbox{if} \quad k\ve \ln N < 1
\end{array} ; \qquad  \hbox{where} \quad \rho := \frac{1}{\ve N}, \quad k \geq 1 .$$For this theoretical example, the parameter-uniform two mesh differences can be explicitly determined.  Note first that $$D_\ve ^N  \leq \Bigl\{ \begin{array}{ll} \frac{k\ln N}{N+k\ln N} \qquad \qquad  \hbox{if} \quad k\ve \ln N \geq 1 \\
N^{-1} \qquad \qquad  \hbox{if} \quad k\ve \ln N < 1
\end{array} .$$ Hence, if $N\geq 4, k \geq 1$, we have $$D^N =\frac{k\ln N}{N+k\ln N} \quad \hbox{and} \quad \lim _{N\rightarrow \infty} \bar p^N =1.$$ Let us now consider the particular values of $\ve = 2^{-4},k=4$ and $N=32$. First observe that $$D_{\ve =2^{-4}} ^N  = \Bigl\{ \begin{array}{ll} \frac{\rho _*}{1+\rho _*} \qquad \qquad  \hbox{if} \quad  \ln N \geq 4 \\
N^{-1}\qquad \qquad  \hbox{if} \quad   \ln N < 4
\end{array} , \qquad \hbox{where} \quad \rho _* = 16N^{-1}.$$ Then, in particular, $$D_{\ve =2^{-4}} ^{32} = 2^{-5}, D_{\ve =2^{-4}} ^{64} =0.2,$$ which yields the order of convergence as $$\bar p_{\ve =2^{-4}} ^{32} \approx -2.68 \quad \hbox{although}\quad  \bar p ^{32} \approx 0.55.$$ Thus, we can have negative local orders of convergence $\bar p^N_\ve$, for particular values of $\ve $ and $N$ and still have positive parameter-uniform orders of convergence $\bar p^N$. This phenomena will appear in the numerical results section in §5. 

In practice, note that the parameter-uniform orders $\bar p^N$ can only be estimated over a finite set $R^J_\ve := \{2^{-j}, j=0,1,\ldots J\}$ of values of the singular perturbation parameter $\ve \in (0,1]$. That is, we define $$\bar p^N_{R^J_\ve} :=  \log_2\left (\frac{D_{R^J_\ve}^N}{D_{R^J_\ve}^{2N}} \right), \quad \hbox{where} \quad  D_{R^J_\ve}^N:= \max _{\ve \in R^J_\ve } D^N_\ve.$$ When a method is known to be parameter-uniform, $J$ is taken sufficiently large so $\bar p^N_{R^K_\ve} = \bar p^N_{R^J_\ve}$, for any $K >J$ and $\bar p^N_{R^J_\ve}$ is taken to be the computed estimate of $\bar p^N$.

In this paper, we construct a numerical method that displays positive orders of convergence $\bar p^N_{R^{20}_\ve}$, for $N$ sufficiently large; i.e., $N \geq N_0$, where $N_0$ is independent of $\ve$, when the numerical method is applied to the Hemker problem and the range of the singular perturbation parameter is $\ve \in R_\ve:= \{2^{-j}, j=0,1,\ldots 20\}$. For smaller values of the parameter ($\ve < 2^{-24}$), we have observed a degradation in the local orders of convergence.  Hence, we cannot claim that the numerical method described in this paper is parameter-uniform. 

To conclude this section, we note that the two-mesh differences are computed to enable the computation of approximate orders of convergence.  To generate approximations to the pointwise errors for any particular value of $N$ and $\ve$, we compare the computed solution for a given number of mesh points $N$ to the computed solution on a fine mesh.  That is, we approximate the nodal error by $$\Vert U_\ve^N-u \Vert _{\Omega ^N_\ve } \approx \Vert  U_\ve^N-\bar U_\ve^{4N_*} \Vert _{\Omega _\ve^N }, \quad \hbox{for any} \quad N \leq N_*;$$ and the global error by $$\Vert \bar U_\ve^N-u \Vert _{\Omega } \approx \Vert \bar U_\ve^N-\bar U_\ve^{4N_*} \Vert _{\Omega _\ve^N \cup  \Omega _\ve^{4N_*}}, \quad \hbox{for any} \quad N \leq N_*.$$

The numerical algorithm
=======================

Polar coordinates are a natural co-ordinate system to employ in the semi-annular region to the left of the line $x=0$ and rectangular co-ordinates are natural to the right of the line $x=0$. To incorporate both co-ordinate systems, we first generate an approximate solution to the problem (\[cont-prob\]) on the sector $$\label{annulus}
\tilde A:= \{(r,\theta) \vert 1< r < R, \frac{\pi}{2} - \tau\leq  \theta  \leq\frac{3\pi}{2} +\tau\}$$ which is a proper subset of the domain $\Omega$ and the parameter $\tau$ is specified below in (\[tau-1\]). 

The continuous problem (\[cont-prob\]), restricted to the sector $\tilde A$, is transformed into the problem: Find a periodic function, $\tilde u(r,\theta) =  u(x,y)$ such that

\[polar\] $$\begin{aligned}
\tilde  L\tilde u:= -\frac{\ve }{r^2}\tilde  u _{\theta, \theta}  -\ve  \tilde u_{rr} +  \bigl(\cos (\theta ) -{\frac{\ve}{r}} ) \tilde u _r  - \frac{\sin (\theta )}{r}  \tilde  u _{\theta} = 0,  \ \hbox{in } \tilde A ;\\
\tilde u(1,\theta) =1, \quad ; \quad  \tilde  u(R,\theta) =0, \quad \frac{\pi}{2} -\tau \leq \theta  \leq \frac{3\pi}{2} +\tau.\end{aligned}$$ The remaining boundary points of this sector $\tilde A$ are internal points within the domain $\Omega$, where the discrete solution has not yet been specified.  In order to generate an initial approximation to the solution $\tilde  u$, we impose homogeneous Neumann conditions at these internal points of the form: $$\begin{aligned}
 u^A_x (x,y)=0, \quad \hbox{if} \quad x^2+y^2=R^2 \quad \hbox{and} \quad x >0 \\
 \tilde u^A_\theta (r,\frac{\pi}{2}-\tau  )= \tilde u^A_\theta (r,\frac{3\pi}{2}+\tau  )=0, \quad \hbox{for} \quad 1 < r < R,\end{aligned}$$

where $u^A(x,y) \approx u(x,y)$ for $x \leq 0$.

The choice of a homogeneous Neumann condition at the outflow is motivated by the following observation: Consider the one-dimensional convection-diffusion problem: Find $z(x), x \in [0,L]$ such that $$-\ve z'' +az'=f(x), \ x \in (0,L);\quad  z(0)=A, z(L)=B; \quad a(x) \geq \alpha >0$$ and the approximate problem: Find $z_A(x), x \in [0,L]$ such that $$-\ve z_A'' +az_A'=f(x), \ x \in (0,L);\quad  z_A(0)=z(0), z'_A(L)=0.$$ Using a comparison principle (as in Theorem \[compar2\]), we can establish the bound $$\vert (z-z_A)(x) \vert \leq \frac{\ve \vert z'(L)\vert}{\alpha } e^{-\frac{\alpha(L-x)}{\ve}}.$$ Hence, $z_A$ is an $O(N^{-1})$-approximation to $z$, at some $O(\ve \ln N)$ distance away from the end-point $x=L$. That is: $$\vert (z-z_A)(x) \vert \leq C\ve \vert z'(L)\vert N^{-1}\leq CN^{-1}; \quad \hbox{if} \quad x \in [0, L-\frac{\ve \ln N}{\alpha}].$$

This problem (\[polar\]) is discretized using simple upwinding on a tensor product piecewise-uniform Shishkin mesh, whose construction is motivated by the bounds (\[left\]), (\[center\]).  Two transition points (where the mesh step changes in magnitude) are used in the radial direction.  The choice of the first transition point $\sigma _1$ is motivated by considering the bound $(\ref{left})$ at some fixed distance $x < - \delta <0$ to the left of $x=0$ and by the theoretical error bounds [@annulus] established for this mesh in the region where $x < - \delta <0$. The choice of this point $x=-\delta$ is arbitrary.  As in [@annulus], we simply take it to correspond to a criticial angle $\theta _*$, such that $$\kappa \cos \theta _* = -\frac{1}{2}.$$ Hence.  we have that $$\tilde u (r,\theta) \leq C e^{-\frac{r-1}{2\ve}}, \quad \hbox{for} \quad \theta _* \leq \theta \leq 2\pi - \theta _*$$ and $$\tilde u (r,\theta) \leq CN^{-1}, \quad \hbox{if} \quad r \geq 1+2\ve \ln N \quad \hbox{and} \quad \theta _* \leq \theta \leq 2\pi - \theta _*.$$ A second transition point $\sigma _2$ is motivated by the bound $(\ref{center})$ applied along the line $x=0$. In the angular direction, the bound $(\ref{center})$ also motivates the inclusion of a transition point $\tau$ in the vicinity of the characteristic points.  See also [@eckhaus pg.  269], [@ilin pg.  188] and [@waechter pg.  1183] for motivation for these scales in the vicinity of the characteristic points $(0,\pm 1)$.

[**The Shishkin mesh $\tilde  \Omega ^{N}_A$**]{}

*The radial domain $[1,R]$ is divided into three subregions.  The radii $
r=1, \ r=1+\sigma _1, \ r=1+ \sigma _1 +\sigma _2, \  r=R,
$ mark the subregion boundaries and the radial transition points are taken to be*

\[fitted-mesh-special\] $$\sigma _1:= \min \{ \frac{R-1}{4},  2\ve  \ln N   \} \quad \hbox{and} \quad
 \sigma _2:= \min \{ \frac{R-1}{4}, 3 \ve ^{2/3} \ln N  \}.$$ The $N$ radial mesh points are distributed in the ratio $N/4:N/4:N/2$ across these three subintervals.  For the angular coordinate, the interval $[\frac{\pi}{2}-\tau,\frac{3\pi}{2}+\tau ]$ is split into three subintervals with the start/end points of each subinterval, respectively, at $$\frac{\pi}{2} -\tau, \frac{\pi}{2} +\tau, \frac{3\pi}{2} -\tau, \frac{3\pi}{2} +\tau$$ and the mesh points are distributed in the ratio $N/4:N/2:N/4$ across the three associated subintervals.  The transition points are determined by $$\label{tau-1}
 \tau := \min \{ \frac{\pi}{6}, \sqrt{6} \ve^{1/3} \ln N   \}.$$

A schematic image of this mesh is presented in Figure \[Fig.1\].  However, in practice, the refinement in the radial direction only becomes apparent to the user for very small values of $\varepsilon$.

At the mesh points on the sector $ \tilde A$, the computed solution will (in polar coordinates) be denoted by $U_A(r_i,\theta_j), \ (r_i,\theta_j) \in \bar \Omega _A^{N}$. This approximation is extended to the global approximation $\bar U_A(r,\theta),\ (r,\theta) \in \bar A$, using simple bilinear interpolation[^4]

[**Stage 1**]{}: [*Numerical method for $U_A$ defined over an annular subregion $\tilde A$*]{}

Find $U_A$ such that [^5]:

\[discrete-problem-A\] $$\begin{aligned}
 \tilde L^{N} U_A=  0,\qquad (r_i,\theta_j) \in  \Omega ^{N}_A; \\ \hbox{where} \quad 
 \tilde L^{N} U:= -\frac{\ve}{r_i^2}  \delta ^2 _{\theta}U   -\ve  \delta ^2 _{r}U  +  (\cos (\theta _j )-\frac{\ve}{r_i} ) D^{\pm}_r U  -\frac{\sin (\theta _j )}{r_i}D^{\pm}_\theta U; 
\\
U_A(1,\theta _j) =1, \qquad U_A(R,\theta _j) =0,  \quad \frac{\pi}{2} \leq \theta _j \leq \frac{3\pi}{2}; \\
(\cos \theta  D^-_r-\frac{\sin \theta}{R}  D^\pm_\theta \bigr) U_A(R,\theta _j) =0,  \quad   \theta _j \in (\frac{\pi}{2} -\tau , \frac{\pi}{2}) \cup  (\frac{3\pi}{2}, \frac{3\pi}{2} +\tau );\\
D^+_\theta U^A (r_i,\frac{\pi}{2}-\tau  )= D^-_\theta U^A (r_i,\frac{3\pi}{2}+\tau )=0, \quad \hbox{for} \quad 1 < r_i < R.\end{aligned}$$

The boundary condition (\[discrete-problem-A\]d) corresponds to applying the Neumann condition $u_x=0$ at these internal points of $\bar \Omega$.

In Table \[Tab1\], we present the results from applying this numerical method to problem (\[cont-prob\]) posed on the sector $A \subset \Omega$. We display the orders of convergence only for the region where $x \leq 0$ and we observe global convergence over the parameter range $\ve \in [2^{-20},1]$.

[|c| c c c c c c c|]{}\
$\ve | N$&N=**[8]{}&**[16]{}&**[32]{}&**[64]{}&**[128]{}&**[256]{}&**[512]{}\
$\ve =1$ & 0.9855 & 0.9966& 0.9996 & 0.9996& 1.0002 & 1.0001 & 1.0000\
$2^{-2}$ & 1.0014 & 0.9115 & 0.9643& 0.9867 & 0.9918 & 0.9958 & 0.9979\
$2^{-4}$ & 0.4693 & 0.6789 & 0.7825 & 0.6875 & 0.6950& 0.7541 & 0.9911\
$2^{-6}$ & 0.5832 & 0.7176& 0.7760 & 0.6887 & 0.7691 & 0.7990 & 0.8255\
$2^{-8}$ & 0.7205& 0.7771 & 0.8441 & 0.6893 & 0.7713& 0.7991& 0.8266\
$2^{-10}$ & 0.3249 & 0.7079& 0.9677& 0.8376 & 0.9761 & 0.9231 & 0.8270\
$2^{-12}$ & 0.1131& 0.4659& 0.8732 & 0.9145& 0.8086& 0.9597& 1.0268\
$2^{-14}$ & 0.1242 & 0.5344 & 0.6122& 0.7242 & 0.8605& 0.7995& 0.9612\
$2^{-16}$ & 0.1669 & 0.4506 & 0.4704 & 0.6033 & 0.7562& 0.8573 & 0.8242\
$2^{-18}$ & 0.2113& 0.2972 & 0.3714& 0.5046 & 0.6604& 0.7978& 0.8561\
$2^{-20}$ & 0.2111& 0.1873& 0.2770 & 0.4167& 0.5625& 0.7077 & 0.8373\
$\bar p^N_{R^{20}_\ve}$& 0.2111& 0.1873& 0.2770 & 0.4167& 0.5625& 0.7077 & 0.8373\
**************

We next introduce a rectangular mesh, which will be aligned to the internal characteristic layers.  We retain the computed solution $\bar U_A$ in the upwind region where $x \leq 0$ and then solve the problem (\[cont-prob\]) over the remaining rectangle $$\label{rectangle}
S:= \{ (x,y) | 0 < x \leq R ,-R \leq y \leq R\}$$ using a piecewise-uniform mesh $ \Omega ^N_S$, whose transition parameters are related to the bounds (\[right\]) on the continuous solution $u$. Note that there is a free parameter $\mu$ in these bounds, which we simply set to be $\mu =1/2$. Then, $$u(x,y) \leq CN^{-1}, \quad \hbox{if} \quad y \geq 1 +2\ve \ln N.$$

[**The Shishkin mesh $ \Omega ^{N}_S$**]{}

*The mesh $ \Omega ^N_S:= \omega _u \times \omega _3$ is a tensor product mesh of a uniform mesh $\omega _u $ in the horizontal direction and a Shishkin mesh $\omega _3$, which refines in the region of the interior characteristic layers.  The mesh $\omega _3$ is generated by splitting the vertical interval $[-R,R]$ into the five subregions*

\[Shish-par-mesh\] $$[-R,-1 -\tau_2 ] \cup [-1 -\tau _2, -1 +\tau _1] \cup [-1 +\tau _1, 1-\tau _1] \cup [1-\tau_1, 1+\tau _2] \cup [1+\tau _2, R ],$$ distributing the mesh elements in the ratio $N/8:N/4:N/4:N/4:N/8$ and $$\tau _1  := \min \{ \frac{1}{2},   2\ve ^{1/2} \ln N  \}; \quad 
\tau _2 :=  \min \{ \frac{R-1}{2}, 2 \ve ^{1/2}  \ln N  \} .$$

[**Stage 2**]{}: [*Numerical method for $U_B$ defined over the downwind region $S$:*]{}

Find $U_B(x_i,y_j)$ such that

\[discrete-problem-B\] $$\begin{aligned}
L^{N,M} U_B:= \Bigl( -\ve  \delta ^2 _{x}   -\ve  \delta ^2 _{y}  + D^{-}_x \Bigr)  U_B(x_i,y_j)     = 0 ,  \quad (x_i,y_j) \in \Omega ^N_S \setminus \bar C_1;\\
U_B(x_i,y_j)\equiv 1, \quad (x_i,y_j) \in \Omega ^N_S \cap  \bar C_1; \quad C_1:=\{ (x,y) | x^2+y^2 \leq 1\}; 
 \end{aligned}$$ with the remaining boundary values computed from the equations $$\begin{aligned}
D^-_x  U_B(R,y_j)=0,  \ -R < y_j < R,\quad U_B(x_i,-R)= U_B(x_i,R)=0, \ x_i \in [0,R]; \\ U_B(0,y_j) = \bar U_A (0,y_j),\quad  y_j \in (-R, R) \setminus [-1,1].\end{aligned}$$

Here $\bar U_A (0,y_j)$ is a linear interpolant of the values $U_A (r_i,\theta _j)$ along the line $x=0$.

The initial computed global approximation $\bar U_1^N$ to the solution of problem (\[cont-prob\]) is: $$\label{initial-approx}
\bar U_1^N (x,y) := \left \{
\begin{array}{ll}
\bar U_A(r,\theta ), & \hbox{for} \quad (r,\theta ) \in \bar \Omega _A \setminus (\{ x\geq 0 \}) \\
\bar U_B (x,y), & \hbox{for} \quad (x,y) \in \bar S \setminus \bar C_1, 
\end{array}\right.$$ where $U_A$ is defined by (\[discrete-problem-A\]) and $U_B$ is defined by (\[discrete-problem-B\]). 

In Table \[Tab2\], we do not observe convergence of these initial numerical approximations $\bar U_1^N$. Hence, although we observe convergence in the annulus up to $x \leq 0$, this does not suffice to generate convergence across the entire domain, even if the mesh is fitted to the characteristic layers.  In Figures \[Spike-N128\], \[Spike-N256\] we plot the error across the entire domain and we observe a spike in the global pointwise error in the vicinity of the characteristic points.  This global error does not decrease when the mesh is refined. 

We next describe the construction of a correction $\bar U_2^N$ to the initial approximation $\bar U_1^N$.

[|c| c c c c c c c|]{}\
$N$&**[8]{}&**[16]{}&**[32]{}&**[64]{}&**[128]{}&**[256]{}&**[512]{}\
$\ve =1$ & 1.1831 & 0.5014& 0.7157 & 0.7719& 0.8259& 0.8743& 0.9100\
$2^{-2}$ & 1.7592 & 0.0305 & 0.4492 & 0.6880 & 0.7791 & 0.8460 &0.8944\
$2^{-4}$ & 1.0764 & 1.0413 & 0.0241 & 0.4368 & 0.6719 & 0.7768 & 0.8486\
$2^{-6}$ & 0.2376 & 0.7503 & 1.7967 & -0.9041 & 0.7151& 0.4173 & 0.7648\
$2^{-8}$ & -0.0726 & 0.3699 & 0.6842 & 0.1372 & 0.8439 & -0.0165 & 0.6041\
$2^{-10}$ & -0.0700 & 0.0504 & 0.4604 & 1.0622 & 0.4357 & 0.2011 & 0.4706\
$2^{-12}$ & -0.1139& -0.0051 & 0.2692 & 0.8506 & 1.0414 & -0.0081 & 0.3631\
$2^{-14}$ &-0.1508 & -0.0405 & 0.1558 & 0.5954 & 1.2268 & 0.1912 & 0.2453\
$2^{-16}$& -0.1680 &-0.0685 & 0.0477 & 0.3999 & 1.0007 & 0.7912 & 0.0780\
$2^{-18}$ & -0.1686 &-0.0873& -0.0159 & 0.2119 & 0.7662 & 1.3778 & -0.0380\
$2^{-20}$ & -0.1690 & -0.0922 & -0.0467 & 0.0795 & 0.5156 & 1.2055 & 0.5538\
**************

![ Approximate global error $ \bar U_1^{128}-\bar U_1^{2048}$ on the mesh $\tilde \Omega _A^{N} \cup  \Omega ^N_S$, (\[fitted-mesh-special\]) and (\[Shish-par-mesh\]), for $\ve =2^{-10}$[]{data-label="Spike-N128"}](Hemker_fig4)

![ Approximate global error $ \bar U_1^{256}-\bar U_1^{2048}$ on the mesh $\tilde \Omega _A^{N} \cup  \Omega ^N_S$, (\[fitted-mesh-special\]) and (\[Shish-par-mesh\]), for $\ve =2^{-10}$[]{data-label="Spike-N256"}](Hemker_fig5)

Let us return to the bounds on the continuous solution given in Theorem \[main\].  From (\[center\]), the bound on the solution remains constant along the parabolic path $y-1=-\frac{1}{2} x^2$. This is the motivation to introduce a third co-ordinate system that is aligned to these parabolic curves.  Under this transformation, a mixed derivative term will appear in the transformed elliptic operator and we are required to restrict the dimensions of the sub-domain (where this new transformation is utilized) in order to preserve inverse-monotonicity of the corresponding discrete operator. 

We introduce a patched region $Q:=Q^+\cup Q^-$, in a neighborhood of the vertical line $x=0$. We discuss the approach on the upper region $$\label{defQ}
Q^+:= \{ (x,y) \vert y=t-x^2/2,\quad  0< x < L,\ 1-3\tau _0  < t < 1+3\delta \};$$ with an analogous definition of the lower region $Q^-$. The width $L$ and height $3(\delta +\tau_0)$ of this strip will be specified later in order to retain stability of the discrete operator. 

A natural coordinate system for this patched region $Q^+$ is $$s=x, \quad t= y+x^2/2, \qquad (y = t-s^2/2).$$ Then let $\hat u(s,t) := u(x,y),  \hat Q^+ := (0,L) \times (1-\tau _0 ,1+3\delta )$ and $$\begin{aligned}
u_x = \hat u _s +s \hat u_t, \ u_y= \hat u _t, \ u_{yy}=\hat u_{tt}, \quad
u_{xx} = \hat u_{ss} +2s \hat u_{ts} + s^2 \hat u_{tt} + \hat u_t.\end{aligned}$$ Under this transformation the problem (\[cont-prob\]) on the patched region can be specified as follows $$\begin{aligned}
\hat L \hat u (s,t) := -\ve (\hat u_{ss} +2s \hat u_{ts} + (1+s^2) \hat u_{tt}) + \hat u _s + (s -\ve)\hat u _t = 0, (s,t) \in \hat Q^+ \setminus \bar C_1\\
\hat u (s,t) = \hat u(s,t),\ (s,t) \in \overline{ \hat Q^+} \setminus \hat Q^+, \qquad  \hat u (s,t) =1, \quad \hbox{if} \quad s^2 + (t-s^2/2)^2 \leq 1.\end{aligned}$$ We will use the initial computed approximation $U_1^N$ (\[initial-approx\]) as boundary values to solve on the patched region $\hat Q^+$ as follows: $$\begin{aligned}
\hat L \hat u := -\ve (\hat u_{ss} +2s \hat u_{ts} + (1+s^2) \hat u_{tt}) +  \hat u _s + ( s -\ve)\hat u _t = 0,\quad  (s,t) \in \hat Q^+, \\
\hat u(s,t) =1,\quad \hbox{if} \quad s^2 +(t-s^2/2)^2 \leq 1, \qquad  \hat u(s,1+3\delta) = \bar U_1(s,1+3\delta), \ s \in [0,L], \\
\hat u(0,t) =\bar U_1(0,t),\qquad  \hat u_s (L,t) = 0,\ t \in (1-\tau _1 ,1+3\delta ).  \end{aligned}$$ We again note the use of a Neumann boundary condition at the artificial internal boundary $s=L$. To ensure that the corner point $(L, 1-\tau_0)$ of this patch $Q^+$ lies within the inner boundary ($ r <1$) of problem (\[cont-prob\]), we require that $$1-\tau _0 < L^2 +(1-\tau _0 - \frac{L^2}{2})^2;$$ which, in turn, requires that $$\tau _0 \geq \frac{L^2}{2}\frac{\sqrt{1+L^2}}{1+\sqrt{1+L^2}}.$$ For simplicity we shall simply take $$\tau _0 = \tau _3 \geq \frac{L^2}{2}.$$

We now specify the next phase of our numerical algorithm, where we correct the initial approximation $\bar U_1$.

We numerically solve the problem on the patch $\hat Q^+$, in the $s,t$ coordinate system, using a Shishkin mesh $ \hat \Omega ^{N,M}_P$ with two transition points located at $t=1+\tau _3$ and $1 + \tau_3 +\tau _4$. The choice for the transition parameters $\tau _3, \tau _4$ is motivated by the bounds (\[center\]) and (\[right\]), written in the $(s,t)$ coordinates combined yield $$\hat u (s,t) \leq C \min \{ e^{-\frac{t-1}{3\ve ^{2/3}}},  e^{-\frac{t-1}{\sqrt{\ve}}} \}, \quad \hbox{for} \quad s \in \ve ^{1/3} [-1,1].$$

[**The Shishkin mesh $\hat \Omega ^{N,M}_P$**]{} [ *We use a uniform mesh in the $s$-direction and a Shishkin mesh in the $t$-direction.  The vertical strip $[1-\tau _3 ,1+3\delta]$ is split into the following four sub-regions $$[1-\tau _3,1] \cup [1,1+\tau_3] \cup [1+\tau_3, 1+(\tau_3+\tau_4)] \cup [1+(\tau_3+\tau_4),1+3\delta]$$ where $$\label{tau-3}
\tau _3 := \min \{ \delta , 3{\ve} ^{2/3} \ln M \}, \qquad \tau _4 := \min \{ \delta , 2\sqrt{\ve}\ln M \},$$ and $M/4$ mesh elements are distributed uniformly within each of these sub-intervals.* ]{} 

Within this patched region, we adopt the following notation $$h = N^{-1}L, \quad k \leq k_j:= t_j-t_{j-1} \leq K; \quad 2\bar k_j = k_j+k_{j+1},$$ where, as $\tau _3 \leq \tau _4$, we have $$4M^{-1}  \tau_3  \leq k \leq K \leq   12 \delta M^{-1}.$$ Note that if we assume that $$\label{M-and-N}
 12\delta \leq MN^{-1},$$ then the maximum mesh step $K$ in the vertical direction will be such that $K \leq N^{-1}$.

[**Stage 3**]{}: [*Numerical method for $U_C$ defined near the characteristic points*]{}

Find $U_C$ such that:

\[discrete-problem-transformed\] $$\begin{aligned}
\hat L^{N,M}  \hat U_C:=
-\ve {\cal{L} }^{N,M}  \hat U_C +   D^{-}_s \hat U_C +(s_i-\ve) D^{\pm}_t \hat U_C   =0,\quad (s_i,t_j) \in  \hat \Omega ^{N,M}_P; \\
\hbox{where} \quad
  {\cal{L}} ^{N,M}  Y := \delta ^2_{ss} Y + 2s\delta _{st} Y +(1+s^2)\delta ^2_{tt} Y;\quad  2\delta _{st} := D^-_tD^-_s +D^+_tD^+_s 
    \label{patch-diff-op}\end{aligned}$$ and for the remaining mesh points $$\begin{aligned}
\hat U_C(s_i,t_j)=1,\quad \hbox{if} \quad s_i^2 +(t_j-s_i^2/2)^2 \leq 1; \\ \hat U_C(s_i,1+5\delta ) = \bar U_1(s_i,1+5\delta ), \ 0 \leq s_i \leq L;\\
\hat U_C(0,t_j)=\bar U_1(0,t_j), \quad D^-_s \hat U_C (L,t_j) = 0, \quad  1-3\delta< t_j < 1+3\delta .\end{aligned}$$

The presence of the mixed derivative term $\ve s \hat u _{st}$ in the transformed problem, creates the danger of loss of stability in the discretization of the differential operator $\hat L$ [@ray09].  However, by restricting the dimensions of this parabolic patch, we are able to preserve an appropriate sign pattern in the system matrix elements, so that the matrix ${\cal{L} }^{N,M}$ is an $M$-matrix. 

If we choose the dimension of the patch $\hat Q^+$ (\[defQ\]) to satisfy $$\label{def-L}
L \leq L_*:= 2\sqrt{NM^{-1}\tau _3} \quad \hbox{and} \quad 12 \delta \leq MN^{-1}$$ then the finite difference operator $\hat L^{N,M}$ (\[discrete-problem-transformed\]) satisfies a discrete comparison principle.  That is, for any mesh function $Z$, $$\begin{aligned}
\hbox{if} \quad Z(s_i,1-\tau _1) \geq 0, Z(s_i,1+3\delta) \geq 0, Z(0, t_j) \geq 0 \quad D^-_sZ(L, t_j) \geq 0, \ \forall (s_i,t_j) \in \bar Q^+ \\
\hbox{and} \quad \hat L^{N,M}Z(s_i, t_j) \geq 0, \qquad \forall (s_i,t_j) \in Q^+\end{aligned}$$ then $Z(s_i, t_j) \geq 0, \quad \forall (s_i,t_j) \in \bar Q^+$.

Let us examine the sign patterns of the second order operator ${\cal{L} }^{N,M}$ (defined in (\[patch-diff-op\])).  From assumption (\[M-and-N\]), we have that $$1-\frac{\bar k_j  s_i }{h} \geq 1- \frac{KL}{h} \geq 1- \frac{K}{2N} \geq 0.$$ Using this, we see that for the internal mesh points $$-{\cal{L} }^{N,M} Y(s_i,t_j) = \sum _{n=j-1}^{j+1} \sum _{k=i-1}^{i+1} a_{kn} Y(s_k,t_n);$$ where the sign of most of the coefficients $a_{kn}$ are easily identified to be $$\begin{aligned}
a_{i+1,j-1} =a_{i-1,j+1} =0;\quad
a_{i-1,j-1} = - \frac{s_i}{k_jh} \leq 0 &;&
a_{i+1,j+1} = -\frac{s_i}{k_{j+1}h} \leq 0; \\
a_{i,j} =\frac{2}{h^2} +\frac{2}{k_jk_{j+1}} \bigl(1+s_i^2 -  
 \frac{\bar k_j \vert s_i \vert)}{h} \bigr) &&>0; \\
a_{i,j-1} = - \frac{(1+s_i^2)}{\bar k_j k_j}+  \frac{s_i }{k_jh}=
 -\frac{1}{\bar k_jk_j} \bigl(1-\frac{\bar k_j  s_i }{h} +s_i^2\bigr) < 0 &;& a_{i,j+1} = -\frac{(1+s_i^2)}{\bar k_jk_{j+1}} +  \frac{ s_i }{hk_{j+1}} < 0.\end{aligned}$$ Finally, we look at the last two terms, $$a_{i-1,j} =a_{i+1,j} =-\frac{1}{h^2} + \frac{s_i}{hk_j}, \quad a_{i+1,j} =-\frac{1}{h^2} + \frac{s_i}{hk_{j+1}}.$$ Observe that $$\max \{a_{i-1,j}, a_{i+1,j} \} \leq  - \frac{1}{h^2} \bigl(1-\frac{Lh}{k}\bigr)
\leq- \frac{1}{h^2} \bigl(1-\frac{ML^2}{ 4N\tau _3}\bigr) \leq 0,$$ if we choose $L$ such that (\[def-L\]) is satisfied.  This sign pattern on the matrix elements insures that the system matrix associated with the finite difference scheme is an M-matrix [@fhmos pg.19], which suffices to establish the result. 

In the final phase, we solve the following discrete problem over the rectangle $$\label{downwind}
 S^*:= (L_*,R) \times (-R, R) \subset S,\qquad L_* :=2 \sqrt {NM^{-1} \tau _3};$$ using the mesh $ \Omega ^N_S$, which was defined in (\[Shish-par-mesh\]). 

[**Stage 4**]{}: [*Numerical method for $U_D$ defined over the downwind region $S^*$*]{}

Find $U_D$ such that

\[final-scheme\] $$\begin{aligned}
L^{N,M}U_D(x_i,y_j)     = 0 ,  \ (x_i,y_j) \in  \Omega ^N_S \setminus \bar C_1; 
\\
 U_D (L_*,y_j) =\left \{ \begin{array}{ll} \bar U_C (L_*,y_j), \  &y_j \in (-1-3\delta, 1+3\delta); \\
 \bar U_1 (L_*,y_j),\ &y_j \in [-R,-1-3\delta] \cup  [1+3\delta, R],\ \delta  \leq M/(12N);
\end{array}\right. \\ 
D^-_x  U_D (R,y_j)=0,  \ -R < y_j < R,\quad  U_D(x_i,-R)= U_D(x_i,R)=0, \ x_i \in [L_*,R] .  \end{aligned}$$

Then our corrected numerical approximation is given by $$\label{corrected-approx}
\bar U_2^N (x,y) := \left \{
\begin{array}{lll}
\bar U_1(x,y), & \hbox{for} \quad (x,y) \in \bar \Omega  \setminus (\{ x  \geq L^* \} \cup Q^+\cup Q^-) \\
\bar U_C(x,y) , & \hbox{for} \quad (x,y) \in  Q^+\cup Q^-, \\
\bar U_D(x,y) , & \hbox{for} \quad (x,y) \in (\{ x  \geq L^* \} \cap \Omega )\setminus (Q^+\cup Q^-)
\end{array}\right.  .$$ In the next section, we present some numerical results to illustrate the convergence properties of this corrected approximation, which is defined across three different coordinate systems.  A schematic image of the composite mesh $\tilde \Omega _A^N \cup \hat \Omega _P^N\cup \Omega _S^N$ is presented in Figure \[Fig.5\]. 

Numerical results
=================

In the previous section, we have seen that the initial approximations $\bar U_1^N$ displayed a lack of convergence, due to the presence of large errors in the neighbourhood of the characteristic points.  The corrected approximations $\bar U_2^N$ incorporate a parabolic patch near these points.  In the numerical experiments in this section, we have taken $M=N,\ \delta =0.05$ and, for ease of generating Tables, we have simply taken $L=  2 \sqrt{\min \{ \delta ,  \ve^{2/3}  \ln 2048 \}} $ in the patched region.  When we include the patch, we observe convergence in Table \[Tab3\] of the corrected approximations over this parabolic patch $\Omega _P^N$ over an extensive range of $\ve $ and $N$. In Table \[Tab4\], the global orders of convergence over the entire domain for the corrected approximation $\bar U_2^N$ are given.  These orders indicate that the corrected approximations are converging for all values of $\ve \in [2^{-20},1]$. In the final Table \[Tab5\], the approximate global errors over the entire domain are displayed for all $\ve \in [2^{-20},1]$. We observe that as $\ve \rightarrow 0$ the global errors continue to grow for each fixed $N$. Hence, the method appears not to be parameter-uniform.  Nevertheless, for any fixed value of $\ve$ we do observe convergence as $N $ increases.  In particular, we see in Figures \[Fig.8\] and \[Fig.9\] that for the corrected approximations $\bar U_2^N$, the approximate global errors $ \bar U_2 -\bar U_2^{2048}$ essentially halve as the number of mesh points are doubled.  This in sharp contrast to the approximate global errors $\bar U_1 -U_1^{2048}$ displayed in Figures \[Spike-N128\], \[Spike-N256\]. 

![ Approximate global error $ \bar U_2^{128}-\bar U_2^{2048}$ in corrected approximation for $\ve =2^{-10}$ []{data-label="Fig.8"}](Hemker_fig7)

![ Approximate global error $ \bar U_2^{256}-\bar U_2^{2048}$ in corrected approximation for $\ve =2^{-10}$ []{data-label="Fig.9"}](Hemker_fig8)

           $\varepsilon| N$     $N=8$        16        32        64      128      256      512
  ------------------------- --------- --------- --------- --------- -------- -------- --------
                        $1$    1.4739    0.3521    0.6916    1.1128   0.7936   0.9664   0.9716
                   $2^{-2}$    1.7825    0.2514    0.4127    1.1106   0.6181   0.9426   0.9508
                   $2^{-4}$    1.4364    0.6690    0.4384    0.8073   0.6223   0.8619   0.4013
                   $2^{-6}$    0.2310    0.7570    1.5302   -0.6245   2.0571   0.5677   0.2180
                   $2^{-8}$   -0.2532    0.3617    0.8478   -0.0286   1.1685   1.2720   0.1293
                  $2^{-10}$   -0.0002    0.0086    0.4198    1.0663   1.6521   1.2344   0.8949
                  $2^{-12}$   -0.0208   -0.0051    0.1917    0.8547   1.4528   1.2501   1.0148
                  $2^{-14}$    0.1888   -0.0509    0.1504    0.5917   1.2201   1.3731   1.0606
                  $2^{-16}$    0.3670   -0.0685    0.0477    0.3999   1.0007   1.6196   1.1642
                  $2^{-18}$    0.3904   -0.0179   -0.0203    0.2101   0.7656   1.4258   1.3736
                  $2^{-20}$    0.7398   -0.2868   -0.0467    0.0795   0.5156   1.2055   1.6702
    $\bar p^N_{R^{20}_\ve}$    0.3593    0.0521   -0.0351    0.1095   0.5156   1.2055   1.6702

  :  Computed double-mesh global orders of convergence $\bar p^{N}_{\epsilon}$ for the corrected approximations $\bar U_C$ (\[discrete-problem-transformed\]) measured over the patched region $Q$ (\[defQ\]), where $L= 2 \sqrt{\min \{ \delta ,  \ve^{2/3}  \ln 2048 \}} $ and $\delta =0.05$.[]{data-label="Tab3"}

           $\varepsilon| N$     $N=8$        16        32        64      128      256       512
  ------------------------- --------- --------- --------- --------- -------- -------- ---------
                        $1$    2.0261    0.6227    0.2352    0.7464   0.8763   0.9448    0.9799
                   $2^{-2}$    2.0506    1.0501    0.3720    0.3373   0.8437   0.9907    0.4472
                   $2^{-4}$    1.9773    0.1282    1.9129    0.8324   0.3435   0.3734   -0.3310
                   $2^{-6}$    0.2310    0.7570    1.5302   -0.6245   2.5573   0.2520    0.0335
                   $2^{-8}$   -0.2532    0.3617    0.8478   -0.0286   1.1685   1.2720    0.1293
                  $2^{-10}$   -0.0002    0.0086    0.4198    1.0663   1.6521   1.2344    0.8949
                  $2^{-12}$   -0.0208   -0.0051    0.1917    0.8547   1.4528   1.2501    1.0148
                  $2^{-14}$    0.1888   -0.0509    0.1504    0.5917   1.2201   1.3731    1.0606
                  $2^{-16}$    0.3670   -0.0685    0.0477    0.3999   1.0007   1.6196    1.1642
                  $2^{-18}$    0.3904   -0.0179   -0.0203    0.2101   0.7656   1.4258    1.3736
                  $2^{-20}$    0.7398   -0.2868   -0.0467    0.0795   0.5156   1.2055    1.6702
    $\bar p^N_{R^{20}_\ve}$    0.3593    0.0521   -0.0351    0.1095   0.5156   1.2055    1.6702

  : Computed double-mesh orders of global convergence $\bar p^{N}_{\epsilon}$ for the corrected approximations $\bar U_2$ (\[corrected-approx\]) measured over the entire domain $\Omega$[]{data-label="Tab4"}

    $\varepsilon| N$    $N=8$       16       32       64      128      256      512
  ------------------ -------- -------- -------- -------- -------- -------- --------
                 $1$   0.1689   0.0828   0.0686   0.0435   0.0245   0.0127   0.0057
            $2^{-2}$   0.2735   0.1434   0.0765   0.0530   0.0311   0.0197   0.0100
            $2^{-4}$   0.4079   0.2272   0.1415   0.0887   0.0515   0.0305   0.0169
            $2^{-6}$   0.5720   0.4206   0.2369   0.1186   0.1109   0.0330   0.0188
            $2^{-8}$   0.6790   0.5657   0.4065   0.2034   0.1852   0.0650   0.0123
           $2^{-10}$   0.7331   0.6325   0.4882   0.3030   0.1344   0.0518   0.0198
           $2^{-12}$   0.7769   0.6916   0.5624   0.3859   0.1971   0.0698   0.0326
           $2^{-14}$   0.8140   0.7414   0.6287   0.4662   0.2718   0.1100   0.0446
           $2^{-16}$   0.8458   0.7847   0.6871   0.5419   0.3462   0.1615   0.0547
           $2^{-18}$   0.8718   0.8213   0.7393   0.6122   0.4329   0.2280   0.0799
           $2^{-20}$   0.8950   0.8519   0.7830   0.6720   0.5095   0.3025   0.1213

  : Approximate global errors $\Vert \bar U^N_2 - \bar U^{2048}_2 \Vert $ over the domain $\Omega$[]{data-label="Tab5"}

The computed orders of convergence in Table \[Tab4\] suggest that this multi-stage numerical method is producing a converging sequence of numerical approximations to the analytical solution of the Hemker problem (\[cont-prob\]) across the entire domain for singular perturbation values $\ve \in [2^{-20},1]$.

Conclusions
===========

Based on parameter explicit pointwise bounds on how the continuous solutions decays away from the circle, a numerical method was constructed for the Hemker problem.  There are no spurious oscillations present in the numerical solutions, as we use simple upwinding in all co-ordinate directions used.  Several layer adapted Shishkin meshes are utilized and these grids are aligned both to the geometry of the domain and to the dominant direction of decay within the boundary/interior layer functions.  Numerical experiments indicate that the method is producing accurate approximations over an extensive range of the singular perturbation parameter. 

[9]{}

M. Augustin, A. Caiazzo, A. Fiebach, J. Fuhrmann, V. John, A. Linke and R. Umla, An assessment of discretizations for convection-dominated convection-diffusion equation, [*Comput.  Methods Appl.  Mech.  Engrg.*]{}, [** 200**]{} (47-48), 2011, 3395–3409. 

G. Barrenechea, V. John, P. Knobloch and R. Rankin, [A unified analysis of algebraic flux correction schemes for convection-diffusion equations]{}, [*SeMA Journal.  Boletin de la Sociedad Espanñola de Matemática Aplicada*]{}, [**75**]{} (4), [2018]{}, [655–685]{}. 

P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan and G. I. Shishkin, [*Robust Computational Techniques for Boundary Layers*]{}, Chapman and Hall/CRC Press, Boca Raton, (2000). 

R. K. Dunne, E. O’ Riordan and G. I. Shishkin, Fitted mesh numerical methods for singularly perturbed elliptic problems with mixed derivatives, [*IMA J. Num.  Anal.*]{}, [** 29**]{}, 2009, 712–730. 

W. Eckhaus, Boundary layers in linear elliptic singular perturbation problems, [*SIAM Rev.*]{}, [**14**]{}, 1972, [225–270]{}. 

B. García-Archilla, [Shishkin mesh simulation: a new stabilization technique for convection-diffusion problems]{}, [*Comput.  Methods Appl.  Mech.  Engrg.*]{}, [** 256**]{}, [2013]{}, [1–16]{}. 

H. Han, Z. Huang and R. B. Kellogg, A tailored finite point method for a singular perturbation problem on an unbounded domain. [* J. Sci.  Comput.*]{}, [** 36**]{} (2), 2008, 243–261. 

E. D. Havik, P. W. Hemker and W. Hoffmann, Application of the over-set grid technique to a model singular perturbation problem, [*Computing*]{}, [**65**]{} (4), [2000]{}, [339–356]{}. 

A. F. Hegarty and E. O’Riordan, A parameter-uniform numerical method for a singularly perturbed convection-diffusion problem posed on an annulus, [*Comput.  Math.  Appl.*]{}, [** 78**]{} (10), 2019, 3329–3344, . 

A. F. Hegarty and E. O’Riordan, A numerical method for the Hemker Problem, [*Boundary and Interior Layers - Computational and Asymptotic Methods*]{}, BAIL 2018, Glasgow, G. R. Barrenechea and J. Mackenzie (Eds.),  Lecture Notes in Computational Science and Engineering, vol.  135, Springer, 2020 (to appear). 

P. W. Hemker, A singularly perturbed model problem for numerical computation, [*J. Comp.  Appl.  Math*]{}, [**76**]{}, 1996, 277–285. 

A. M. Il’in, [*Matching of asymptotic expansions of solutions of boundary value problems*]{}, Mathematical Monographs, [**102**]{}, American Mathematical Society, (1992). 

V. John and L.  Schumacher, [A study of isogeometric analysis for scalar convection-diffusion equations]{}, [*Appl.  Math.  Lett.*]{}, [** 27**]{}, [2014]{}, [43–48]{}. 

P. A.  Lagerstrom, [*Matched asymptotic expansions: Ideas and techniques*]{}, Applied Mathematical Sciences, [**76**]{}, [Springer-Verlag, New York]{}, (1988). 

J. J. H. Miller, E. O’Riordan and G. I. Shishkin, [*Fitted numerical methods for singular perturbation problems*]{}, World-Scientific, Revised edition, (2012). 

E. O’Riordan and G. I. Shishkin, A technique to prove parameter–uniform convergence for a singularly perturbed convection–diffusion equation, [*J. Comp.  Appl.  Math*]{}, [**206**]{}, 2007, 136–145.  E. O’Riordan and G. I. Shishkin, Parameter uniform numerical methods for singularly perturbed elliptic problems with parabolic boundary layers, [*Appl.  Numer.  Math.*]{}, [** 58**]{}, 2008, 1761–1772. 

M. H. Protter and H. F. Weinberger, [*Maximum Principles in Differential Equations*]{}.  Springer–Verlag, New York, (1984). 

G. I. Shishkin, [*Discrete approximation of singularly perturbed elliptic and parabolic equations*]{}, Russian Academy of Sciences, Ural section, Ekaterinburg, (1992). ( in Russian)

R. T. Waechter, Steady longitudinal motion of an insulating cylinder in a conducting fluid, [*Proc.  Camb.  Philos.  Soc.*]{}, [** 64**]{}, 1968, 1165–1201. 

[**Appendix 1: Bounds on the continuous solution**]{}

In this appendix, a barrier function $B(x,y)$ is constructed, which involves the construction of four separate sub-components $B_i(x,y) \geq 0 , i=1,2,3,4$ with the properties that $$LB_i (x,y) \geq 0, \ (x,y) \in \Omega; \quad  B_i(x,y) \geq 1, \ \hbox{if} \ x^2+y^2 =1  \quad \hbox{and} \quad  (B_i)_x(R,y) \geq 0 .$$ Each sub-component incorporates a smooth cut-off function $\psi _i$, which limits the support of each $B_i(x,y) $ to a proper subset $\Omega ^*_i \subset \bar \Omega $ of the domain.  In addition, for each sub-component we will define a proper subset $\Omega _i \subset \Omega ^*_i$ over which the cut-off function $\psi _i$ is a non-zero constant.  In fact, we will construct these functions so that $$\psi _i \equiv 1 \quad \hbox{on} \quad \Omega _i\qquad \hbox{and} \quad \psi _i \equiv 0 \quad \hbox{on} \quad \Omega \setminus \Omega ^*_i.$$ Moreover, $\bar \Omega = \cup _{i=1}^4 \Omega _i$ and $B(x,y) := \sum _{i=1}^4 B_i(x,y)$.

1.   As in [@annulus appendix] we can construct a smooth cut-off function $ \psi _1(\theta ): [0,\pi] \rightarrow [0,1] $, with the following properties:

    \[cut-off\] $$\begin{aligned}
    - \ve ^2 \psi _1'' - \ve   \psi _1' + C_* \ve ^{2/3} \psi _1 &=&0, \ \theta \in [0,\pi]; \quad \psi _1 '(\theta ) \geq 0,  \forall \ \theta \in (0,\pi ); \\ 
    \psi  _1 \equiv 0, \ \theta \in    [0, \frac{\pi}{2} +\theta _* ], && \psi _1 \equiv 1,  \theta \in    [\frac{\pi }{2} +2\theta _*, \pi ], \\
    \ve ^{2/3}C_* \leq (1-\kappa)\kappa\cos ^2(\frac{\pi}{2} +\theta _*) &,&  \sin (2\theta _*) = \ve ^{1/3},\end{aligned}$$

    where $0<\kappa<1$ is arbitrary.  Observe that this cut-off function is only non-constant over an interval of magnitude $O(\ve ^{1/3})$. For all $\theta \in [0,\pi]$, consider the following barrier function $$\label{barrier}
    \tilde B^-(r,\theta ) := \psi _1(\theta) E_-(r,\theta),  \quad E_-(r,\theta) :=e^{\frac{ \kappa \cos (\theta ) (r-1)}{\ve}},\ 0< \kappa < 1.$$ Note the following expressions for the partial derivatives of this barrier function: $$\begin{aligned}
    \tilde B^-_\theta &=& -\frac{ \kappa \sin (\theta ) }{\ve} (r-1) B^- + \psi _1'(\theta) E_-; \\
    \tilde B^-_{\theta \theta} 
    &=& -\bigl(  \frac{\kappa \cos (\theta) }{\ve} (r-1) - \frac{ \kappa ^2 \sin ^2 (\theta) }{\ve ^2} (r-1)^2 \bigr) B^-  + (\psi _1''(\theta) -2\psi _1'(\theta ) \frac{ \kappa \sin(\theta) }{\ve} (r-1)) 
    E_-; \\
    \tilde B^-_r &=& \frac{\kappa \cos (\theta) }{\ve}  B^-,   \qquad
    \tilde B^-_{rr} = \frac{\kappa ^2 \cos ^2 (\theta) }{\ve ^2}  B^-.\end{aligned}$$ Combining these expressions, we can deduce (see [@annulus appendix] for details) that $$\begin{aligned}
    \tilde L \tilde B^- %&\geq& \frac{\kappa \tilde B^-}{\ve}\bigl( \cos ^2 (\theta) (1-\kappa ) -\ve \vert \cos \theta\vert \bigr) + \bigl(-\ve \psi _*'' -  \psi _* ' \bigr)  E_- \\
    &\geq&
    \frac{-\ve ^2\psi _1'' -  \ve \psi _1 '+  (1-\kappa)\kappa\cos ^2\theta    \psi _1}{\ve }   E_-, \qquad \forall \theta \in [0, \pi].\end{aligned}$$ Hence, by our design of the cut-off function $\psi _1 (\theta)$ $$\tilde L \tilde B^-(r,\theta) \geq 0 \quad \hbox{and} \quad \tilde B^-(1,\theta ) =1 \qquad \hbox{if} \quad  \frac{\pi }{2} +2\theta _*  \leq  \theta \leq \pi.$$ This barrier function can be extended to the region $\theta \in [\pi,2\pi]$ by defining $$\begin{aligned}
    \tilde B^-(r,\theta ) :=  \psi _1(2\pi-\theta) e^{\frac{ \kappa \cos (\theta ) (r-1)}{\ve}}, \quad \theta \in [\pi, 2\pi].\end{aligned}$$ In the rectangular coordinates, we now define the first sub-component to be: $$\begin{aligned}
    B_1(x,y) &:=& \tilde B^-(r,\theta ), \quad \forall (x,y) \in \Omega _1^*, \quad \hbox{where} \\
    \bar \Omega _1^*  &:=&  \{ (x,y) \vert x \leq -\arcsin(\theta _*), -R \leq y \leq R \} \cap \Omega  \quad \hbox{and} \\
    \bar \Omega _1  &:=&  \{ (x,y) \vert x \leq -\ve ^{1/3}, -R \leq y \leq R \} \cap \Omega .\end{aligned}$$ Hence, $LB_1 = \tilde L \tilde  B^- \geq 0$ and $ \ B_1(x,y) \geq 0$, for all $ (x,y) \in \bar \Omega _1^*.$ Moreover, $$B_1(x,y) \geq 1, (x,y) \in \bar \Omega_1 \cap \partial \Omega _I \quad \hbox{and} \quad \frac{\partial B_1}{\partial x}(x,y) \geq 0, (x,y) \in \bar \Omega \cap \partial \Omega _N.$$

2.   Consider the following barrier function $$B(x,y) := \bigl(1+\frac{\alpha x}{\ve ^{1/3}} \bigr)(1-\alpha C)^{-1}e^{-\frac{\kappa (0.5x^2+y-1)}{\ve ^{2/3}}},\quad -C\ve ^{1/3} < x < C\ve ^{1/3}, \quad \alpha C < 1;$$ where the possible ranges for the parameters $\alpha, \kappa, C$ will be specified below.  Note the following expressions for the partial derivatives of this barrier function: $$\begin{aligned}
     B_x &=& \frac{(1-\alpha C)^{-1}}{\ve ^{1/3}}\bigl( \alpha - \frac{\kappa x}{\ve ^{{1/3}}} (1+\frac{\alpha x}{\ve ^{1/3}}) \bigr)e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}}; \\
    \ve  B_{yy} &=& \frac{\kappa ^2(1-\alpha C)^{-1}}{\ve ^{1/3}}\bigl(1+\frac{\alpha x}{\ve ^{1/3}}\bigr)e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}}; \\
    \tilde B_{xx} 
    &=&  \frac{1}{\ve } O(\ve ^{1/3})e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}},\quad -C\ve ^{1/3} < x < C\ve ^{1/3}.\end{aligned}$$ Let us introduce $\xi := x \ve ^{-1/3}$, then for $\xi \in C (-1,1)$, $0< \alpha C < 1$, $$\begin{aligned}
    -\ve  B_{yy} + B_x &=&\frac{(1-\alpha C)^{-1}}{\ve ^{1/3}}\bigl(\alpha - \kappa \xi (1+\alpha \xi) - \kappa ^2 (1+\alpha \xi )\bigr)e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}} \\
    &\geq &\frac{1}{\ve ^{1/3}}\bigl(\frac{\alpha}{(1+\alpha C)}  - C\kappa - \kappa ^2 \bigr)e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}} \\
    %&\geq &\frac{(1+\alpha C)^2}{\ve ^{1/3}}\bigl(\frac{4\alpha+C^2(1+\alpha C)}{4(1+\alpha C)}  - (\kappa +\frac{C}{2})^2\bigr)e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}} \end{aligned}$$ In order that $ LB(x,y) \geq 0, x \in C\ve ^{1/3} (-1,1)$ (and recalling that $\alpha C < 1$) we can impose the following constraints: $$\begin{aligned}
    \label{constraints}
    %\alpha C <1, && \hbox{so that} \quad \tilde B(-C\ve ^{1/3},y) > 0, \qquad \hbox{and   }\\
    %\hbox{in order that}\quad  \tilde L \tilde B \geq 0, &&  \hbox{we impose} \quad
     \kappa  \leq  \frac{C}{3}, \quad \alpha >\frac{8}{9} C^2, \quad  C^3 < \frac{9}{8}.\end{aligned}$$ Note also that over the set of points $$S_1:= \{(x,y) \vert x^2+y^2 \geq 1\} \cap \{(x,y) \vert-C\ve ^{1/3} \leq x \leq C\ve ^{1/3}\} \cap\{(x,y) \vert 0.5x^2+y \leq 1 \}$$ we have that $$\begin{aligned}
     B(x,y) \geq  (1+\frac{\alpha x}{\ve ^{1/3}})(1-\alpha C)^{-1} \geq 1 >0, \quad  (x,y) \in S_1.\end{aligned}$$ In addition, $$\hbox{if} \quad
    x=\ve ^{1/3}, \hbox{then on the boundary}\quad  x^2+y^2=1, \quad y= 1- \frac{\ve ^{2/3}}{\sqrt{1-\ve^{2/3}}+1}.$$ Hence, $B(x,y)$ is uniformly bounded on the domain $\Omega  \cap \{ \vert x \vert \leq C \ve ^{1/3} \} $. 

    Now for some positive parameter $\mu$, we wish to construct a smooth cut-off function $$\psi _2 (x): C(1+\mu) \ve ^{1/3}[-1, 1] \rightarrow [0,1],$$ with certain properties.  First, for any function $ \psi  (x)$, $$\begin{aligned}
     L ( \psi (x)  B) &=&  \psi  L (  B) +  \psi '  B- 2 \ve \psi 'B_x - \ve \psi '' B \\
    &=&  \psi  L (  B) + 
    \bigl( (- \ve \psi ''+  \psi ')(1+ \alpha \xi )- 2 \ve ^{1/3} \psi ' \bigl( \ve ^{1/3}\alpha - \kappa \xi (1+\alpha \xi ) \bigr) e^{-\frac{ \kappa (0.5x^2+y-1)}{\ve ^{{2/3}}}}.\end{aligned}$$ We construct the cut-off function $ \psi _2 (x)$ so that, for some $\mu >0$,

    \[cut-off2\] $$\begin{aligned}
    - \ve   \psi _2'' +   0.5  \psi _2'   =0, \  \psi _2 '(x) \geq 0,  \forall \ x \in C(1+\mu) \ve ^{1/3}(-1, 0); \\ 
    - \ve   \psi _2'' +   1.5  \psi _2'   =0,  \psi _2 '(x) \leq 0,  \forall \ x \in C(1+\mu) \ve ^{1/3}(0,1); \\ 
     \psi _2  \equiv 0, \ x \in    C\ve ^{1/3}([-(1+\mu),-(1+0.5\mu)] \cup [1+0.5\mu, 1+\mu]), \\
     \psi _2 \equiv 1,  x \in   C\ve ^{1/3}[-1, 1].  \end{aligned}$$

    Then, if we have the strict inequality $\alpha C <1$ and for $\ve$ sufficiently small, we have that $$\begin{aligned}
     L (\psi _2 B) \geq 0, \quad \quad \hbox{if} \quad x \in C\ve ^{1/3}(1+\mu) [-1,1] \quad \hbox{ and} \quad x^2+y^2 \geq 1,\ y \geq 0; 
    \\ 
     \psi _2 B \geq 1, \quad \quad \hbox{if} \quad x \in C\ve ^{1/3}[-1,1] \quad \hbox{ and} \quad x^2+y^2 = 1, y \geq 0.  \end{aligned}$$ To match this barrier function with the first barrier function $B_1(x,y)$ and to also satisfy the constraints (\[constraints\]), we shall take $$C =1, \ \alpha =0.9  \quad \hbox{and} \quad  \kappa =\frac{1}{3} .$$ With these choices for the parameters, we define our second sub-component to be $$\begin{aligned}
    B_2(x,y) &:=&  \psi _2 (x) \bigl(10+\frac{9 x}{\ve ^{1/3}} \bigr)e^{-\frac{ (0.5x^2+\vert y \vert -1)}{3\ve ^{2/3}}}, \quad \forall (x,y) \in \Omega _2^*, \quad \hbox{where} \\
    \bar \Omega _2^*    &:=&  \{ (x,y) \vert -\ve ^{1/3} (1+\mu)  \leq  x \leq (1+\mu) \ve ^{1/3}, 0<\mu <1, \ -R \leq y \leq R \} \cap \Omega    \\
    \bar \Omega _2  &:=&  \{ (x,y) \vert -\ve ^{1/3} \leq x \leq \ve ^{1/3}, -R \leq y \leq R \} \cap \Omega .\end{aligned}$$

3.   Consider the following function $$B^+(x,y) := e^{C_1x} e^{-\frac{\mu (y-1)}{\sqrt{\ve}} },\quad x\geq 0, \ -C_3\sqrt{\ve} \leq y-1 ;\quad  C_1 \geq 2+\mu ^2.$$ Note the following expressions for the partial derivatives of this function: $$B^+_x = C_1B^+\quad \hbox{and} \quad  -\ve  B^+_{yy} = -\mu ^2  B^+$$ and so, for $\ve$ sufficiently small, $$L  B^+ 
    =  \bigl(C_1- \ve C_1^2 - \mu ^2 \bigr)B^+ \geq 0, \quad \forall  y.$$ As before we will confine the influence of this barrier function via a cut-off function $$\psi _3 (y): [-R, R] \rightarrow [0,1] .$$ Note that $$\begin{aligned}
     L ( \psi (y)  B^+) = \psi  L (  B^+) - 2 \ve   \psi  ' B^+_y - \ve  \psi  '' B^+ 
    = \psi  L (  B ^+) + 
    (- \ve  \psi ''+ 2\sqrt{\ve}  \psi  ')  B ^+.\end{aligned}$$ Then, we construct a smooth cut-off function $\psi _3  (y)$ so that

    \[cut-off3\] $$\begin{aligned}
    &&  \psi _3  '(y) \geq 0,  \forall \ y \in (-R, R); \\
    && \psi _3  \equiv 0, \ y \in  (-R,1- 2C_3\sqrt{\ve}], \quad \psi _3 \equiv 1,  y \in  [1- C_3\sqrt{\ve}, R];\\
    &&- \ve   \psi _3'' +2\sqrt{\ve} \psi _3 '  +  \psi _3= 0, \forall \ y \in (1-2C_3\sqrt{\ve}, 1- C_3\sqrt{\ve}).\end{aligned}$$

    Observe that $ \psi _3 (y) B^+(x,y) \geq 1,$ if $ 1- C_3\sqrt{\ve} \leq y \leq 1$ and $\psi _3 (y)B _x(R,y) \geq 0$. We now define the third subcomponent for the barrier function $$\begin{aligned}
    B_3(x,y) &:=& e^{C_1x} \Bigl(   \psi _3 (y) e^{-\frac{(y-1)}{\sqrt{\ve}} } +(1-  \psi _3 (y+2))e^{\frac{(y+1)}{\sqrt{\ve}} }\Bigr), \quad \forall (x,y) \in \Omega _3^*, \quad \hbox{where} \\
    \bar \Omega _3^*    &:=&  \{ (x,y) \vert x \geq \ve ^{1/3},  y \in [- R, -1+2C_3\sqrt{\ve}] \cup  [1-2C_3\sqrt{\ve},R] \} \cap \Omega   \quad \hbox{and} \\
    \bar \Omega _3  &:=&  \{ (x,y) \vert  x \geq \ve ^{1/3},  y \in [- R, -1+C_3\sqrt{\ve}] \cup  [1-C_3\sqrt{\ve},R] \} \cap \Omega ;\end{aligned}$$ where the constant $C_3$ is an arbitrary positive constant.  Note that $$B_3(x,y) \geq 1 , \ \quad \hbox{if} \quad x^2+y^2 =1.$$

4.   Finally the fourth sub-component of the barrier function is defined by $$\begin{aligned}
    B_4(x,y) &:=&  e^{2x} \psi _4(y), \quad \forall (x,y) \in \Omega _4^*, \quad \hbox{where} \\
    \bar \Omega _4^*    &:=&  \{ (x,y) \vert x \geq \ve ^{1/3}, -1+0.5C_3\sqrt{\ve} \leq y \leq 1-0.5C_3\sqrt{\ve} \} \cap \Omega   \quad \hbox{and} \\
    \bar \Omega _4 &:=&  \{ (x,y) \vert  x \geq \ve ^{1/3},   -1+C_3\sqrt{\ve} \leq y \leq 1-C_3\sqrt{\ve} \} \cap \Omega .\end{aligned}$$ Note that, for $\ve$ sufficiently small, $$LB_4= e^{2x}(-\ve \psi ''_4 +2(1-2\ve)\psi _4) \geq 0,$$ when the cut-off function $\psi _4   (y): [-R, R] \rightarrow [0,1] $ is an even function constructed so that

    \[cut-off4\] $$\begin{aligned}
    &&  \psi _4  '(y) \leq 0,  \ \psi _4  (-y) = \psi _4 (y) \quad \forall \ y \in (0, R);\\
    && \psi _4  \equiv 0, \ y \in  [1-0.5C_3\sqrt{\ve}, R]; \quad 
     \psi _4 \equiv 1,  y \in  [0, 1- C_3\sqrt{\ve}];\\
    &&- \ve   \psi _4''   +  \psi _4= 0, \forall \ y \in (0, 1-0.5C_3\sqrt{\ve}).\end{aligned}$$

We can now apply Theorem \[compar2\], across the entire domain $ \bar \Omega = \cup _{i=1}^4  \bar \Omega _i$, to deduce that $$0 \leq  u(x,y) \leq \sum _{i=1}^4  B_i(x,y), \quad (x,y) \in \bar \Omega.$$

[^1]: The first author acknowledges the support of MACSI, the Mathematics Applications Consortium for Science and Industry (www.macsi.ul.ie), funded by the Science Foundation Ireland Investigator Award 12/IA/1683. 

[^2]: Department of Mathematics and Statistics, University of Limerick, Ireland. 

[^3]: School of Mathematical Sciences, Dublin City University, Dublin 9, Ireland. 

[^4]: Over any computational cell $\Omega _{i,j}:= (x_{i-1},x_i)\times (y_{j-1},y_j)$, where $h_i:=x_i-x_{i-1}, k_j:=y_j-y_{j-1}$, we denote the bilinear interpolant of any function $g(x,y)$ by $\bar g$. For any smooth function $v$, where $\Vert v_{xx} \Vert + \Vert v_{yy} \Vert \leq C$, then $
    \Vert v -\bar v \Vert _{\Omega _{i,j}, \infty}\leq C h_i^2 +C k_j^2. 
    $ For a layer function of the form $w(x,y) = e^{-\alpha x/\ve ^p}$, $$\begin{aligned}
    \Vert w -\bar w \Vert _{\Omega _{i,j}, \infty}
    \leq C \Vert w \Vert _{\Omega _{i,j}, \infty}\leq CN^{-1}, \quad \hbox{if} \quad x_i \geq \frac{\ve ^p \ln N}{\alpha}, \quad \hbox{and}  \\
    \Vert w -\bar w \Vert _{\Omega _{i,j}, \infty}\leq C h_i\Vert w _x\Vert _{\Omega _{i,j}, \infty}
    \leq CN^{-1}\ln N, \quad \hbox{if} \quad h_i \leq C\ve ^p N^{-1} \ln N.\end{aligned}$$

[^5]: The finite difference operators $D^+_r,D^-_r, D^\pm_r, \delta^2_r$ are, respectively, defined by $$\begin{aligned}
    D^+_rZ(r_i,\theta _j) :=\frac{Z (r_{i+1},\theta _j)-Z(r_i,\theta _j)}{r_{i+1}-r_i},\quad  D_r^-Z(r_i,\theta _j) :=\frac{Z(r_{i},\theta _j)-Z (r_{i-1},\theta _j)}{r_{i}-r_{i-1}}; \\
    2(bD_r^\pm)Z := (b-\vert b \vert) D_r^+Z + (b+\vert b \vert) D_r^-Z; \quad \delta^2_r Z(r_i,\theta _j) :=\frac{D^+_rZ(r_{i},\theta _j)-D_r^-Z(r_i,\theta _j)}{(r_{i+1}-r_{i-1})/2}.\end{aligned}$$
