A Uniformly Convergent Scheme for Singularly Perturbed Unsteady Reaction–Diffusion Problems

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The majority of real-life events in applied science and engineering problems are particular situations of mathematical models of differential equations that depend on parameters. Singularly perturbed problems (SPPs) are those whose highest order derivative is multiplied by a small parameter $\epsilon ,0<\epsilon \ll 1$ and whose solutions or the derivatives approach the discontinuous limit as $\epsilon$ approaches zero as discussed by Roos [1], Boffi [2], and Sigamani et al. [3]. The most interesting feature about SPPs is that the solutions contain boundary and/or interior layers when the value of the perturbation parameter $\epsilon$ approaches zero. These features make the classical numerical methods incompatible with accurately resolving these types of problems. This emphasises the need for parameter-independent numerical methods, as explained in Miller, O’riordan, and Shishkin [4]; Linß [5]; Rao, Kumar, and Singh [6]; and Kumar et al. [7]. Consequently, parameter-uniform convergent algorithms have been developed. The fitted operator approach, the fitted mesh approach, the domain decomposition approach, the grid equidistribution approach, and the discrete Schwarz waveform relaxation method are some of the examples of such parameter-independent uniformly convergent numerical methods.

Reaction–diffusion equations (RDEs) are a set of partial differential equations (PDEs) that describe the variations in species, densities, or concentrations of at least two interacting variables in both space and time. For instance, in chemistry, physics, biology, and fluid dynamics, RDEs can describe a vast array of very complex systems according to Ji and Fenton [8]. Lately, singularly perturbed (SP) time-dependent reaction–diffusion problems (RDPs) have been active research areas for the aforementioned disciplines.

Many researchers have treated a class of SPPs using finite difference method (FDM) on nonuniform meshes for the spatial domain and uniform meshes for the temporal domain discretization. For example, S. Kumar and M. Kumar [9] investigated a coupled system of first-order singularly perturbed quasilinear differential equations. In Kadalbajoo and Awasthi [10], a uniformly convergent FDM is designed to treat a SP time-dependent convection–diffusion problem (CDP) in one space dimension. A higher-order FDM with the pattern expanding method is applied for solving a coupled system of SPP first-order ordinary differential equations (see in Rao and Kumar [11]). The authors further proved that the proposed method is second-order accurate in the spatial direction and first-order in the temporal direction. In Kumar and Vigo-Aguiar [12], a first-order uniformly convergent scheme is used to treat a class of SPP degenerate parabolic CDPs on a rectangular region. Clavero and Gracia [13] developed a numerical approach for solving one-dimensional parabolic-type SP RDPs. Furthermore, the authors used the uniform stability analysis to prove that the formulated scheme is parameter-uniform convergent.

Due to the simplicity of the analysis of uniform meshes, other scholars have treated SP parabolic-type RDPs using uniform meshes for the purpose of domain discretization in the spatial and temporal directions. The following are a few to mention: Kumar and Ramos [14] investigated a parameter-independent numerical method based on equidistributed meshes for a class of SP parabolic-type RDPs with Robin boundary conditions. The considered problem domain is discretized using the modified Euler method for the temporal direction on uniform meshes and a central difference method for the spatial direction via the equidistribution of a suitably chosen monitor function. The various error analyses and discussions have proven that the proposed method is a parameter-uniform convergence method with second-order accuracy in the space and first-order accuracy in the time directions. Munyakazi and Patidar [15] employed a fitted numerical technique to treat time-dependent, SP RDPs. The continuous domain of the considered problem is discretized using the backward Euler method in the temporal direction and the nonstandard finite difference method (NSFDM) in the spatial directions on uniform meshes. The convergence analysis of the scheme proved that the suggested approach was found to be parameter-uniform convergent.

Formulating and examining a second-order accurate uniformly convergent numerical scheme for solving SP unsteady RDPs on uniform meshes are the ultimate goal of the present work. The Crank–Nicolson method for time derivatives and NSFDM for space derivative are used to discretize the proposed problem domain. The formulated scheme is second-order accurate parameter-uniformly convergent, according to the paper’s extensive analyses.

Notations. Throughout this work, $C$ and its subscripts denote generic positive constants independent of the perturbation parameter $\epsilon$ and step sizes $h$ and $k$. Also, $.$ denotes the standard supremum norm, defined as $f={\sup }_{x,t\in {\overline{D}}_{xt}}fx,t$ for a function $f$ defined on some domain ${\overline{\mathcal{D}}}_{xt}$. The function $\Psi x,t$ is a solution of the proposed problem and ${\Psi }^{i,j}x,t$ is its $i^{\text{th}}$ spatial derivative at $j^{\text{th}}$ time level. $D_l$ and $D_r$ are left and right Dirichlet boundary conditions, respectively.

2. Statement of the Problem

In this work, we consider a singularly perturbed parabolic-type unsteady RDP with initial- and interval-boundary conditions of the following form: $$\begin{array}{}\text{(1)}& \begin{array}{c}{\mathcal{L}}_{\epsilon }\Psi :\equiv {\Psi }_{t}x,t-\epsilon {\Psi }_{xx}x,t+ax,t\Psi x,t=fx,t,\forall x,t\in {\mathcal{D}}_{xt}\hfill \\ \Psi x,0=S_{0}x\forall x\in {\mathcal{D}}_x\hfill \\ \Psi 0,t=S_{l}t\forall t\in {\mathcal{D}}_t\hfill \\ \Psi 1,t=S_{r}1,t\forall t\in {\mathcal{D}}_t\hfill \end{array}\end{array}$$where ${\mathcal{L}}_{\epsilon }\Psi x,t={\Psi }_{t}x,t-\epsilon {\Psi }_{xx}x,t+ax,t\Psi x,t$ is the differential operator, $\epsilon 0<\epsilon \ll 1$ is the perturbation parameter, ${\mathcal{D}}_x=0,1$, ${\mathcal{D}}_t=0,T$, and ${\mathcal{D}}_{xt}={\mathcal{D}}_x\times {\mathcal{D}}_t$. The functions $fx,t$ and $ax,t$ over ${\overline{\mathcal{D}}}_{xt}$, $S_{0}x$ over ${\overline{\mathcal{D}}}_x$, $S_{l}t$, and $S_{r}1,t$) over ${\overline{\mathcal{D}}}_t$ are sufficiently smooth, and the bounded function $ax,t$ satisfies $ax,t\ge A^{\ast }>0$.

2.1. Bounds on the Solution and Its Derivatives

Checking the bounds of the analytical solution $\Psi x,t$ of Problem (1) and its derivatives using Hȯlder’s continuity with appropriate compatibility conditions as described in [16] is the primary focus of this section. Additionally, we assumed that $0<\epsilon \ll 1$ and $a\ge A^{\ast }>0\forall x,y\in {\overline{\mathcal{D}}}_{xt}$. In this case, the data $a,f\in C^{4,2}{\overline{\mathcal{D}}}_{xt}$ are sufficiently smooth, and the exact solution $\Psi x,t\in C^{4,2}{\overline{D}}_{xt}$ is also sufficiently smooth. The required compatibility conditions at the corner points $0,0$ and $1,0$ are defined as $S_{0}0,0=S_{l}0$ and $S_{0}1,0=S_{r}1$ and satisfies $$\begin{array}{}\text{(2)}& \begin{array}{c}\frac{\partial S_{l}0}{\partial t}-\epsilon \frac{{\partial }^2{S}_{0}0}{\partial x^2}+a0,0S_{0}0=f0,0\hfill \\ \frac{\partial S_{r}1}{\partial t}-\epsilon \frac{{\partial }^2{S}_{0}1}{\partial x^2}+a1,0S_{0}1=f1,0\hfill \end{array}\end{array}$$

