Numerical Solution of a Singularly Perturbed Problem on a Circular Domain

. We consider a singularly perturbed elliptic problem, of convection-diﬀusion type, posed on a circular domain. Using polar coordinates, simple upwinding and a piecewise-uniform Shishkin mesh in the radial direction, we construct a numerical method that is monotone, pointwise accurate and parameter-uniform under certain compatibility constraints. Numerical results are presented to illustrate the performance of the numerical method when these constraints are not imposed on the data.


Introduction
There have been many recent publications on parameter-uniform numerical methods [9,2] for singularly perturbed linear problems of the convection-diffusion form where the domain is the unit square.In contrast, there have been few publications for the case when the domain is non-rectangular and the problem is of convection-diffusion type.In the case of elliptic problems posed on arbitrary convex domains, several theoretical difficulties arise in designing a parameter-uniform numerical method [1].An invertible transformation can sometimes be designed to map a non-rectangular domain to the unit square.However, in general, the differential operator in the transformed variables will contain a mixed second order derivative term.The construction of a parameteruniformly stable discretization of a mixed second order derivative on a layer-adapted highly anisotropic mesh remains an open question.In this paper, we examine a circular domain for which the standard transformation using polar coordinates can be utilized throughout the entire domain.For this particular geometry, there are no mixed derivative terms present in the transformed problem.
The issue of compatibility conditions at the characteristic points and the construction of an ε-uniform asymptotic expansions have been studied by Jung and Teman [6].In [3], we impose more stringent constraints on the problem data in a neighbourhood of the characteristic points, in order to exclude the potential presence of additional singularities near these two points.Under these data constraints, a parameter-uniform numerical method can be constructed [3], which captures the boundary layer at the outflow boundary.Here, we present numerical results for a problem which does not satisfy the minimal compatibility conditions.These experimental results suggest that some positive order of uniform convergence may be retained in the case of no compatibility.However, a theoretical justification for this conjecture remains an open question.
Notation: Our interest lies in designing parameter-uniform numerical methods and so throughout this paper, C denotes a generic constant that is independent of the singular perturbation parameter ε and of all discretization parameters.We will always use the pointwise maximum norm, which we denote throughout by • .

Continuous problem
Consider the singularly perturbed elliptic problem: Assume that the data are sufficiently smooth so that ũ ∈ C 3,α ( Ω).For this problem, a boundary layer will form in the vicinity of the outflow boundary } and there will be no layer present in the vicinity of the inflow boundary Polar coordinates are a natural co-ordinate system to employ for this problem.In these coordinates, the continuous problem (1) is transformed into the problem: Find The solution u can be decomposed into the sum of a regular component v and a singular component w such that This type of decomposition of the solution was first introduced by Shishkin [9] for a large class of singularly perturbed partial differential equations.The decomposition is related to an asymptotic expansion [8], but note that there is no explicit identification of a remainder term.Moreover, the components v and w are not explicitly identified.
From a numerical analysis perspective, the advantage of this type of decomposition is that parameter explicit pointwise bounds on the partial derivatives (up to and including third order) of these two components can be established.These bounds on the derivatives of the components are central to establishing informative pointwise bounds on the truncation error associated with any proposed numerical method.
The reduced problem is defined as: Find ṽ0 such that As identified in [6], singularities appear in the vicinity of the points (±1, 0) unless compatibility conditions of some level m are imposed on the data.The first order correction to the reduced solution is given by Requiring a certain regularity on ṽ1 places additional regularity requirements on ṽ0 .Explicit compatibility conditions on the derivatives of the regular components ṽ0 , ṽ1 are given in [5,Lemma 2.2] to ensure any desired level of regularity of these components.The outflow boundary conditions for the regular component are taken to be v 0 + εv 1 in order that the partial derivatives (up to second order) of the regular component are bounded independently of ε.For example, compatibility conditions of level m = 9 suffice for ṽ0 ∈ C 5 (Ω) and ṽ1 ∈ C 5 (Ω).However, to obtain pointwise bounds on the boundary layer component w, additional constraints were imposed in [3]: Assumption Assume that there exists a 0.5 < δ < 1 such that Note that this assumption supersedes the compatibility conditions (3) of any order.
Theorem 1. [3] The solution u of problem (1), ( 4) can be decomposed into the sum u = v + w, where the derivatives of the regular component v satisfy the bounds and the boundary layer component satisfies (for some positive γ)

