A NSFD Discretization of Two-Dimensional Singularly Perturbed Semilinear Convection-Diffusion Problems

Despite the availability of an abundant literature on singularly perturbed problems, interest toward non-linear problems has been limited. In particular, parameter-uniform methods for singularly perturbed semilinear problems are quasi-non-existent. In this article, we study a two-dimensional semilinear singularly perturbed convection-diffusion problems. Our approach requires linearization of the continuous semilinear problem using the quasilinearization technique. We then discretize the resulting linear problems in the framework of non-standard finite difference methods. A rigorous convergence analysis is conducted showing that the proposed method is first-order parameter-uniform convergent. Finally, two test examples are used to validate the theoretical findings.

1. Introduction

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

$\begin{array}{ll}-\epsilon (u_{xx}+u_{yy})+a_1(x,y)u_x+a_2(x,y)u_y=-f(x,y,u(x,y)),(x,y)\in \Omega :=(0,1)\times (0,1),& (1.1)\end{array}$

subject to boundary conditions

$\begin{array}{ll}u(x,y)=u_0(x,y),(x,y)\in \partial \Omega ,& (1.2)\end{array}$

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

$\begin{array}{ll}{a}_1(x,y)⩾{\alpha }_1>0,a_2(x,y)⩾{\alpha }_2>0,\forall (x,y)\in \bar{\Omega },& (1.3)\end{array}$

$\begin{array}{ll}{f}_u(x,y,u)⩾\beta >0,\forall (x,y)\in \bar{\Omega },& (1.4)\end{array}$

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.

2. Quasilinearization

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

$\begin{array}{ll}f(x,y,u^{(1)})=f(x,y,u^{(0)})+(u^{(1)}-u^{(0)}){\left(\frac{\partial f}{\partial u}\right)}_{(x,y,u^{(0)})}+\cdots .& (2.1)\end{array}$

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

$\begin{array}{ll}f(x,y,u^{(r+1)})=f(x,y,u^{(r)})+(u^{(r+1)}-u^{(r)}){\left(\frac{\partial f}{\partial u}\right)}_{(x,y,u^{(r)})}+\cdots ,& (2.2)\end{array}$

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

$\begin{array}{ll}\begin{array}{c}{L}u(x,y)\equiv -\epsilon (u_{xx}+u_{yy})+a_1(x,y)u_x+a_2(x,y)u_y\\ +b(x,y)u=z(x,y),(x,y)\in \bar{\Omega },\hfill \end{array}& (2.3)\end{array}$

$\begin{array}{ll}u(x,y)=u_0(x,y),& (2.4)\end{array}$

where

$\begin{array}{ll}b(x,y)=\frac{\partial f}{\partial u},\text{and}z(x,y)=f(x,y,u^{(r)})-u^{(r)}\frac{\partial f}{\partial u}.& (2.5)\end{array}$

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

$\Vert u^{r+1}-u^r\Vert <Tol,$

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}\nu (x,y)⩾0,\forall x\in \Omega$, implies that $\nu (x,y)⩾0,\forall (x,y)\in \bar{\Omega }$.

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

$\begin{array}{ll}\parallel u(x,y)\parallel ⩽{\alpha }^{-1}\parallel z\parallel ,\forall x\in \Omega ,& (2.6)\end{array}$

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

$\begin{array}{ll}\left|u^{\left(i,j\right)}\left(x,y\right)\right|⩽C\left(1+{\epsilon }^{-\left(i+j\right)}exp\left(-\frac{{\alpha }_1\left(1-x\right)}{\epsilon }\right)exp\left(-\frac{{\alpha }_2\left(1-y\right)}{\epsilon }\right)\right),\left(x,y\right)\in \bar{\Omega },& (2.7)\end{array}$

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

3. Scheme for the Problem

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