Under the above suitable conditions on the initial data and the initial-boundary value, Problem (1) has unique solution $\Psi x,t$ which exhibits twin boundary layers of width $O\sqrt{\epsilon }$ at $x=0$ and $x=1$. According to the approach in [7], these conditions guarantee existence of a constant $C$ independent of the parameter $\epsilon$ such that $$\begin{array}{}\text{(3)}& \begin{array}{c}\Psi x,t-S_{0}x\le Ct,x,t\in {\mathcal{D}}_x\hfill \\ \Psi x,t\le C,x,t\in {\overline{\mathcal{D}}}_{xt}\hfill \end{array}\end{array}$$

Lemma 2.1 (continuous maximum principle).

Let $\varphi x,t$ be any sufficiently smooth function satisfying $\varphi x,t\ge 0,\forall x,t\in \partial {\mathcal{D}}_{xt}$ and ${\mathcal{L}}_{\epsilon }\varphi x,t\ge 0,\forall x,t\in {\mathcal{D}}_{xt}$. Then, $\varphi x,t\ge 0,\forall x,t\in {\overline{\mathcal{D}}}_{xt}$ [17].

Proof 2.1.

Suppose that $\eta ,\delta \in {\mathcal{D}}_{xt}$ such that $\varphi \eta ,\delta ={\text{min}}_{x,t\in {\overline{\mathcal{D}}}_{xt}}\varphi x,t<0$.

From this condition, it is seen that ${\varphi }_t\eta ,\delta =0,$${\varphi }_x\eta ,\delta =0$, and ${\varphi }_{xx}\eta ,\delta \ge 0$. Having these, all the differential operator $\mathcal{L}$ can also be defined in the following manner: $$\begin{array}{}& {\mathcal{L}}_{\epsilon }\varphi \eta ,\delta ={\varphi }_t\eta ,\delta -\epsilon {\varphi }_{xx}\eta ,\delta +a\eta ,\delta \varphi \eta ,\delta =-\epsilon {\varphi }_{xx}\eta ,\delta +a\eta ,\delta \varphi \eta ,\delta \le 0,\text{since},{\varphi }_{xx}\ge 0,\varphi \le 0\text{and}a>0\text{at}\eta ,\delta \end{array}$$

This implies ${\mathcal{L}}_{\epsilon }\varphi \eta ,\delta \le 0$ which contradicts to our assumption $L_{\epsilon }\psi \eta ,\delta \ge 0$.

Therefore, we can conclude that the minimum of $\varphi x,t$ is nonnegative.

Lemma 2.2 (stability estimate).

Let $\Psi x,t\in C^{2,1}{\overline{\mathcal{D}}}_{xt}$ be the solution of the continuous problem of Equation (1). Then, it satisfies the bound [18]: $$\begin{array}{}& \Psi x,t\le \max S_{0}x,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\end{array}$$

Proof 2.2.

To prove this lemma, we defined two barrier functions ${\Omega }^{\pm }$ as follows: $$\begin{array}{}& {\Omega }^{\pm }x,t=\max S_{0}x,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm \Psi x,t\end{array}$$

When $t=0$, $$\begin{array}{}& {\Omega }^{\pm }x,0=\max S_{0}x,S_{l}0,S_{r}0+\frac{1}{A^{\ast }}f\pm \Psi x,0=\max S_{0}x,S_{l}0,S_{r}0+\frac{1}{A^{\ast }}f\pm S_{0}x\ge 0\end{array}$$

When $x=0$, $$\begin{array}{}& {\Omega }^{\pm }0,t=\max S_{0}0,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm \Psi 0,t=\max S_{0}0,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm S_{l}t\ge 0\end{array}$$

When $x=1$, $$\begin{array}{}& {\Omega }^{\pm }1,t=\max S_{0}1,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm \Psi 1,t=\max S_{0}0,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm S_{r}t\ge 0\end{array}$$

Combining all these on the domain ${\mathcal{D}}_{xt}$ results in $$\begin{array}{}& {\mathcal{L}}_{\epsilon }{\Omega }^{\pm }x,t=\frac{\partial {\Omega }^{\pm }x,t}{\partial t}-\epsilon \frac{{\partial }^2{\Omega }^{\pm }x,t}{\partial x^2}+ax,t{\Omega }^{\pm }x,t=ax,t\frac{1}{A^{\ast }}f+bx,t\times \max S_{0}x,S_{l}t,S_{r}t,+\frac{1}{A^{\ast }}f\pm L_{\epsilon }\Psi x,t=bx,t\times \max S_{0}x,S_{l}t,S_{r}t+\frac{1}{A^{\ast }}f\pm fx,t\ge 0\end{array}$$

Hence, from Lemma (2.1), it follows that ${\Omega }^{\pm }x,t\ge 0,\forall x,t\in {\overline{\mathcal{D}}}_{xt}$.

Lemma 2.3.

Let $\Psi x,t$ be the solution of Problem (1). Then, we have $$\begin{array}{}\text{(4)}& {\Psi }^{i,j}x,t\le C1+{\epsilon }^{-i/2}\exp -\sqrt{\frac{A^{\ast }}{\epsilon }}x+\exp -\sqrt{\frac{A^{\ast }}{\epsilon }}1-x\end{array}$$where $0\le i+2j\le 4$.

The proof of this lemma is found in [8].

3. Domain Discretization

The main concern of this section is to discretize the problem domain using the concept of semidiscretization for the temporal and spatial directions separately. Finally, we combined the two discretized domains together and found the fully discrete solution of the proposed problem defined in Equation (1).

3.1. Temporal Discretization

Divide the domain $0,T$ into equal number of $J$ mesh intervals of length $k$: $$\begin{array}{}\text{(5)}& {\overline{\mathcal{D}}}_t^J=t_j:t_j=jk,j=0,1,2,3,\cdots J\text{and}k=\frac{T}{J}\end{array}$$

Following the approach in [19] and using the Crank–Nicolson method, the problem in Equation (1) is discretized as follows: $$\begin{array}{}\text{(6)}& \begin{array}{c}\Psi x,0=S_{0}x,\forall x\in {\mathcal{D}}_x\hfill \\ \frac{{\Psi }^{j+1}-{\Psi }^j}{k}=\epsilon \frac{{\Psi }_{xx}^{j+1}+{\Psi }_{xx}^j}{2}-\frac{{a\Psi }^{j+1}+{a\Psi }^j}{2}+\frac{f^{j+1}+f^j}{2}\hfill \\ \Psi 0,t_{j+1}=S_l{t}_{j+1},0\le j\le J-1\hfill \\ \Psi 1,t_{j+1}=S_{r}1,t_{j+1},0\le j\le J-1\hfill \end{array}\end{array}$$

Again, separating the $j^{\text{th}}$ and ${j+1}^{\text{th}}$ time-level equations and rearranging gives $$\begin{array}{}\text{(7)}& \begin{array}{c}\Psi x,0=S_{0}x,\forall x\in {\mathcal{D}}_x\hfill \\ -\epsilon {\Psi }_{xx}^{j+1}+2/k+a^{j+1}{\Psi }^{j+1}=\epsilon {\Psi }_{xx}^j+2/k-a^j{\Psi }^j+f^{j+1}+f^j\hfill \\ \Psi 0,t_{j+1}=S_l{t}_{j+1},0\le j\le J-1\hfill \\ \Psi 1,t_{j+1}=S_{r}1,t_{j+1},0\le j\le J-1\hfill \end{array}\end{array}$$

