The study of singularly perturbed problems has flourished since the publication of Prandtl's seminal work in 1904 on “boundary layers” [1]. Many researchers have paid attention to the theoretical and computational aspects of those problems. The usual task has been to provide means of dealing with the challenges that come with the perturbation parameter and its impact on the solution behavior. While countless successes have been recorded in the case of linear singularly perturbed problems [see for example [2–7]], little attention has been paid to the non-linear case.

In this article, we study the two-dimensional singularly perturbed semilinear convection-diffusion problems

subject to boundary conditions

where ε is the perturbation parameter with 0 < ε ≪ 1. The semilinear source term f(x, y, u(x, y)) and the coefficient functions a₁(x, y), a₂(x, y) are assumed to be sufficiently smooth and satisfy

where α₁, α₂ and β are constants and ∂Ω is the boundary of Ω. Under these conditions (1.1)-(1.4) has a unique solution which displays boundary layers at x = 1 and y = 1 when ε approaches zero.

Problems such as (1.1)-(1.4) are encountered in diverse areas of applied mathematics and engineering such as aerodynamics, liquid crystal modeling, chemical reactor theory, magnetohydrodynamics, oceanography, fluid mechanics, heat conduction, quantum mechanics [see [8–13]]. The difficulty with such problems is that researchers have to deal with both the perturbation parameter and the complexity due to the semilinearity, besides the higher-dimensional aspect. Perhaps, that is the reason why only few people have shown some interest in them.

Sirotkin and Tarvainen [14] proposed the parallel two-level Schwarz methods and studied their convergence properties. Boglaev proposed a number of methods. In [15], he constructed a blocked domain decomposition algorithm. He achieved a first-order rate of convergence on both meshes. In [16], he proposed a uniform monotone iterative method on layer adapted meshes. In [17], he developed a monotone Schwarz algorithm. Boglaev and Duoba [18] designed a multi-domain decomposition algorithm to solve a singularly perturbed advection-diffusion problem with a parabolic layer. The authors achieved a first-order convergence result. Kopteva [19] and Stynes [20] propose finite element methods. Also, Newton and Picard methods were described as the numerical solver for the concerned problems by Vulkov and Zadorin in [21].

All the methods above are based on the use of non-uniform meshes and are essentially first order accurate. Due to the design of the mesh-grid and hence, that of the methods, the order of convergence is usually affected adversely by a logarithmic factor. In this article, we propose a method based on the non-standard finite difference rules of Mickens [22]. It is worth mentioning that these methods are designed on uniform grids. To the best of our awareness, this is the first time that such methods are used on elliptic singularly perturbed semilinear problems in two dimensions. These methods were used in [23, 24] for linear elliptic reaction-diffusion and reaction-convection-diffusion problems in two dimensions, respectively.

We adopt the quasilinearization approach to convert the semilinear problem into a sequence of linear problems. Then, we design a fitted operator numerical method on the converted problems. We show that the method is first order uniformly convergent in both x and y variables with respect to the perturbation parameter. Numerical experiments corroborate the theoretical results.

The rest of the article is structured as follows: In Section 2, we use the quasilinearization technique to linearize the concerned problem and present some qualitative properties of the solution and its derivatives. In Section 3, we present the proposed fitted operator finite difference method while in Section 4, we perform the convergence analysis. In Section 5, we provide some test models to show the efficiency of the presented scheme as well as to validate the theoretical result. The article ends with a brief conclusion in Section 6.

We transform the semilinear equation (1.1) using the quasilinearization approach. We choose a reliable initial guess u(0)(x,y)=u0(x,y)≡u(0). Then, we consider a truncated Taylor series expansion of f(x, y, u) about the initial approximation as follows.

We then derive the following iterates through the process by deriving the steps that involve u⁽²⁾(x, y), u⁽³⁾(x, y), and so on. Assuming that this process converges, we obtain the recurrence relations

where r is the iteration number (or iteration index) with r = 0, 1, ⋯  .

Substituting (2.2) into (1.1) results in a 2D linear singularly perturbed convection diffusion problem of the form

where

We solve the linear problem (2.3)–(2.4) using fitted operator finite difference scheme. The successive iteration of the 2D linear equations(2.3)–(2.4), with iteration function (2.5) converges to the solution of the semilinear problem (1.1)–(1.4). We take the convergence stopping criteria as

where Tol is the tolerance.

The solution of (2.3)-(2.4) enjoys the properties below [25].

Lemma 2.1. (Continuous maximum principle) Assume that ν(x, y) is sufficiently smooth function which satisfy ν(x, y) ⩾ 0, ∀ (x, y) ∈ ∂ Ω. Then Lν(x,y)⩾0,∀x∈Ω, implies that ν(x,y)⩾0,∀(x,y)∈Ω̄.

Lemma 2.2. (Uniform stability estimate) Let u(x,y) be the solution of (2.3)-(2.4) then we have

where α = min{α₁, α₂} is independent of ε.