$\begin{array}{ll}\begin{array}{c}{{L}}^{h,k}(U_{i,j})\equiv -\epsilon \left[\frac{U_{i+1,j}-2U_{ij}+U_{i-1,j}}{{({\varphi }_{ij})}_h^2}\right]-\epsilon \left[\frac{U_{i,j+1}-2U_{ij}+U_{i,j-1}}{{({\varphi }_{ij})}_k^2}\right]+a_{1ij}\frac{U_{ij}-U_{i-1,j}}{h}\hfill \\ +a_{2ij}\frac{U_{ij}-U_{i,j-1}}{k}+b_{ij}{U}_{ij}=z_{ij},i=1,\dots ,n-1,j=1,\dots ,m-1,\end{array}& (3.1)\end{array}$

with boundary conditions of the four sides as

$\begin{array}{ll}{U}_{i,0}=U_{0,j}=U_{i,m}=U_{n,j}=u_{0ij},i=0,1,\dots ,n,j=0,1,\dots ,m.& (3.2)\end{array}$

The denominator functions ${\varphi }_{ij}^2$ are given by

$\begin{array}{ll}{\left({\varphi }_{ij}\right)}_h^2=\frac{\epsilon h}{a_{1ij}}\left(exp\left(\frac{a_{1ij}h}{\epsilon }\right)-1\right)=h^2+{O}\left(\frac{h^3}{\epsilon }\right),& (3.3)\end{array}$

and

$\begin{array}{ll}{\left({\varphi }_{ij}\right)}_k^2=\frac{\epsilon k}{a_{2ij}}\left(exp\left(\frac{a_{2ij}k}{\epsilon }\right)-1\right)=k^2+{O}\left(\frac{k^3}{\epsilon }\right).& (3.4)\end{array}$

We rewrite (3.1) in five term recurrence relation as

$\begin{array}{ll}\begin{array}{c}-r_l^{h+}{U}_{i+1,j}-r_l^{h-}{U}_{i-1,j}+r_l^c{U}_{ij}-r_l^{k+}{U}_{i,j+1}-r_l^{k-}{U}_{i,j-1}=z_{ij},\\ i=1(1)n-1,j=1(1)m-1,\hfill \end{array}& (3.5)\end{array}$

where

$\begin{array}{l}\begin{array}{c}{r}_l^{h+}=\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_h^2},r_l^{h-}=\left(\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_h^2}+\frac{a_{1ij}}{h}\right),r_l^{k+}=\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_k^2},r_l^{k-}=\left(\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_k^2}+\frac{a_{2ij}}{k}\right),\hfill \\ r_l^c=\frac{2\epsilon }{{\left({\varphi }_{ij}\right)}_h^2}+\frac{2\epsilon }{{\left({\varphi }_{ij}\right)}_k^2}+\frac{a_{1ij}}{h}+\frac{a_{2ij}}{k}+b_{ij},\text{and}l=\left(i-1\right)\left(m-1\right)+j.\hfill \end{array}\end{array}$

We form a linear system

$AU=G,$

where $\text{U}={\left[U_{i0},\dots U_{n-1,0};U_{1,1}\dots U_{n-1,1};\dots ;U_{1,j}\dots U_{n-1,m-1}\right]}^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.

$\begin{array}{ll}\begin{array}{c}{A}_{l,l+1}=-r_l^{k+},i=1(1)n-1,j=1(1)m-2,\hfill \\ A_{l,l-1}=-r_l^{k-},i=1(1)n-1,j=2(1)m-1,\hfill \\ A_{l,l}=r_l^c,i=1(1)n-1,j=1(1)m-1,\hfill \\ A_{l,l+(n-1)}=-r_l^{h+},i=1(1)n-2,j=1(1)m-1,\hfill \\ A_{l,l-(n-1)}=-r_l^{h-},i=2(1)n-1,j=1(1)m-1,\hfill \end{array}& (3.6)\end{array}$

and