To check the semidiscrete maximum principle and stability estimate of the temporal semidiscretization of the proposed scheme, we first rewrite Equation (7) at ${j+1}^{\text{th}}$ time level in terms of spatial variables as follows: $$\begin{array}{}\text{(8)}& \begin{array}{c}{\mathcal{L}}^{\ast }\Gamma :\equiv -\epsilon {\Gamma }^{\prime \prime }x+B^{\ast }x\Gamma x=g^{\ast }x,\forall x\in {\overline{\mathcal{D}}}_x\hfill \\ \Gamma 0=S_l{t}_{j+1},\Gamma 1=S_r{t}_{j+1},\forall t_{j+1}\in {\overline{\mathcal{D}}}_t^J\hfill \end{array}\end{array}$$where ${\mathcal{L}}^{\ast }$ is a discrete operator, $\Gamma x={\Psi }^{j+1}x,B^{\ast }x=2/k+a^{j+1}x$, and $g^{\ast }x=\epsilon {\Psi }_{xx}^j+2/k-a^j{\Psi }^j+f^{j+1}+f^j$.

Lemma 3.1 (semidiscrete maximum principle).

If ${\varphi }^{j+1}x\in C^2{\overline{\mathcal{D}}}_x$ is the solution of semidiscrete Problem (7) such that ${\varphi }^{j+1}0\ge 0,{\varphi }^{j+1}1\ge 0$ on $\partial {\overline{\mathcal{D}}}_x$ and ${\mathcal{L}}^{\ast }{\varphi }^{j+1}x\ge 0,\forall x\in {\mathcal{D}}_x$, then ${\varphi }^{j+1}x\ge 0,\forall x\in {\overline{\mathcal{D}}}_x$.

Proof 3.1.

Let $\eta \in {\overline{\mathcal{D}}}_x$ such that ${\varphi }^{j+1}\eta =\text{mi}{\text{n}}_{x\in {\overline{D}}_x}{\varphi }^{j+1}x<0$ and since ${\varphi }^{j+1}0\ge 0,{\varphi }^{j+1}1\ge 0$ then $\eta \notin 0,1$, that is, $\eta \in {\mathcal{D}}_x$. Moreover, since ${\varphi }_x^{j+1}\eta =0$ and ${\varphi }_{xx}^{j+1}\eta \ge 0$, we have that $$\begin{array}{}& {\mathcal{L}}^{\ast }{\varphi }^{j+1}\eta =-\epsilon {\varphi }_{xx}^{j+1}\eta +B^{\ast }\eta {\varphi }^{j+1}\eta \le 0\end{array}$$

This shows ${\mathcal{L}}^{\ast }{\varphi }^{j+1}x\le 0$, which contradicts to our assumption ${\mathcal{L}}^{\ast }{\varphi }^{j+1}x\ge 0$.

Consequently, we conclude that the minimum of ${\varphi }^{j+1}x$ is nonnegative.

Lemma 3.2 (stability estimate).

Let $\Gamma x$ be the solution of the semidiscrete problem (10). Then, we have $$\begin{array}{}\text{(9)}& \Gamma x\ge \frac{1}{\alpha }{\mathcal{L}}^{\ast }\Gamma x+{\Gamma x}_{{\mathcal{D}}_x},B^{\ast }x\ge \alpha >0\end{array}$$

The reader may refer the proof of this lemma in Lemma 2.2 of Section 2.

The local truncation error of semidiscrete Problem (7) is given by $e^{j+1}\equiv \Psi x,t_j-{\Psi }^{j+1}$, where ${\Psi }^{j+1}$ is the computed solution of the proposed boundary value problem. Moreover, this error measures the contribution of each time step to the global error of the temporal semidiscretization.

Lemma 3.3 (local error estimate).

Suppose that $$\begin{array}{}& \frac{{\partial }^i\Psi x,t}{\partial t^i}\le C,\forall x,t\in {\overline{\mathcal{D}}}_x\text{and}0\le i\le 2\end{array}$$

Then, the local error estimate in the temporal direction is given by $$\begin{array}{}\text{(10)}& {e^{j+1}}_{\infty }\le C_1{k}^3\end{array}$$

Proof 3.2.

To prove this, we follow the approach in [19] and use Taylor series expansion of the functions ${\Psi }^j$ and ${\Psi }^{j+1}$ about $x,t_{j+1/2}$: $$\begin{array}{}\text{(11)}& {\Psi }^{j+1}=\Psi x,t_{j+1/2}+\frac{k}{2}{\Psi }_{t}x,t_{j+1/2}+\frac{k^2}{8}{\Psi }_{tt}x,t_{j+1/2}+O{k}^3\text{(12)}& {\Psi }^j=\Psi x,t_{j+1/2}-\frac{k}{2}{\Psi }_{t}x,t_{j+1/2}+\frac{k^2}{8}{\Psi }_{tt}x,t_{j+1/2}+O{k}^3\end{array}$$

Subtracting Equation (11) from Equation (12) gives $$\begin{array}{}\text{(13)}& \frac{{\Psi }^{j+1}-{\Psi }^j}{k}={\Psi }_{t}x,t_{j+1/2}+O{k}^2\end{array}$$

Substituting Equation (13) in Equation (1), we get $$\begin{array}{}\text{(14)}& \frac{{\Psi }^{j+1}-{\Psi }^j}{k}=\epsilon {\Psi }_{xx}x,t_j+\frac{k}{2}-a\Psi x,t_j+\frac{k}{2}+fx,t_m+\frac{k}{2}+O{k}^2\end{array}$$where $fx,t_{j+1/2}=fx,t_{j+1}+fx,t_{j}fx,t_{j+1}+fx,t_j/2\Psi x,t_{j+1/2}=\Psi x,t_{j+1}+\Psi x,t_j\Psi x,t_{j+1}+\Psi x,t_j/2$ and $ax,t_{j+1/2}=ax,t_{j+1}+ax,t_{j}ax,t_{j+1}+ax,t_j/2$.

The solution to the local error is given by $$\begin{array}{}\text{(15)}& \begin{array}{c}{\mathcal{L}}^{\ast }{e}^{j+1}=O{k}^3\hfill \\ e^{j+1}0=0,e^{j+1}1=0\hfill \end{array}\end{array}$$

Thus, using the maximum principle in Lemma (3.1), the required result is satisfied.

Hence, ${e^{j+1}}_{\infty }\le C_1{k}^3$ as required.

Lemma 3.4 (global error estimate).

Based on the assumption of Lemma (3.3), the global error estimate of the semidiscretization at $j^{\text{th}}$ time level is given by $$\begin{array}{}\text{(16)}& {E^j}_{\infty }\le C_2{k}^2\end{array}$$

Proof 3.3.

Since the global error estimate of the $j^{\text{th}}$ time level is given by $$\begin{array}{}\text{(17)}& {E^j}_{\infty }={\sum _{s=1}^j{T}^s}_{\infty },\text{since}j\le T/k\le {T^1}_{\infty }+{T^2}_{\infty }+\cdots +{T^s}_{\infty }\le C_{1}j{k}^3,\text{using}\text{Lemma}3.3\le C_{1}T{k}^2,\text{since}jk\le T\le C_2{k}^2,C_2=C_{1}T\end{array}$$where $C_2$ is a constant independent of the parameter $\epsilon$ and time step length $k$.

From this, it is seen that the Crank–Nicolson method is of second-order accurate.

3.2. Spatial Discretization

Let the domain $0,1$ be divided into equal number of $K$ mesh points each of width $h$. Then, the domain ${\overline{\mathcal{D}}}_x$ is discretized as follows: $$\begin{array}{}\text{(18)}& {\overline{\mathcal{D}}}_x^K=x_i=ih,x_K=1,i=1,2,3,\cdots K-1\text{and}h=\frac{1}{K}\end{array}$$