Discrete problem
A popular assumption within the literature on numerical methods for singularly perturbed problems is to assume that ε ≤ CN −1 .
In this case a classical finite element method will generate large oscillatory solutions and thereby fail to capture accurately any layers present in the solution.This is one reason why this case is often viewed as the case of most interest.For certain classes of singularly perturbed problems, one can generate a uniformly valid asymptotic expansion for the continuous solution.In addition, there are problem classes for which a combination of a simple numerical method and the analytical expression for the leading order term for the layer function can be utilized to generate an approximation, to the continuous solution, of O(N −1 ).For example, in the case of problem (1), ( 4), we observe that Coupling this asymptotic expansion with the restriction on the parameters of ε ≤ CN −1 , means that we only require a sufficiently accurate approximation to the reduced solution ṽ0 (x, y).This can be easily generated by simply discretizing the first order problem defining the reduced problem and then interpolating these nodal values to produce a global approximation.However, these mixed numerical/asymptotic approaches only have validity under the assumption ε ≤ CN −1 .In [3], we have designed a parameter-uniform numerical method, which is valid for all values of 0 < ε ≤ 1.
We discretize problem (2), (4) using simple upwinding on a tensor product mesh, with M mesh elements uniformly distributed in the angular direction and N mesh elements in the radial direction distributed across a piecewise uniform Shishkin mesh [9].The mesh points (r i , θ j ) are defined by: The numerical method on this mesh is: For where 2(aD The finite difference operators D + , D − are the standard forward and backward first order difference operators and δ 2 denotes the standard discrete approximation of the second derivative.Once the nodal values have been determined, a global approximation Ū can be generated using the bilinear interpolant where φ i (r), ψ j (θ) are piecewise linear basis functions, defined by the nodal values of If u is the solution of the continuous problem (1), ( 4) and Ū is the bilinear interpolant of the discrete solution U , then

Numerical results
In this section, we examine the performance of the numerical method applied to two particular problems.However, these problems do not satisfy the constraint (4) imposed in [3].Observe that the choice of the transition parameter σ in (5c) depends on (4).For the problems examined in this section, we have simply replaced (5c) by Approximations to the uniform order of convergence are estimated using the double mesh method [2].The numerical method was applied for ε ∈ {2 −j } 20 0 and N ∈ {2 j } 10 3 .The maximum pointwise two-mesh differences D N ε and the parameter-uniform maximum pointwise two-mesh differences D N , are computed from Approximations p N ε to the local order of convergence and approximations p N to the parameter-uniform order of local convergence are subsequently computed from Example 1 (Compatible problem) Motivated by the test example from [4, equation (61)], we consider the following problem This problem is not covered by the theory presented in [3].However, it is highly compatible at the characteristic points.In Table 1, we present a corresponding table of computed orders of convergence and in Figure 1 we present a sample computed solution and a comparison between the computed solution and a fine mesh solution.We observe In Figure 2, we see some significant difference (between the computed solution and a fine mesh solution) in the vicinity of the two characteristic points.In Table 2, we present a corresponding table of computed orders of convergence.We again observe convergence in the Table.This suggests that the lack of compatibility in the reduced problem does not appear to have a significant detrimental effect on the performance of the method in the outflow region.These numerical results indicate that the question of whether the above numerical method is parameter-uniform or not, for incompatible problems like (8), warrants further investigation.
Fig 1. Plots of numerical solution and approximate error for Example 1 with ε = 2 −15 .