$\begin{array}{ll}\begin{array}{c}{G}_l=z_l+r_l^{h-}\times U(0,y_1)+r_l^{k-}\times U(x_1,0),i=1,j=1,\hfill \\ G_l=z_l+r_l^{h-}\times U(0,y_j),i=1,j=2(1)m-2,\hfill \\ G_l=z_l+r_l^{h-}\times U(0,y_{m-1})+r_l^{k+}\times U(x_1,1),i=1,j=m-1,\hfill \\ G_l=z_l+r_l^{k-}\times U(x_i,0),i=2(1)n-2,j=1,\hfill \\ G_l=z_l,i=2(1)n-2,j=2(1)m-2,\hfill \\ G_l=z_l+r_l^{k+}\times U(x_i,1),i=2(1)n-2,j=m-1,\hfill \\ G_l=z_l+r_l^{h+}\times U(1,y_1)+r_l^{k-}\times U(x_{n-1},0),i=n-1,j=1,\hfill \\ G_l=z_l+r_l^{h+}\times U(1,y_j),i=n-1,j=2(1)m-2,\hfill \\ G_l=z_l+r_l^{h+}\times U(1,y_{m-1})+r_l^{k+}\times U(x_{n-1},1),i=n-1,j=m-1.\hfill \end{array}\hfill & (3.7)\end{array}$

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 ${{L}}^{h,k}{\vartheta }_{i,j}\le 0,\forall 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

$\begin{array}{ll}\left|{\mu }_{l,s}\right|⩽\frac{1}{\alpha }{\displaystyle \underset{1\le j\le n-1,1\le j\le m-1}{max}}\left|{{L}}^{h,k}{\mu }_{i,j}\right|\forall l=0(1)n,s=0(1)m,& (3.8)\end{array}$

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

4. Convergence Analysis

The truncation error of the scheme presented in Section 3 is

$\begin{array}{l}\begin{array}{c}{{L}}^{h,k}{\left(u-U\right)}_{ij}=\left({L}-{{L}}^{h,k}\right)u_{ij}\hfill \\ =-\epsilon {\left(u_{xx}\right)}_{ij}-\epsilon {\left(u_{yy}\right)}_{ij}+a_{1ij}{\left(u_x\right)}_{ij}+a_{2ij}{\left(u_y\right)}_{ij}\hfill \\ +\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_h^2}\left(u_{i+1,j}-2u_{i}j+u_{i-1,j}\right)+\frac{\epsilon }{{\left({\varphi }_{ij}\right)}_k^2}\left(u_{i,j+1}-2u_{i}j+u_{i,j-1}\right)\hfill \\ -\frac{a_{1ij}}{h}\left(u_{ij}-u_{i-1,j}\right)-\frac{a_{2ij}}{k}\left(u_{ij}-u_{i,j-1}\right)\hfill \\ =-\epsilon {\left(u_{xx}\right)}_{ij}-\epsilon {\left(u_{yy}\right)}_{ij}+\left(\frac{\epsilon }{h^2}-\frac{a_{1ij}}{h}+\frac{a_{1ij}^2}{\epsilon }-\frac{a_{1ij}^3}{{\epsilon }^2}\right)\hfill \\ \times \left(h^2{\left(u_{xx}\right)}_{ij}+\frac{h^4}{12}{\left(u_{xxxx}\right)}_{ij}{\eta }_{1ij}\right)+\left(\frac{\epsilon }{k^2}-\frac{a_{2ij}}{k}+\frac{a_{2ij}^2}{\epsilon }-\frac{a_{2ij}^3}{{\epsilon }^2}\right)\hfill \\ \times \left(k^2{\left(u_{yy}\right)}_{ij}+\frac{h^4}{12}{\left(u_{yyyy}\right)}_{ij}{\eta }_{2ij}\right)+\left(\frac{a_{1ij}h}{2}{\left(u_{xx}\right)}_{ij}+\frac{a_{1ij}{h}^2}{6}{\left(u_{xxx}\right)}_{ij}\right)\hfill \\ +\left(\frac{a_{2ij}k}{2}{\left(u_{yy}\right)}_{ij}+\frac{a_{2ij}{k}^2}{6}{\left(u_{yyy}\right)}_{ij}\right)\hfill \\ \text{where }{\eta }_{1ij}\in \left(u_{i+1,j},u_{i-1,j}\right),{\eta }_{2ij}\in \left(u_{i,j+1},u_{i,j-1}\right)\hfill \\ =-\frac{a_{1ij}h}{2}{\left(u_{xx}\right)}_{ij}+\left(\frac{a_{1ij}^2}{\epsilon }{\left(u_{xx}\right)}_{ij}-\frac{a_{1ij}^2}{6}{\left(u_{xxx}\right)}_{ij}+\frac{\epsilon }{12}{\left(u_{xxxx}\right)}_{ij}{\eta }_{1ij}\right)h^2\hfill \\ +\left(\frac{a_{1ij}}{12}{\left(u_{xxxx}\right)}_{ij}{\eta }_{1ij}-\frac{a_{1ij}^3}{{\epsilon }^3}{\left(u_{xx}\right)}_{ij}{\eta }_{1ij}\right)h^3+\left(\frac{a_{1ij}^2}{12\epsilon }{\left(u_{xxx}\right)}_{ij}\right)h^4\hfill \\ -\left(\frac{a_{1ij}^3}{12{\epsilon }^2}{\left(u_{xxxx}\right)}_{ij}{\eta }_{1ij}\right)h^5-\frac{a_{2ij}k}{2}{\left(u_{yy}\right)}_{ij}\hfill \\ +\left(\frac{a_{2ij}^2}{\epsilon }{\left(u_{yy}\right)}_{ij}-\frac{a_{2ij}^2}{6}{\left(u_{yyy}\right)}_{ij}+\frac{\epsilon }{12}{\left(u_{yyyy}\right)}_{ij}{\eta }_{2ij}\right)k^2\hfill \\ +\left(\frac{a_{2ij}}{12}{\left(u_{yyyy}\right)}_{ij}{\eta }_{2ij}-\frac{a_{2ij}^3}{{\epsilon }^3}{\left(u_{yy}\right)}_{ij}{\eta }_{2ij}\right)k^3+\left(\frac{a_{2ij}^2}{12\epsilon }{\left(u_{yyy}\right)}_{ij}\right)k^4\hfill \\ -\left(\frac{a_{2ij}^3}{12{\epsilon }^2}{\left(u_{yyyy}\right)}_{ij}{\eta }_{2ij}\right)k^5\hfill \end{array}\end{array}$

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