Using this discretized domain, the NSFDM is employed in solving Problem (7) in the subsequent sections.

3.2.1. Nonstandard Finite Difference Scheme

To construct the exact finite difference approach, we considered the following constant coefficient homogeneous differential equations, corresponding to Problem (7) as in [20, 21] and [18] as follows: $$\begin{array}{}\text{(19)}& -\epsilon \frac{{\Psi }_{i+1}^{j+1}-2{\Psi }_i^{j+1}+{\Psi }_{i-1}^{j+1}}{{\mu }_{i,j+1}^2}+B_i{\Psi }_i^{j+1}=\epsilon \frac{{\Psi }_{i+1}^j-2{\Psi }_i^j+{\Psi }_{i-1}^j}{{\mu }_{i,j}^2}-C_i{\Psi }_i^j+f_i^{j+1}+f_i^j\end{array}$$with the discrete initial and boundary conditions $$\begin{array}{}\text{(20)}& \begin{array}{c}{\Psi }_i^0=0,i=0,1,2,\cdots ,K\hfill \\ {\Psi }_0^{j+1}={\Psi }_K^{j+1}=0,j=0,1,2,\cdots ,K-1\hfill \end{array}\end{array}$$where $B_i=2/k+a_i^{j+1}$, $C_i=a_i^j-2/k$, and the denominator function ${\mu }_{i,j+1}^2$ is to be determined from the following steps.

Consider the homogeneous part of Equation (19) with constant coefficients, where the constant coefficients are the lower bounds of the coefficient functions as follows: $$\begin{array}{}\text{(21)}& -\epsilon \frac{d}{d{x}^2}{\Psi }^{j+1}+B_i^{\ast }{\Psi }^{j+1}=0\end{array}$$where $B_i{x}_i,t_{j+1}\ge B_i^{\ast }>0$ and the corresponding characteristic equation is given by $-\epsilon {\theta }^2+B_i^{\ast }=0$. So, the two linearly independent solutions of Equation (21) are given by $$\begin{array}{}\text{(22)}& \exp {\theta }_{1}x\text{and}\exp {\theta }_{1}x,{\theta }_1,{\theta }_2=\pm \sqrt{\frac{B_i^{\ast }}{\epsilon }}\end{array}$$

Hence, ${\Psi }_i^{j+1}=c_1\exp {\theta }_1{x}_i+c_2\exp {\theta }_2{x}_i$ is the solution of Equation (21). Now, by taking the successive points and applying the theory of difference equations as in [15, 22], the following set of equations can be obtained: $$\begin{array}{}\text{(23)}& \begin{array}{ccc}{\Psi }_{i-1}^{j+1}& \exp {\theta }_1{x}_{i-1}& \exp {\theta }_2{x}_{i-1}\\ {\Psi }_i^{j+1}& \exp {\theta }_1{x}_i& \exp {\theta }_2{x}_i\\ {\Psi }_{i+1}^{j+1}& \exp {\theta }_1{x}_{i+1}& \exp {\theta }_2{x}_{i+1}\end{array}=0\end{array}$$

Rearranging Equation (23), we get $$\begin{array}{}\text{(24)}& \begin{array}{ccc}{\Psi }_{i-1}^{j+1}& \exp {\theta }_1{x}_i-h& \exp {\theta }_2{x}_i-h\\ {\Psi }_i^{j+1}& \exp {\theta }_1{x}_i& \exp {\theta }_2{x}_i\\ {\Psi }_{i+1}^{j+1}& \exp {\theta }_1{x}_i+h& \exp {\theta }_2{x}_i+h\end{array}=0\end{array}$$

Factorizing the above equation results in $$\begin{array}{}\text{(25)}& \exp {\theta }_1{x}_i+{\theta }_2{x}_i\begin{array}{ccc}{\Psi }_{i-1}^{j+1}& \exp -{\theta }_{1}h& \exp -{\theta }_{2}h\\ {\Psi }_i^{j+1}& 1& 1\\ {\Psi }_{i+1}^{j+1}& \exp {\theta }_{1}h& \exp {\theta }_{2}h\end{array}=0\end{array}$$

Solving the determinant of Equation (25) gives $$\begin{array}{}\text{(26)}& \exp {\theta }_1+{\theta }_2{x}_i\exp {\theta }_{2}h-\exp {\theta }_{1}h{\Psi }_{i-1}^{j+1}-\exp {\theta }_2-{\theta }_{1}h-\exp {\theta }_1-{\theta }_{2}h{\Psi }_i^{j+1}+\exp -{\theta }_{1}h-\exp -{\theta }_{2}h{\Psi }_{i+1}^{j+1}\end{array}$$

Replacing ${\theta }_1=-\sqrt{B_i^{\ast }/\epsilon }$ and ${\theta }_2=\sqrt{B_i^{\ast }/\epsilon }$ in Equation (26) and simplifying the terms using the double-angle formula for hyperbolic functions, we get the exact finite difference scheme for Equation (25) as follows: $$\begin{array}{}\text{(27)}& {\Psi }_{i-1}^{j+1}-2\cosh \sqrt{\frac{B_i^{\ast }}{\epsilon }}{\Psi }_i^{j+1}+{\Psi }_{i+1}^{j+1}=0\end{array}$$

After some rearrangements, we get $$\begin{array}{}\text{(28)}& -\epsilon \frac{{\Psi }_{i-1}^{j+1}-2{\Psi }_i^{j+1}+{\Psi }_{i+1}^{j+1}}{4\epsilon /B_i^{\ast }{\sinh }^{2}1/2h\sqrt{B_i^{\ast }/\epsilon }}+B_i^{\ast }{\Psi }_{i+1}^{j+1}=0\end{array}$$for $i=1,2,\cdots ,K-1$. As a result, the denominator function for the second-order derivative approximation can be found as $$\begin{array}{}\text{(29)}& {\mu }^2\epsilon ,h=4\frac{\epsilon }{B_i^{\ast }}{\sinh }^2\frac{1}{2}h\sqrt{\frac{B_i^{\ast }}{\epsilon }}\end{array}$$

For $i=1,2,3,\cdots ,K-1,\text{and}i=0,1,2,\cdots ,J-1$, we defined the general variable coefficient form of the denominator function as follows: $$\begin{array}{}\text{(30)}& {\mu }_{i,j+1}^2=4\frac{\epsilon }{B_i{x}_i,t_{j+1}}{\sinh }^2\frac{1}{2}h\sqrt{\frac{B_i{x}_i,t_{j+1}}{\epsilon }}\end{array}$$

Similarly considering the homogeneous differential equation of Equation (7) of the $j^{\text{th}}$ time level can be written as $$\begin{array}{}\text{(31)}& \epsilon \frac{d}{d{x}^2}{\Psi }^j-C_i{\Psi }^j=0\end{array}$$

The linearly independent solutions of Equation (31) are defined as $\exp {\theta }_{1}x$ and $\exp {\theta }_{2}x$, where ${\theta }_1,{\theta }_2=\pm \sqrt{C_i/\epsilon }$.

Following the same procedure, the $j^{\text{th}}$-level denominator function can be obtained as follows: $$\begin{array}{}\text{(32)}& {\mu }_{i,j}^2=4\frac{\epsilon }{C_i}{\sinh }^2\frac{1}{2}h\sqrt{\frac{C_i}{\epsilon }}\end{array}$$

3.3. The Fully Discrete Problem