Lemma 2.3. Let u(x,y) be the solution of (2.3)-(2.4) and a(x), b(x), z(x) be smooth functions. Then

where α₁, α₂ and C positive constant independent of ε.

Let n and m be two positive integers, we partition the domain Ω: = [0, 1] × [0, 1] into n and m equal intervals so the step sizes are h = 1/n and k = 1/m, we obtain the nodes as x_i = x₀ + ih, i = 1, …, n − 1 and y_j = y₀ + jk, j = 1, …, m − 1 where x₀ = y₀ = 0 and x_n = y_n = 1. We denote the approximation of u(x_i, y_j) at the grid points of x_i and y_j by the unknown U_ij.

We write the discrete version of (2.3)–(2.4) as

with boundary conditions of the four sides as

The denominator functions ϕij2 are given by

and

We rewrite (3.1) in five term recurrence relation as

where

We form a linear system

where U=[Ui0,…Un-1,0;U1,1…Un-1,1;…;U1,j…Un-1,m-1]T.

A is pentadiagonal matrix of size (n − 1)(m − 1) × (n − 1)(m − 1) and G is a column vector of size (n − 1)(m − 1) with their entries respectively described as follows.

and

We provide some results that we will use to prove the convergence of the proposed method. These results are similar to those presented in [24] and can be proven in a similar manner.

Lemma 3.1. (Discrete maximum principle). Let ϑ_i,j be a discrete function define on Ω satisfying ϑ_0,j ⩾ 0, ϑ_n,j ⩾ 0, j = 1(1)m − 1, ϑ_{i, 0} ⩾ 0, ϑ_i,m ⩾ 0, i = 1(1)n − 1 and Lh,kϑi,j≤0,∀i=1(1)n-1,j=1(1)m-1 then ϑ_i,j ⩾ 0 ∀ i = 0(1)n, j = 0(1)m.

Lemma 3.2. (Uniform stability estimate) if μ_i,j is any mesh function such that μ_i,j = 0 on ∂Ω^(N,M). Then

where α = min{α₁, α₂}.

The truncation error of the scheme presented in Section 3 is

Applying the bound on the solution and its derivatives in Lemma 2.3 and by Lemma 5.2 of [26], we obtain

Now, using Lemma (3.2) we have

The analysis above is the proof of the following theorem:

Theorem 4.1. Let u(x, y) be the solution of (2.3)-(2.4) and U(x, y) be the numerical approximation of u(x, y) using the scheme (3.1)-(3.2). If a₁(x, y), a₂(x, y), b(x, y) and z(x, y) are sufficiently smooth functions, then there exists a constant C independent of ε, h and k such that

In order to validate and confirm our theoretical results, we present two test examples. Since the exact solutions of our models are not available; we oblige to use the double mesh principle [27] to evaluate the maximum pointwise error and the ε-uniform error as follows

where Ui,jn,m is the discrete solution on the mesh Ω^{n, m} and Ui,j2n,2m is the discrete solution on the mesh Ω^{2n, 2m}. The corresponding rate of convergence and the ε-uniform rate of convergence are formulated as

We define the iteration stopping criterion as

Example 5.1. Boglaev [16], Consider the following singularly perturbed semilinear problem

where a1(x,y)=a2(x,y)=0.1,f(x,y,u)=u-45-u.

Example 5.2. Boglaev [15], Consider the following singularly perturbed semilinear problem

where a₁(x, y) = a₂(x, y) = 1, f(x, y, u) = 1−exp(−u).

To demonstrate the efficiency of the proposed scheme, we tabulate the maximum pointwise errors and the corresponding order of convergence. For the sake of simplicity, we considered same values of m and n as shown in Tables 1, 2. These tables indicate a first-order uniform rate of convergence that conforms to the theoretical findings in Section 4. In producing our tables, we were limited by the software used as it could not handle large matrices. Had we been able to produce the tables for larger values of n and m (say, 64, 128, 512, etc.), we would have seen that the rate of convergence is one for Example 5.2 as well.

Figures 1, 2 are plots of the numerical solution of examples 5.1 and 5.2, respectively, for n = m = 32 and different values of ε. These plots exhibit the layer behavior of the numerical solution as ε approaches zero.

We wished to compare our results with those existing in the literature however we noticed that authors that published work on this problem focused more on the number of iterations while our focus is on maximum nodal errors and rates of convergence.

In this article, we constructed a fitted operator finite difference method to solve two-dimensional semilinear singularly perturbed convection-diffusion problems. First, we converted the semilinear problems into a sequence of linear two-dimensional singularly perturbed convection-diffusion problems via the quasilinearization technique. Next, we discretized the problem using the presented non-standard numerical scheme. Then, we performed the error analysis of the method and found that it is first order uniformly convergent in both x and y variables with respect to the perturbation parameter ε. We used two test examples to illustrate the robustness of the method and to validate the theoretical findings.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Plan was prepared by JM and AA. Computations and analysis were done by OK under JM's supervision. Write up was done by all authors. All authors agree to be accountable for the content of this article. All authors contributed to the article and approved submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.