$\begin{array}{l}\left|{{L}}^{h,k}{(u-U)}_{ij}\right|⩽C(h+k).\end{array}$

Now, using Lemma (3.2) we have

$\begin{array}{ll}{\displaystyle \underset{0\le i\le n,0\le j\le m}{max}}\left|{(u-U)}_{ij}\right|⩽C(h+k).& (4.1)\end{array}$

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

$\begin{array}{ll}{\displaystyle \underset{0\le i\le n,0\le j\le m}{max}}\left|{(u-U)}_{ij}\right|⩽C(h+k).& (4.2)\end{array}$

5. Numerical Results

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

$\begin{array}{ll}{E}_{\epsilon }^{n,m}={\displaystyle \underset{(x_i,y_j)\in {\bar{\Omega }}^{n,m}}{max}}|U_{i,j}^{n,m}-U_{i,j}^{2n,2m}|,E^{n,m}={\displaystyle \underset{\epsilon }{max}}{E}_{\epsilon }^{n,m},& (5.1)\end{array}$

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

$\begin{array}{ll}{P}_{\epsilon }^{n,m}={log}_2\left(\frac{E_{\epsilon }^{n,m}}{E_{\epsilon }^{2n,2m}}\right),P^{n,m}={\displaystyle \underset{\epsilon }{max}}{P}_{\epsilon }^{n,m}.& (5.2)\end{array}$

We define the iteration stopping criterion as

$\begin{array}{ll}\parallel U^{(r+1)}-U^{(r)}\parallel ⩽10^{-8},r=1,2,\dots .& (5.3)\end{array}$

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

$\begin{array}{l}\begin{array}{c}-\epsilon (u_{xx}+u_{yy})+a_1(x,y)u_x+a_2(x,y)u_y+f(x,y,u)=0,\\ (x,y)\in \Omega :={(0,1)}^2,\end{array}\end{array}$