In this section, we combined the temporal and spatial semidiscretized results by substituting Equation (30) and Equation (32) into Equation (7) as follows: $$\begin{array}{}\text{(33)}& -\epsilon \frac{{\Psi }_{i+1}^{j+1}-2{\Psi }_i^{j+1}+{\Psi }_{i-1}^{j+1}}{{\mu }_{i,j+1}^2}+B_i^{\ast }{\Psi }_i^{j+1}=\epsilon \frac{{\Psi }_{i+1}^j-2{\Psi }_i^j+{\Psi }_{i-1}^j}{{\mu }_{i,j}^2}-C_i^{\ast }{\Psi }_i^j+H_i^j\end{array}$$

Simplifying these equations, we get the following tridiagonal systems of equations for the proposed problem in (1) with its initial-boundary conditions: $$\begin{array}{}\text{(34)}& \begin{array}{c}{\mathcal{L}}_{\epsilon }^{K,J}:\equiv R_i^{-}{\Psi }_{i-1}^{j+1}+R_i^c{\Psi }_i^{j+1}+R_i^{+}{\Psi }_{i+1}^{j+1}=r_i^{-}{\Psi }_{i-1}^j+r_i^c{\Psi }_i^j+r_i^{+}{\Psi }_{i+1}^j+H_i^j\hfill \\ {\Psi }_i^0=S_0^0{x}_i,i=0,1,2,\cdots ,K-1\hfill \\ {\Psi }_0^{j+1}=S_l^{j+1}0,j=0,1,2,\cdots ,J-1\hfill \\ {\Psi }_K^{j+1}=S_r^{j+1}1,j=0,1,2,\cdots ,J-1\hfill \end{array}\end{array}$$where the coefficients of the tridiagonal system are determined for $i=1,2,3,\cdots ,K-1,\text{and}j=0,1,2,\cdots ,J-1$ as follows: $$\begin{array}{}\text{(35)}& \begin{array}{c}{R}_i^{-}=\frac{-\epsilon }{{\mu }_{i,j+1}^2}\hfill \\ R_i^c=\frac{2\epsilon }{{\mu }_{i,j+1}^2}+\frac{2}{k}+a_i^{j+1}\hfill \\ R_i^{+}=\frac{-\epsilon }{{\mu }_{i,j+1}^2}\hfill \\ H_i^j=f_i^{j+1}+f_i^j\hfill \end{array}\text{and}\begin{array}{c}{r}_i^{-}=\frac{\epsilon }{{\mu }_{i,j}^2}\hfill \\ r_i^c=\frac{-2\epsilon }{{\mu }_{i,j}^2}+\frac{2}{k}-a_i^j\hfill \\ r_i^{+}=\frac{\epsilon }{{\mu }_{i,j}^2}\hfill \end{array}\end{array}$$

Since the boundary values of ${\Psi }_0^{j+1}$ and ${\Psi }_K^{j+1}$ are given from Equation (1), we found the remaining values of ${\Psi }_i^{j+1},\text{for}i=1,2,3,\cdots ,K-1$, from the matrix below by means of the inverse matrix approach: $$\begin{array}{}\text{(36)}& \begin{array}{cccccc}{R}_0^c& R_0^{+}& 0& 0& \cdots & 0\\ R_1^{-}& R_1^c& R_1^{+}& 0& \cdots & 0\\ 0& R_2^{-}& R_2^c& R_2^{+}& 0& 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & \cdots & R_{N-2}^{-}& R_{K-2}^c& R_{N-2}^{+}\\ 0& \cdots & \cdots & \cdots & R_{K-1}^{-}& R_{K-1}^c\end{array}\begin{array}{c}{\Psi }_1^{j+1}\\ {\Psi }_2^{j+1}\\ {\Psi }_3^{j+1}\\ ⋮\\ {\Psi }_{N-2}^{j+1}\\ {\Psi }_{N-1}^{j+1}\end{array}=\begin{array}{cccccc}{r}_0^c& r_0^{+}& 0& 0& \cdots & 0\\ r_1^{-}& r_1^c& r_1^{+}& 0& \cdots & 0\\ 0& r_2^{-}& r_2^c& r_2^{+}& 0& 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & \cdots & r_{K-2}^{-}& r_{K-2}^c& r_{K-2}^{+}\\ 0& \cdots & \cdots & \cdots & r_{K-1}^{-}& r_{K-1}^c\end{array}\begin{array}{c}{\Psi }_1^j\\ {\Psi }_2^j\\ {\Psi }_3^j\\ ⋮\\ {\Psi }_{K-2}^j\\ {\Psi }_{K-1}^j\end{array}+\begin{array}{c}{H}_1^j\\ H_2^j\\ H_3^j\\ ⋮\\ H_{K-2}^j\\ H_{K-1}^j\end{array}\end{array}$$

Moreover, for the sake of simplicity, the matrix is rearranged to its compact form as follows: $$\begin{array}{}\text{(37)}& {\Psi }_i^{j+1}=R^{-1}{G}_i^j\end{array}$$where $G_i^j=r{\Psi }_i^j+H_i^j$.

Lemma 3.5 (discrete maximum principle).

Let ${\varphi }_i^{j+1}$ be any mesh function defined on ${\overline{\mathcal{D}}}_{xt}^{K,J}$ such that ${\varphi }_0^{j+1}\ge 0$ and ${\varphi }_K^{j+1}\ge 0$ [20]. Then, ${\mathcal{L}}^{K,J}{\psi }_i^{J+1}\ge 0$, for $\forall i=1,2,3,\cdots ,K-1$ implies that ${\varphi }_i^{j+1}\ge 0,\forall i=1,2,3,\cdots ,K$.

Proof 3.4.

Suppose that there exists a mesh point $i^{\ast },j+1$ for $i^{\ast }\in 1,2,3,\cdots ,K-1$ such that ${\varphi }_{i^{\ast }}^{j+1}=\text{mi}{\text{n}}_{0\le i\le K}{\varphi }_i^{j+1}$ and assume that ${\varphi }_{i^{\ast }}^{j+1}<0,\varphi i^{\ast }+1^{j+1}-{\varphi }_{i^{\ast }}^{j+1}>0,{\varphi }_{i^{\ast }-1}^{j+1}-{\varphi }_{i^{\ast }}^{j+1}>0$ and ${\mathcal{L}}^{K,J}{\varphi }_i^{J+1}\ge 0$. Since $i^{\ast }\notin 0,K$, then for $i^{\ast }\in 1,2,3,\cdots ,K-1$, we have that $$\begin{array}{}& {\mathcal{L}}^{K,J}{\psi }_{i^{\ast }}^{J+1}=-\epsilon \frac{{\varphi }_{i^{\ast }+1}^{j+1}-2{\varphi }_{i^{\ast }}^{j+1}+{\varphi }_{i^{\ast }-1}^{j+1}}{{\mu }_{i^{\ast }}^2}+B^{\ast }{\varphi }_{i^{\ast }}^{j+1}=-\epsilon \frac{{\varphi }_{i^{\ast }+1}^{j+1}-{\varphi }_{i^{\ast }}^{j+1}+{\varphi }_{i^{\ast }-1}^{j+1}-{\varphi }_{i^{\ast }}^{j+1}}{{\mu }_{i^{\ast }}^2}+B^{\ast }{\varphi }_{i^{\ast }}^{j+1}\le 0\end{array}$$

This contradicts to the assumption ${\mathcal{L}}^{K,J}{\psi }_i^{K+1}\ge 0$.

Hence, ${\varphi }_i^{J+1}\ge 0,\forall i=1,2,3,\cdots ,K$ as required.

Lemma 3.6 (uniform stability estimate).

Suppose that ${\varphi }_i^{j+1}$ be any mesh function at the ${j+1}^{\text{th}}$ level such that ${\varphi }_0^{j+1}=0={\varphi }_K^{j+1}$ [15]. Then, $$\begin{array}{}& {\varphi }_i^{j+1}\le \frac{1}{\beta }\underset{1\le j\le K-1}{\max }{\mathcal{L}}_{\epsilon }^{K,J}{\varphi }_i^{j+1},\text{for}1\le i\le K-1\end{array}$$

Proof 3.5.

Let $W=1/\beta {\max }_{1\le i\le K-1}{\mathcal{L}}_{\epsilon }^{K,J}{\varphi }_i^{j+1}$ and ${{\Omega }^{\pm }}_i^{j+1}$ be the barrier functions defined by $$\begin{array}{}& {{\Omega }^{\pm }}_i^{j+1}=W\pm {\varphi }^{j+1}\end{array}$$

Since ${{\Omega }^{\pm }}_0^j={{\Omega }^{\pm }}_K^j=W>0$ and $1\le i\le K-1$, it gives $$\begin{array}{}& {\mathcal{L}}_{\epsilon }^{K,J}{{\Omega }^{\pm }}_i^{j+1}=-\epsilon \frac{W\pm {\varphi }_{i+1}^{j+1}-2W\pm {\varphi }_i^{j+1}+W\pm {\varphi }_{i-1}^{j+1}}{{\mu }_i^2}+\frac{1}{k}+a_i^{j+1}W\pm {\varphi }_i^{j+1}=\frac{1}{k}+a_i^{j+1}W\pm {\mathcal{L}}_{\epsilon }^{K,J}{\varphi }_i^{j+1}=\frac{1}{\beta }\frac{1}{k}+a_i^{j+1}\max {\mathcal{L}}_{\epsilon }^{K,J}{\varphi }_i^{j+1}\pm {\mathcal{L}}_{\epsilon }^{K,J}{\varphi }_i^{j+1}\end{array}$$

Since $1/k+a_i^{j+1}>\beta$, we have ${\mathcal{L}}_{\epsilon }^{K,J}{{\Omega }^{\pm }}_i^{j+1}\ge 0$. Thus, by discrete maximum principle in Lemma (3.5), we obtain ${{\Omega }^{\pm }}_i^{j+1}\ge 0,0\le i\le K.$

Hence, the uniform stability estimate of the discrete scheme at the ${j+1}^{\text{th}}$ time level is proven.

Lemma 3.7.

For a fixed mesh and for all integers $i$, we have that $$\begin{array}{}& \underset{\epsilon ⟶0}{\lim }\underset{1\le i\le K-1}{\max }\frac{\exp -C{x}_i/\sqrt{\epsilon }}{{\epsilon }^{i/2}}=0\text{and}\underset{\epsilon ⟶0}{\lim }\underset{1\le i\le K-1}{\max }\frac{\exp -C1-x_{i}C1-x_i/\sqrt{\epsilon }}{{\epsilon }^{i/2}}=0\end{array}$$

For the detailed proof of this lemma, interested readers can refer to Lemma 5.3 in [16].

4. Convergence Analysis

Let ${\Psi }_i$ and $\psi x_i$ be the approximate and exact solutions of Problem (1); then, the local truncation error of the proposed NSFDM from Equation (33) is defined as follows: $$\begin{array}{}\text{(38)}& {\mathcal{L}}_{\epsilon }^N{\Psi }_i-{\psi }_i={\mathcal{L}}_{x,\epsilon }-{\mathcal{L}}_{x,\epsilon }^K{\psi }_i\end{array}$$

Using the NSFDM results in $$\begin{array}{}\text{(39)}& {\mathcal{L}}_{\epsilon }^{K,J}{\Psi }_i-\psi x_i=-\epsilon {\psi }_i^{\prime \prime }+\epsilon \frac{{\psi }_{i+1}-2{\psi }_i+{\psi }_{i-1}}{{\mu }_{i,j}^2}\end{array}$$

By using the Taylor series expansions of ${\psi }_{i-1}$ and ${\psi }_{i+1}$ from Equation (39), we get $$\begin{array}{}\text{(40)}& {\mathcal{L}}_{\epsilon }^{K,J}{\Psi }_i-{\psi }_i=-\epsilon {\psi }_i^{\prime \prime }+\frac{\epsilon }{{\mu }_{i,j}^2}{h}^2{\psi }^{\prime \prime }{\zeta }_i+\frac{h^4}{12}{\psi }^4{\zeta }_i{\zeta }_i\in x_{i-1},x_{i+1}\end{array}$$

Furthermore, the truncated Taylor series expansion of $1/{\mu }_{i.j}^2$ is given by $$\begin{array}{}\text{(41)}& \frac{1}{{\mu }_{i,j}^2}=\frac{1}{h^2}-\frac{B_i}{12\epsilon }+\frac{B_i{h}^2}{240{\epsilon }^2},\text{where},B_i=2/k+a_i^j\end{array}$$

Substituting Equation (41) into Equation (40) and using Lemma (2.3) and (3.7) give $$\begin{array}{}& {\mathcal{L}}_{\epsilon }^{K,J}{\Psi }_i-{\psi }_i\le -\frac{B_i}{12}{h}^2+\frac{B_i^2}{240}{h}^4+O{h}^8\le C{h}^2\end{array}$$

Since $h$ is assumed to be very small, we have the relations $h^2>h^4>h^8>\cdots$.

Hence, the absolute error is defined by $$\begin{array}{}\text{(42)}& {\Psi }_i-\psi x_i\le C{h}^2\end{array}$$where $C$ is independent of the perturbation parameter $\epsilon$ and the step size $h$. Based on the above results, we stated the following theorems.

Theorem 4.1.

Let ${\Psi }_i^j$ be the numerical solution of Equation (34) and $\psi x,t_j$ be the solution of the problem in Equation (6) both at time level $j$; then, $$\begin{array}{}\text{(43)}& \underset{0\le i\le K,0\le j\le J}{\max }{\Psi }_i^j-\psi x,t_j\le C{h}^2\end{array}$$

Therefore, the NSFDM is a second-order uniformly convergent scheme in space variable. The subsequent theorem is the combined results of Lemma 3.2 and Theorem 4.1.

Theorem 4.2.

Let ${\Psi }_i^j$ be the $j^{\text{th}}$ time-level approximate solution of Equation (34) and $\psi x,t$ be the solution of Equation (1). Then, $$\begin{array}{}\text{(44)}& \underset{0\le i\le K,0\le j\le J}{\max }{\Psi }_i^j-\psi x,t_j\le C{h}^2+k^2\end{array}$$

From the above theorem, it is concluded that the proposed nonstandard finite difference and Crank–Nicolson schemes are second-order uniformly convergent methods.

5. Numerical Illustrations and Discussions

The absence of exact solutions is a major issue in numerical computations and can be caused on by a number of characteristics, including larger dimension, nonlinearity, and other model problem parameters. As a result, this condition poses a challenge to the investigation of the error estimate of the presented model problems generally. To get around these problems, we estimate the maximum point-wise errors (MPWEs) and corresponding order of convergences (OCs) of the approximate solutions through the doubling mesh principle. The proposed method was confirmed through simulation and discussion using three model examples.

Example 5.1.

Consider a SP unsteady parabolic-type RDP of the form in [14] $$\begin{array}{}& \begin{array}{c}{\Psi }_t-\epsilon {\Psi }_{xx}+\frac{1+x^2}{2}\Psi =t^3,x,t\in 0,1\times 0,1\hfill \\ \Psi x,0=0,x\in 0,1\hfill \\ D_l\Psi 0,t=-\frac{128}{35}{\pi }^{-1/2}{t}^{7/2},\text{and}{D}_r\Psi 1,t=-\frac{128}{35}{\pi }^{-1/2}{t}^{7/2},t\in 0,1\hfill \end{array}\end{array}$$

Example 5.2.