$u=1on\partial \Omega ,$

where $a_1(x,y)=a_2(x,y)=0.1,f(x,y,u)=\frac{u-4}{5-u}.$

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

$\begin{array}{l}\begin{array}{c}-\epsilon (u_{xx}+u_{yy})+a_1(x,y)u_x+a_2(x,y)u_y+f(x,y,u)=0,\\ (x,y)\in \Omega :={(0,1)}^2,\end{array}\end{array}$

$u=1on\partial \Omega ,$

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.

TABLE 1

Table 1. Maximum pointwise errors and rate of convergence for Example 5.1 when n = m = {4, 8, 16, 32, 64, 128, 512, 1, 024}.

TABLE 2

Table 2. Maximum pointwise errors and rate of convergence for Example 5.2 when n = m = {4, 8, 16, 32, 64, 128, 512, 1, 024}.

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.

FIGURE 1

Figure 1. Plots of the approximate solution of Example 5.1 with n = m = 32. (A) ε = 1. (B) ε = 10⁻¹. (C) ε = 10⁻². (D) ε = 10⁻⁴.

FIGURE 2

Figure 2. Plots of the approximate solution of Example 5.2 with n = m = 32. (A) ε = 1. (B) ε = 10⁻¹. (C) ε = 10⁻². (D) ε = 10⁻⁴.

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.

6. Conclusion

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.

Data Availability Statement

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

Author Contributions

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.

Conflict of Interest

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.

Publisher's Note

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.

Acknowledgments

The research of JM and AA was supported by the National Research Foundation of South Africa. The authors wish to thank the reviewers for their valuable comments and suggestions that helped improve the quality of this paper.

References

1. Prandtl L. Über Flussigkeitsbewegung bei sehr kleiner Reibung. In: Proc. III Intern. Congr Math. Heildeberg (1904).

2. He X, Wang K. New finite difference methods for singularly perturbed convection-diffusion equations. Taiwanese J Math. (2018) 22:949–78. doi: 10.11650/tjm/171002

CrossRef Full Text | Google Scholar

3. Kadalbajoo MK, Patidar KC. Spline approximation method for solving self-adjoint singular perturbation problems on nonuniform grids. J Comput Anal Appl. (2003) 5:425–51. doi: 10.1155/S0161171203112094

CrossRef Full Text | Google Scholar

4. Nhan TA, Vulanovic R. The Bakhvalov mesh: a complete finite-difference analysis of two-dimensional singularly perturbed convection-diffusion problems. Num Algorith. (2021) 87:203–21. doi: 10.1007/s11075-020-00964-z

CrossRef Full Text | Google Scholar

5. Rao SCS, Chaturvedi AK. Parameter-uniform numerical method for a two-dimensional singularly perturbed convection-reaction-diffusion problem with interior and boundary layers. Math Comput Simul. (2021) 187:656–86. doi: 10.1016/j.matcom.2021.03.016

CrossRef Full Text | Google Scholar

6. Roos HG, Line T. Gradient recovery for singularly perturbed boundary value problems II: Two-dimensional convection-diffusion. Math Models Methods Appl Sci. (2001) 11:1169–79. doi: 10.1142/S0218202501001288

CrossRef Full Text | Google Scholar

7. Stynes M, Tobiska L. Necessary L2-uniform convergence conditions for difference schemes for two-dimensional convection-diffusion problems. Comput Math Appl. (1995) 29:45–53. doi: 10.1016/0898-1221(94)00237-F

CrossRef Full Text | Google Scholar

8. Alexandrov VM, Pozharskii DA, Tielking JT. Three-dimensional contact problems. Solid Mech Appl. (2002) 55:B73–4. doi: 10.1115/1.1483357

CrossRef Full Text | Google Scholar

9. Beebe NH. A Complete Bibliography of Lecture Notes in Computational Science and Engineering. (2017).

10. Gryanik V, Hartmann J. Towards a Solution of the Closure Problem for Convective Atmospheric Boundary-Layer Turbulence. Trieste: The Abdus Salam International Centre for Theoretical Physics (2017).