Consider a SP unsteady parabolic-type RDP of the form in [15] $$\begin{array}{}& \begin{array}{c}{\Psi }_t-\epsilon {\Psi }_{xx}+\frac{1+x^2}{2}\Psi ={1+t}^3,x,t\in 0,1\times 0,1\hfill \\ \Psi x,0=0,x\in 0,1\hfill \\ \Psi 0,,t=\Psi 1,t=0,t\in 0,1\hfill \end{array}\end{array}$$

Example 5.3.

Consider a SP unsteady parabolic-type RDP of the form in [23] $$\begin{array}{}& \begin{array}{c}{\Psi }_t-\epsilon {\Psi }_{xx}+1+x^2+t^2{e}^t\Psi =e^t-1+\sin \pi x,x,t\in 0,1\times 0,1\hfill \\ \Psi x,0=0,x\in 0,1\hfill \\ \Psi 0,,t=U1,t=0,t\in 0,1\hfill \end{array}\end{array}$$

Since the exact solution of the model problem is unknown, we let ${\Psi }^{K,J}{x}_i,t_j$ and ${\Psi }^{2K,2J}{x}_i,t_j$ be the approximate solutions with mesh lengths $h,k\ne 0$ and $h/2,k/2\ne 0$, respectively; then, we used the following formulas:

• The MPWE is given by

$$\begin{array}{}\text{(45)}& E^{K,J}=\underset{0\le i,j\le K,J}{\max }{\Psi }^{K,J}{x}_i,t_j-{\Psi }^{2K,2J}{x}_{2i},t_{2j}\end{array}$$

• The $\epsilon -$uniform MPWE is given by

$$\begin{array}{}\text{(46)}& E_{\epsilon }^{K,J}=\underset{0\le i,j\le K,J}{\max }{E}^{K,J}\end{array}$$

• The $\epsilon -$uniform rate of convergence is given by

$$\begin{array}{}\text{(47)}& R_{\epsilon }^{K,J}={\log }_2\frac{E_{\epsilon }^{K,J}}{E_{\epsilon }^{2K,2J}}\end{array}$$

6. Discussions

To validate the legitimacy of the proposed method, three model examples were considered. The MPWEs and the corresponding OCs are presented in Tables 1, 2, and 3. For a fixed value of the perturbation parameter, $\epsilon$, and decreasing the values of the step sizes $h$ and $k$ across a row, it is seen that the MPWEs decrease. This is to mean that the accuracy of the method increases. On the other hand, fixing the values of step sizes $h$ and $k$ and a decrease in the values of the perturbation parameter $\epsilon$ down a column, we observed that the values of the MPWEs are constant in a uniform pattern in each column and the OC approaching to two. These results confirmed that the proposed method is a second-order uniformly convergent scheme. In addition to this, the 3D and 2D plots depicted in Example 5.1 (Figures 1 and 2), Example 5.2 (Figures 3 and 4), and Example 5.3 (Figures 5 and 6) elucidated that the decrease in the values of the perturbation parameter $\epsilon$ leads the proposed problem to have twin boundary layers near $x=0$ and $x=1$ as discoursed in the analytical aspects of the proposed problem. The comparison Tables 4 and 5 exhibit that the present method is well performed as compared with the methods in [14, 15]. The overall discussions presented in this work confirmed that the theoretical analysis coincides with the results of the practical examples.

Table 1

MPWEs $E_{\epsilon }^{K,J}$ and order of convergence $R_{\epsilon }^{K,J}$ for Example 5.1.

$\mathit{\epsilon }\downarrow \mathbf{K},\mathbf{J}⟶$	32, 8	64, 16	128, 32	256, 64	512, 128
$10^{-6}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$10^{-7}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$10^{-8}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$10^{-10}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$10^{-50}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$E_{\epsilon }^{K,J}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
$R_{\epsilon }^{K,J}$	1.9988	1.9997	1.9999	2.0000

Table 2

MPWEs $E_{\epsilon }^{K,J}$ and order of convergence $R_{\epsilon }^{K,J}$ for Example 5.2.

$\mathit{\epsilon }\downarrow \mathbf{K},\mathbf{J}⟶$	32, 10	64, 40	128, 160	256, 640
$2^{-18}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0797\text{e}-06$
	1.9995	1.9999	2.0000
$2^{-20}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0797\text{e}-06$
	1.9995	1.9999	2.0000
$2^{-22}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0798\text{e}-05$
	1.9995	1.9999	2.0000
$2^{-24}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0798\text{e}-05$
	1.9995	1.9999	2.0000
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$2^{-50}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0798\text{e}-05$
	1.9995	1.9999	2.0000
$E_{\epsilon }^{K,J}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0798\text{e}-05$
$R_{\epsilon }^{K,J}$	1.9995	1.9999	2.0000

Table 3

MPWEs $E_{\epsilon }^{K,J}$ and order of convergence $R_{\epsilon }^{K,J}$ for Example 5.3.

$\mathit{\epsilon }\downarrow \mathbf{K},\mathbf{J}⟶$	32, 10	64, 40	128, 160	256, 640	512, 2560
$2^{-16}$	$9.4881e-04$	$2.3940e-04$	$5.9878e-05$	$3.7372e-06$	$1.2584-06$
	1.9867	1.5426	−1.9513	2.0000
$2^{-18}$	$9.4882e-04$	$2.3957e-04$	$5.9887e-05$	$3.7465e-06$	$1.2673-06$
	1.9857	1.9986	2.0000	2.0000
$2^{-20}$	$9.4882e-04$	$2.3957e-04$	$5.9887e-05$	$3.7465e-06$	$1.2673-06$
	1.9857	1.9986	2.0000	2.0000
$2^{-22}$	$9.4882e-04$	$2.3957e-04$	$5.9887e-05$	$3.7465e-06$	$1.2673-06$
	1.9857	1.9986	2.0000	2.0000
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$2^{-50}$	$9.4882e-04$	$2.3957e-04$	$5.9887e-05$	$3.7465e-06$	$1.2673-06$
	1.9857	1.9986	2.0000	2.0000
$E_{\epsilon }^{K,J}$	$9.4882e-04$	$2.3957e-04$	$5.9887e-05$	$3.7465e-06$	$1.2673-06$
$R_{\epsilon }^{K,J}$	1.9857	1.9986	2.0000	2.0000

[figure(s) omitted; refer to PDF]

Table 4

Comparison of the MPWEs and corresponding OCs for Example 5.1.

$\mathit{\epsilon }\downarrow \mathbf{K},\mathbf{J}⟶$	32, 8	64, 16	128, 32	256, 64	512, 128
Our method
$10^{-8}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
	1.9988	1.9997	1.9999	2.0000
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$E_{\epsilon }^{K,J}$	$2.1119\text{e}-03$	$5.2840\text{e}-04$	$1.3213\text{e}-04$	$3.3033\text{e}-05$	$8.2584\text{e}-06$
$R_{\epsilon }^{K,J}$	1.9988	1.9997	1.9999	2.0000
Results in [14]
$10^{-8}$	$2.5497\text{e}-02$	$1.2508\text{e}-02$	$6.1896\text{e}-03$	$3.0782\text{e}-03$	$1.5348\text{e}-03$
	1.0275	1.0149	1.0078	1.0040
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$E_{\epsilon }^{K,J}$	$2.5497\text{e}-02$	$1.2508\text{e}-02$	$6.1896\text{e}-03$	$3.0782\text{e}-03$	$1.5348\text{e}-03$
$R_{\epsilon }^{K,J}$	1.0275	1.0149	1.0078	1.0040

Table 5

Comparison of MPWEs and corresponding OCs for Example 5.2.