11. Gutkowski W, Kowalewski TA. Mechanics of the 21st century. In: Proceedings of the 21st International Congress of Theoretical and Applied Mechanics. Warsaw (2006). doi: 10.1007/1-4020-3559-4

CrossRef Full Text | Google Scholar

12. Sivertsen OI, Andrianov I. Virtual testing of mechanical systems: theories and techniques. Adv Eng. (2002) 55:B66–7. doi: 10.1115/1.1483349

CrossRef Full Text | Google Scholar

13. Ventsel E, Krauthammer T, Carrera E. Thin plates and shells: theory, analysis, applications. Appl Mech Rev. (2002) 55:B72–3. doi: 10.1115/1.1483356

CrossRef Full Text | Google Scholar

14. Sirotkin V, Tarvainen P. Parallel Schwarz methods for convection-dominated semilinear diffusion problems. J Comput Appl Math. (2002) 145:189–211. doi: 10.1016/S0377-0427(01)00575-1

CrossRef Full Text | Google Scholar

15. Boglaev I. A block monotone domain decomposition algorithm for semilinear convection-diffusion problem. J Comput Appl Math. (2005) 173:259–77. doi: 10.1016/j.cam.2004.03.011

CrossRef Full Text | Google Scholar

16. Boglaev I. Uniform convergence of monotone iterative methods for semilinear singularly perturbed problems of elliptic and parabolic types. Electron Trans Num Anal. (2005) 20:86–103.

17. Boglaev I. A monotone Schwarz algorithm for a semilinear convection-diffusion problem. J Num Math. (2004) 12:169–91. doi: 10.1515/1569395041931455

CrossRef Full Text | Google Scholar

18. Boglaev I, Duoba V. Domain decomposition for an advection-diffusion problem with parabolic layers. Appl Math Comput. (2003) 146:27–53. doi: 10.1016/S0096-3003(02)00514-3

CrossRef Full Text | Google Scholar

19. Kopteva N. How accurate is the streamline-diffusion FEM inside characteristic (boundary and interior) layers. Comput Methods Appl Mech Eng. (2004) 193:4875–89. doi: 10.1016/j.cma.2004.05.008

CrossRef Full Text | Google Scholar

20. Stynes M. Convection-diffusion-reaction problems, SDFEM/SUPG and a priori meshes. Int J Comput Sci Math. (2007) 1:412–31. doi: 10.1504/IJCSM.2007.016543

CrossRef Full Text | Google Scholar

21. Vulkov L, Zadorin AI. Two-grid algorithms for the solution of 2D semilinear singularly-perturbed convection-diffusion equations using an exponential finite difference scheme. In: AIP Conference Proceedings. American Institute of Physics (2009). p. 371–379. doi: 10.1063/1.3265351

CrossRef Full Text | Google Scholar

22. Mickens RE. Nonstandard Finite Difference Models of Differential Equations. World Scientific (1994). doi: 10.1142/2081

CrossRef Full Text | Google Scholar

23. Munyakazi JB, Patidar KC. Higher order numerical methods for singularly perturbed elliptic problems. Neural Parallel Sci Comput. (2010) 18:75–88.

24. Munyakazi JB, Patidar KC. Novel fitted operator finite difference methods for singularly perturbed elliptic convection-diffusion problems in two dimensions. J Diff Equat Appl. (2012) 18:799–813. doi: 10.1080/10236198.2010.513330

CrossRef Full Text | Google Scholar

25. Miller JJ, O'riordan E, Shishkin GI. Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions. World Scientific (1996). doi: 10.1142/2933

CrossRef Full Text | Google Scholar

26. Kadalbajoo MK, Patidar KC. ε-uniformly convergent fitted mesh finite difference methods for general singular perturbation problems. Appl Math Comput. (2006) 179:248–66. doi: 10.1016/j.amc.2005.11.096

CrossRef Full Text | Google Scholar

27. Doolan EP, Miller JJ, Schilders WH. Uniform Numerical Methods for Problems With Initial and Boundary Layers. Eindhoven: Boole Press (1980).