$\mathit{\epsilon }\downarrow \mathbf{K},\mathbf{J}⟶$	32, 10	64, 40	128, 160	256, 640
Our method
$2^{-20}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0797\text{e}-06$
	1.9995	1.9999	2.0000
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$E_{\epsilon }^{N,M}$	$3.1521\text{e}-03$	$7.8832\text{e}-04$	$4.9275\text{e}-05$	$3.0798\text{e}-05$
$R_{\epsilon }^{N,M}$	1.9995	1.9999	2.0000
Results in [15]
$2^{-20}$	$1.6700\text{e}-01$	$4.1800\text{e}-02$	$1.05\text{e}-02$	$0.2600\text{e}-02$
	1.9983	1.9931	2.0138
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$
$E_{\epsilon }^{K,J}$	$1.6700\text{e}-01$	$4.1800\text{e}-02$	$1.05\text{e}-02$	$0.2600\text{e}-02$
$R_{\epsilon }^{K,J}$	1.9983	1.9931	2.0138

7. Conclusion

In the present work, we formulated and analyzed a uniformly convergent scheme for the numerical treatment of SP unsteady RDPs. The continuous problem domain is semidiscretized using the Crank–Nicolson method in the temporal direction and NSFDM in the spatial direction on uniform meshes. The various convergence analyses prove that the proposed scheme is uniformly convergent with second-order accuracy in the temporal and spatial directions. From the numerical illustrations and discussions, it has been observed that the theoretical findings and the results of the practical model problems are in agreement and the formulated method well performs compared with some methods in the literature. As further work of the research, it is recommended to apply the method for solving nonlinear singularly perturbed PDEs.

Funding

This research work does not have any financial support.

References

[1] H. G. Roos, Robust Numerical Methods for Singularly Perturbed Differential Equations, 2008.

[2] D. Boffi, "Convection-diffusion problems. An introduction to their analysis and numerical solution," SIAM Review, vol. 62, pp. 987-988, 2020.

[3] V. Sigamani, J. J. Miller, R. Narasimhan, P. Mathiazhagan, F. Victor, Differential Equations and Numerical Analysis,DOI: 10.1007/978-81-322-3598-9, 2016.

[4] J. J. Miller, E. O’riordan, G. I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions,DOI: 10.1142/2933, 1996.

[5] T. Linß, "Sufficient conditions for uniform convergence on layer-adapted grids," Applied Numerical Mathematics, vol. 37 no. 1-2, pp. 241-255, DOI: 10.1016/S0168-9274(00)00043-X, 2001.

[6] S. C. S. Rao, S. Kumar, J. Singh, "A discrete Schwarz waveform relaxation method of higher order for singularly perturbed parabolic reaction–diffusion problems," Journal of Mathematical Chemistry, vol. 58 no. 3, pp. 574-594, DOI: 10.1007/s10910-019-01086-1, 2020.

[7] S. Kumar, J. S. Aakansha, H. Ramos, "Parameter-uniform convergence analysis of a domain decomposition method for singularly perturbed parabolic problems with Robin boundary conditions," Journal of Applied Mathematics and Computing, vol. 69 no. 2, pp. 2239-2261, DOI: 10.1007/s12190-022-01832-w, 2023.

[8] Y. C. Ji, F. H. Fenton, "Numerical solutions of reaction-diffusion equations: application to neural and cardiac models," American Journal of Physics, vol. 84 no. 8, pp. 626-638, DOI: 10.1119/1.4953167, 2016.

[9] S. Kumar, M. Kumar, "Analysis of some numerical methods on layer adapted meshes for singularly perturbed quasi-linear systems," Numerical Algorithms, vol. 71 no. 1, pp. 139-150, DOI: 10.1007/s11075-015-9989-2, 2016.

[10] M. K. Kadalbajoo, A. Awasthi, "Crank–Nicolson finite difference method based on a midpoint upwind scheme on a non-uniform mesh for time-dependent singularly perturbed convection–diffusion equations," International Journal of Computer Mathematics, vol. 85 no. 5, pp. 771-790, DOI: 10.1080/00207160701459672, 2008.

[11] S. C. S. Rao, S. Kumar, "Second order global uniformly convergent numerical method for a coupled system of singularly perturbed initial value problems," Applied Mathematics and Computation, vol. 219 no. 8, pp. 3740-3753, DOI: 10.1016/j.amc.2012.09.075, 2012.

[12] S. Kumar, J. Vigo-Aguiar, "A parameter-uniform grid equidistribution method for singularly perturbed degenerate parabolic convection–diffusion problems," Journal of Computational and Applied Mathematics, vol. 404,DOI: 10.1016/j.cam.2020.113273, 2022.

[13] C. Clavero, J. L. Gracia, "On the uniform convergence of a finite difference scheme for time dependent singularly perturbed reaction-diffusion problems," Applied Mathematics and Computation, vol. 216 no. 5, pp. 1478-1488, DOI: 10.1016/j.amc.2010.02.050, 2010.

[14] S. Kumar, H. Ramos, "Parameter-uniform approximation on equidistributed meshes for singularly perturbed parabolic reaction-diffusion problems with Robin boundary conditions," Applied Mathematics and Computation, vol. 392,DOI: 10.1016/j.amc.2020.125677, 2021.

[15] J. B. Munyakazi, K. C. Patidar, "A fitted numerical method for singularly perturbed parabolic reaction-diffusion problems," Computational and Applied Mathematics, vol. 32 no. 3, pp. 509-519, DOI: 10.1007/s40314-013-0033-7, 2013.

[16] A. Andargie Tiruneh, G. Adamu Derese, M. Amsalu Ayele, "Singularly perturbed reaction diffusion problem with large spatial delay via non-standard fitted operator method," Research in Mathematics, vol. 10 no. 1,DOI: 10.1080/27684830.2023.2171698, 2023.

[17] D. M. Tefera, A. A. Tiruneh, G. A. Derese, "Numerical treatment on parabolic singularly perturbed differential difference equation via fitted operator scheme," Abstract and Applied Analysis, vol. 2021 no. 1,DOI: 10.1155/2021/1661661, 2021.

[18] T. B. Mekonnen, G. F. Duressa, "Computational method for singularly perturbed two-parameter parabolic convection-diffusion problems," Cogent Mathematics & Statistics, vol. 7 no. 1,DOI: 10.1080/25742558.2020.1829277, 2020.

[19] M. Amsalu Ayele, A. Andargie Tiruneh, G. Adamu Derese, "Fitted scheme for singularly perturbed time delay reaction-diffusion problems," Iranian Journal of Numerical Analysis and Optimization, vol. 13 no. 2, pp. 262-284, DOI: 10.22067/ijnao.2022.77453.1161, 2023.

[20] R. E. Mickens, "Nonstandard finite difference schemes for differential equations," Journal of Difference Equations and Applications, vol. 8 no. 9, pp. 823-847, DOI: 10.1080/1023619021000000807, 2002.

[21] T. B. Mekonnen, G. F. Duressa, "Uniformly convergent numerical method for two-parametric singularly perturbed parabolic convection-diffusion problems," Journal of Applied and Computational Mechanics, vol. 7 no. 1, pp. 535-545, DOI: 10.22055/jacm.2020.35193.2596, 2020.

[22] M. M. Woldaregay, G. F. Duressa, "Boundary layer resolving exact difference scheme for solving singularly perturbed convection-diffusion-reaction equation," Mathematical Problems in Engineering, vol. 2022,DOI: 10.1155/2022/2043323, 2022.

[23] F. W. Gelu, G. F. Duressa, G. D. Kebede, "A robust higher-order fitted mesh numerical method for solving singularly perturbed parabolic reaction diffusion problems," Results in Applied Mathematics, vol. 20,DOI: 10.1016/j.rinam.2023.100405, 2023.