Full article: Crank–Nicolson method for solving time-fractional singularly perturbed delay partial differential equations

Formulae display: $MathJax Logo$ ?Mathematical formulae have been encoded as MathML and are displayed in this HTML version using MathJax in order to improve their display. Uncheck the box to turn MathJax off. This feature requires Javascript. Click on a formula to zoom.

ABSTRACT

In this paper, we propose a uniformly convergent numerical scheme for a class of time-fractional singularly perturbed delay partial differential equations exhibiting a right regular boundary layer. The time-fractional derivative is considered in the Caputo sense with order $β \in (0, 1)$ . The domain is discretized with a uniform mesh in both time and space directions. The numerical scheme comprises the discretization technique given by the Crank–Nicolson method in the temporal direction and the non-standard finite difference method in the spatial direction. To show the parameter uniform convergence of the proposed method, the boundedness of truncation error and stability analysis are performed. It is shown that the proposed scheme is second-order convergent in the temporal direction and first-order convergent in the spatial direction. Numerical examples are presented, and their numerical simulations confirm the theoretical findings.

KEYWORDS:

1. Introduction

Delay partial differential equations (DPDEs) give a more realistic model of physical problems by considering both the past history and the present state of the system (Ansari et al., Citation2007). Often, the delay accounts for the time it takes for the process to occur. For example, gestation times, incubation periods, and transport delays. DPDEs have many applications in science and engineering (Kumar & Kumari, Citation2019). When the highest spatial derivative of the DPDE is multiplied by a small parameter $(ϵ)$ , it is said to be singularly perturbed delay partial differential equation (SPDPDE), and ϵ is said to be a perturbation parameter. As the perturbation parameter ϵ → 0, the solution exhibits a boundary or an interior layer. Due to the existence of a boundary or an interior layer on the solution of SPDPDE, when we apply standard numerical methods like the finite difference method, the finite element method, and the finite volume method to a uniform mesh, large oscillations occur in the computed solution, leading to numerical instability and inaccurate solutions (Gowrisankar & Natesan, Citation2017). To overcome this difficulty, ϵ - uniformly convergent numerical methods, that is, numerical methods whose convergence property, accuracy, and computational cost do not depend on the perturbation parameter and only depend on the number of mesh points, must be developed (Miller et al., Citation1996).

In general, there are two kinds of ϵ - uniform numerical methods for solving singularly perturbed problems in the context of finite difference methods. Namely, fitted operator methods (FOMs) and fitted mesh methods (FMMs). In FOMs, exponential fitting factor or positive denominator function can be used to control the rapid change of the solution in the boundary or interior layer region(s). When we use exponential fitting factor and denominator function on a term containing perturbation parameter, the methods are said to be the exponentially fitted operator method (EFOM) and the non-standard finite difference method (NSFDM), respectively. Whereas, FMMs employ non-uniform meshes that will be fine in the boundary layer region and coarse outside the boundary layer region(s).

In the past few decades, fractional calculus has gained more attention because of its wide applications in science and engineering (Sun et al., Citation2018). It is a generalization of classical calculus to arbitrary-order derivatives and integrals and can be used in the oil industry, neural networks, mass transport in fluids, control theory, signal processing, viscoelasticity, rheology, finance, dispersion of pollutants in the atmosphere, and so on (Kilbas et al., Citation2006). Due to the non-local properties of fractional derivative and integral operators, they are used to model memory effects and hereditary properties of physical problems more accurately than integer-order derivative and integral operators, which are neglected in the integer-order integral and derivative operators. Fractional partial differential equations (FPDEs) that are described by non-integer-order derivatives model many complex systems such as anomalous diffusion, viscoelastic materials, etc. (Dai et al., Citation2015). Due to the widespread application of FPDEs in engineering and science, this area has been receiving significant attention, and as a result, numerous analytical (Ata & Kıymaz, Citation2023; Atangana, Citation2014; Chen et al., Citation2007), semi-analytical (Daraghmeh et al., Citation2020; Odibat & Momani, Citation2006), and numerical methods (Karatay & Bayramoglu, Citation2012; Kumar et al., Citation2022; Lin & Xu, Citation2007) were introduced to solve FPDEs.

In this paper, we consider the following time-fractional singularly perturbed delay partial differential equation (TF-SPDPDE):

(1)

\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}

(1)

with the interval and boundary conditions:

(2)

\begin{array}{l} \begin{array}{l} u (x, t) & = φ_{b} (x, t), \forall (x, t) \in [0, 1] \times [- η, 0], \\ u (0, t) & = φ_{l} (t), ∀t \in (0, T), \\ u (1, t) & = φ_{r} (t), ∀t \in (0, T), \end{array} \end{array}

(2)

where $0 < ϵ < < 1$ denotes the perturbation parameter, η > 0 denotes the delay parameter $(η = T / K for some integer K > 1)$ . Here, the domain be $Ω = Ω_{x} \times Ω_{t} = (0, 1) \times (0, T),$ and $Γ = Γ_{l} \cup Γ_{b} \cup Γ_{r}$ is its boundary, where $Γ_{l}$ and $Γ_{r}$ are the left and right sides of $Ω$ corresponding to x = 0 and 1, respectively, and $Γ_{b} = [0, 1] \times [- η, 0]$ . Also, $\frac{\partial^{β} u (x, t)}{\partial t^{β}}$ is the Caputo time-fractional derivative of order $β \in (0, 1)$ , usually denoted by $_{0} D_{t}^{β} u (x, t)$ and defined by

_{0} D_{t}^{β} u (x, t) = \frac{1}{Γ (1 - β)} \int_{0}^{t} (t - s)^{- β} \frac{∂u (x, s)}{∂s} ds,

where Γ represents the gamma function defined as for $z \in C$ with $Re (z) > 0,$

Γ (z) = \int_{0}^{\in f t y} x^{z - 1} \exp (- x) dx .

If $0 < β < 1$ , the problem given in EquationEquation (1)(1) $\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}$ (1) is a TF-SPDPDE, which is the generalization of classical SPDPDE. This type of problem is widely used in oil reservoir simulations, transport of mass and energy, dispersion of chemicals in reactors, etc. (Choudhary et al. (Citation2022).

Suppose that $p (x) \geq ρ > 0, q (x, t) \geq α > 0, r (x, t) \break \geq γ > 0,$ and $g (x, t)$ are bounded smooth functions in the rectangular domain $\overset{―}{Ω} = [0, 1] \times [0, T]$ . In addition, let us assume that the initial boundary functions given are sufficiently smooth and satisfy the following compatibility condition at the corner points of $Ω$ (i.e. at $(0, 0)$ and $(1, 0)$ ). This means

\begin{array}{l} φ_{b} (0, 0) = φ_{l} (0), φ_{b} (1, 0) = φ_{r} (0), \end{array}

\begin{array}{l} _{0} D_{t}^{β} φ_{l} (0) - ϵ \frac{\partial^{2} φ_{b} (0, 0)}{\partial x^{2}} + p (0) \frac{\partial φ_{b} (0, 0)}{∂x} + q (0, 0) φ_{b} (0, 0) \\ = - r (0, 0) φ_{b} (0, - η) + g (0, 0), \end{array}

\begin{array}{l} _{0} D_{t}^{β} φ_{r} (0) - ϵ \frac{\partial^{2} φ_{b} (1, 0)}{\partial x^{2}} + p (1) \frac{\partial φ_{b} (1, 0)}{∂x} + q (1, 0) φ_{b} (1, 0) \\ = - r (1, 0) φ_{b} (1, - η) + g (1, 0) . \end{array}

Under these sufficiently smooth assumptions and compatibility conditions imposed on the initial and boundary data, the interval boundary value problem (1)-(2) admits a unique solution $u (x, t)$ that exhibits a regular boundary layer of width $O (ϵ)$ at $x = 1 of Ω$ , as the perturbation parameter becomes very small $(ϵ \to 0)$ . Due to the non-locality of the fractional derivative operator and the involvement of delay and perturbation parameters, it is not simple to solve this problem effectively and more accurately. Hence, designing the ϵ - uniform numerical methods for such type of problems is an important task and is the objective of the present work.

There have been several investigations into the cases of numerical treatment of integer-order singularly perturbed delay partial differential equations. (e.g. see (Ansari et al., Citation2007; Bashier & Patidar, Citation2011; Duressa et al., Citation2023; Govindarao & Mohapatra, Citation2018; Kumar & Kumari, Citation2020; Negero & Duressa, Citation2022; Tirunehi et al., Citation2023) and the references therein). Moreover, FPDEs that are not singularly perturbed have been widely studied up to now. For example, Rihan (Citation2010) presented fully implicit $θ$ - methods to solve linear and non-linear time-fractional delay partial differential equations. Qiu et al. (Citation2019) developed a numerical method comprising the L1 scheme and central finite difference approximation formula to solve one-dimensional time-fractional Burgers equation. The method was shown unconditional stable and convergent of order 2−α, α ∈ (0, 1) in the temporal direction and of order two in spatial direction. Choudhary et al. (Citation2022) proposed a second-order convergent numerical scheme for a class of time-fractional partial differential equations with delay in time. The authors used a numerical method comprising the Crank–Nicolson method for time discretization and Spline method with a tension factor for spatial discretization on a uniform mesh. However, the problems considered in the works above are not singularly perturbed. Relatively, if the problems are not singularly perturbed (in the absence of the perturbation parameter), then we may not face numerical difficulties to find their numerical solution

However, few studies are conducted on the numerical treatment of time-fractional singularly perturbed delay partial differential equations. For instance, Liu et al. (Citation2021) designed a stable numerical scheme to solve a one- and two-dimensional time-fractional singularly perturbed partial differential equation with variable coefficients. They employed the Tailored finite point method to approximate the space derivative and the Grunwald-Letnikov and L1 scheme to approximate the time-fractional derivative. The scholars (Kumar & Vigo-Aguiar, Citation2021) developed a uniformly convergent numerical scheme to obtain an approximate solution of singularly perturbed time-fractional partial differential equations with delay in time. They used the L1 scheme on a uniform mesh to approximate the time-fractional derivative and stable finite difference method with piecewise uniform meshes (which are adequate for FMMs) for approximating the spacial derivatives. Motivated by the study of (Kumar, Kamalesh and Vigo-Aguiar, J, 2021), the purpose of this paper is to design $\in$ - uniformly convergent numerical method to solve time-fractional singularly perturbed partial differential equations with a large time delay given in Equationequations (1)(1) $\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}$ (1) -(Equation2(2) $\begin{array}{l} \begin{array}{l} u (x, t) & = φ_{b} (x, t), \forall (x, t) \in [0, 1] \times [- η, 0], \\ u (0, t) & = φ_{l} (t), ∀t \in (0, T), \\ u (1, t) & = φ_{r} (t), ∀t \in (0, T), \end{array} \end{array}$ (2) ). We apply the Crank–Nicolson method and the NSFDM to approximate the time-fractional and space derivatives, respectively. Since NSFDM is one of a class of FOM, it has the advantage over FMM that it does not require any prior information about the location and width of the boundary layer(s).

1.1. Properties of the continuous problem

In this section, we present some properties of the continuous problem given in EquationEquation (1)(1) $\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}$ (1) -EquationEquation (2)(2) $\begin{array}{l} \begin{array}{l} u (x, t) & = φ_{b} (x, t), \forall (x, t) \in [0, 1] \times [- η, 0], \\ u (0, t) & = φ_{l} (t), ∀t \in (0, T), \\ u (1, t) & = φ_{r} (t), ∀t \in (0, T), \end{array} \end{array}$ (2) , which are important to analyze the semi-discrete and the full-discrete schemes of the continuous problem that we will discuss in the next sections.

Lemma 1.1.

The Continuous Maximum principle Let $Ψ (x, t) \in C^{2} (Ω) \cap C^{0} (\overset{―}{Ω})$ satisfying $Ψ (x, t) \geq 0, \forall (x, t) \in Γ$ and $L Ψ (x, t) \geq 0, \forall (x, t) \in Ω$ . Then $Ψ (x, t) \geq 0, \break∀ (x, t) \in \overset{―}{Ω}$ .

Proof.

Let us assume that there exists $(ν, ξ) \in \overset{―}{Ω}$ such that

\begin{array}{l} (ν, ξ) = min_{(x, t) \in \overset{―}{Ω}} Ψ (x, t) \end{array}

and $Ψ (ν, ξ) < 0$ . From the given hypothesis, we have $(ν, ξ) \notin Γ as Ψ \geq 0 on Γ$ and $(ν, ξ) \in Ω$ . As $Ψ$ attains minimum at $(ν, ξ),$ we have $\frac{∂Ψ (ν, ξ)}{∂x} = 0, \frac{\partial^{2} Ψ (ν, ξ)}{\partial x^{2}} \geq 0$ and $_{0} D_{t}^{β} Ψ (ν, ξ) \leq 0$ and since $q (ν, ξ) \geq 0$ , we obtain the following:

\begin{array}{l} L Ψ (ν, ξ) = & _{0} D_{t}^{β} Ψ (ν, ξ) - ϵ \frac{\partial^{2} Ψ (ν, ξ)}{\partial x^{2}} + p (ν) \frac{∂Ψ (ν, ξ)}{∂x} \\ + q (ν, ξ) Ψ (ν, ξ) < 0, \end{array}

which contradicts the assumption $L Ψ (x, t) \geq 0$ in $Ω$ . Hence, $Ψ (ν, ξ) \geq 0$ . Therefore, we conclude that $Ψ (x, t) \geq 0, \forall (x, t) \in \overset{―}{Ω}$ .

Lemma 1.2.

The solution of the continuous problem (1)-(2) satisfy the following estimate

| u (x, t) - φ_{b} (x, 0) | \leq Ct, (x, t) \in \overset{―}{Ω} .

Proof.

Using the compatibility conditions at the corner points (i.e. $(0, 0) and (1, 0)$ ) and applying Lemma 1.1 on the barrier functions defined by

Ψ^{\pm} (x, t) = {\begin{cases} Ct \pm (u (x, t) - φ_{b} (x, 0)), 0 \leq t \leq T, \\ Ct \pm (u (x, t) - φ_{b} (x, t)), - η \leq t \leq 0, \end{cases}

gives the required estimate. To show this, we have

\begin{array}{l} Ψ^{\pm} (0, t) = & Ct \pm (u (0, t) - φ_{b} (0, 0)) = Ct \pm (φ_{l} (t) - φ_{l} (0)) \\ = & Ct \pm Ct \geq 0, \end{array}

\begin{array}{l} Ψ^{\pm} (1, t) = & Ct \pm (u (1, t) - φ_{b} (1, 0)) = Ct \pm (φ_{r} (t) - φ_{r} (1)) \\ = & Ct \pm Ct \geq 0. \end{array}

On the boundary $Γ_{b},$

\begin{array}{l} Ψ^{\pm} (x, t) = & Ct \pm (u (x, t) - φ_{b} (x, t)) \\ = & Ct \pm (φ_{b} (x, t) - φ_{b} (x, t)) = Ct \geq 0. \end{array}

Additionally, for $(x, t) \in Ω$

\begin{array}{l} L Ψ^{\pm} (x, t) = & L (Ct \pm (u (x, t) - φ_{b} (x, t))) \\ = & C L (t) \pm L (u (x, t) - φ_{b} (x, t)) \geq 0, \end{array}

since C can be chosen sufficiently large.

Therefore, using the maximum principle (Lemma 1.1), we get

\begin{array}{l} Ψ^{\pm} (x, t) \geq 0, \forall (x, t) \in \overset{―}{Ω} . \end{array}

This implies that

| u (x, t) - φ_{b} (x, 0) | \leq Ct, (x, t) \in \overset{―}{Ω} .

Lemma 1.3.

The solution $u (x, t)$ of the problem (1)-(2) satisfy the following bound:

\begin{array}{l} | u (x, t) | \leq C, \forall (x, t) \in \overset{―}{Ω} . \end{array}

Proof.

For all $(x, t) \in \overset{―}{Ω}$ , we have

\begin{array}{l} | u (x, t) | = & | u (x, t) - φ_{b} (x, 0) + φ_{b} (x, 0) | \leq | u (x, t) - φ_{b} (x, 0) | \\ + | φ_{b} (x, 0) | . \end{array}

Using Lemma 1.2 and the fact that $| φ_{b} (x, 0) |$ is bounded, we get the required result.

Lemma 1.4.

(Stability Estimate of the continuous problem) Let $u (x, t)$ be the solution of problem (1)-(2). Then, we have

\begin{array}{l} ‖ u ‖_{\overset{―}{Ω}} \leq ‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} . \end{array}

Proof.

Define the barrier functions

\begin{array}{l} Φ^{\pm} (x, t) = ‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} \pm u (x, t) . \end{array}

Thus we have

\begin{array}{l} Φ^{\pm} (0, t) = ‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} \pm u (0, t) \geq ‖ u ‖_{Γ} \pm u (0, t) \geq 0, \end{array}

\begin{array}{l} Φ^{\pm} (1, t) = ‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} \pm u (1, t) \geq ‖ u ‖_{Γ} \pm u (1, t) \geq 0, \end{array}

Also, for $(x, t) \in Γ_{b},$

\begin{array}{l} Φ^{\pm} (x, t) = ‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} \pm u (x, t) \geq ‖ u ‖_{Γ} \pm u (x, t) \geq 0, \end{array}

Moreover, for $(x, t) \in Ω,$

\begin{array}{l} \begin{array}{l} L Φ^{\pm} (x, t) & = L (‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}} \pm u (x, t)) \\ = q (‖ u ‖_{Γ} + α^{- 1} ‖ L u ‖_{\overset{―}{Ω}}) \pm L u (x, t) \\ \geq α ‖ u ‖_{Γ} + ‖ L u ‖_{\overset{―}{Ω}} \pm L u (x, t) \\ \geq ‖ L u ‖_{\overset{―}{Ω}} \pm L u (x, t) \geq 0 \end{array} \end{array}

Therefore, by using Lemma 1.1, we get the required result.

2. Numerical scheme

In this section, we present parameter uniformly convergent numerical scheme to approximate the problem given in (1)-(2) by discretizing both the time and space variables with a uniform mesh. To design the scheme, we first discretize the temporal variable using the Crank–Nicolson method. Next, the semi-discretized problem is numerically approximated by using the non-standard finite difference method with a simple upwind finite difference method in the spatial direction.

2.1. Time discretization

The use of Taylor’s series expansion for the term containing the large delay η in EquationEquation (1)(1) $\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}$ (1) is invalid because of Taylor’s series approximation gives better results when η is small. To address this challenge, first we discretize the time domain $[0, T]$ into the N subinterval with a uniform step size $Δt = \frac{T}{N}$ and by choosing K in such a way that $η = K Δt$ for some positive integer $K \in (0, N) .$ Thus, we have

Ω_{t}^{N} = {0 = t_{0} < t_{1} < \dots < t_{K} = η < \dots < t_{N - 1} < t_{N} = T},

Ω_{η}^{K} = {- η = t_{- K} < t_{- K + 1} < \dots < t_{- 1} < t_{0} = 0} .

Now we semi-discretize the problem (1) on $Ω_{t}^{N}$ as follows:

(3)

\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) - ϵ \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial x^{2}} + p (x) \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂x} \\ + q (x, t_{j + \frac{1}{2}}) u (x, t_{j + \frac{1}{2}}) = - r (x, t_{j + \frac{1}{2}}) u (x, t_{j + \frac{1}{2} - K}) + g (x, t_{j + \frac{1}{2}}) \end{array}

(3)

By Taylor series expansion with center at the point $(x, t_{j + \frac{1}{2}}),$ we have

(4)

\begin{array}{l} u (x, t_{j + 1}) = & u (x, t_{j + \frac{1}{2}}) + \frac{Δt}{2} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} \\ + {(\frac{Δt}{2})}^{2} \frac{1}{2!} \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial t^{2}} + O ((Δt)^{3}) \end{array}

(4)

and

(5)

\begin{array}{l} u (x, t_{j}) = & u (x, t_{j + \frac{1}{2}}) - \frac{Δt}{2} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} \\ + {(\frac{Δt}{2})}^{2} \frac{1}{2!} \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial t^{2}} - O ((Δt)^{3}) \end{array}

(5)

Subtracting EquationEquation (5)(5) $\begin{array}{l} u (x, t_{j}) = & u (x, t_{j + \frac{1}{2}}) - \frac{Δt}{2} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} \\ + {(\frac{Δt}{2})}^{2} \frac{1}{2!} \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial t^{2}} - O ((Δt)^{3}) \end{array}$ (5) from EquationEquation (4)(4) $\begin{array}{l} u (x, t_{j + 1}) = & u (x, t_{j + \frac{1}{2}}) + \frac{Δt}{2} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} \\ + {(\frac{Δt}{2})}^{2} \frac{1}{2!} \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial t^{2}} + O ((Δt)^{3}) \end{array}$ (4) , we get the following:

(6)

\frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} = \frac{u (x, t_{j + 1}) - u (x, t_{j})}{Δt} + TE, where

(6)

\begin{array}{l} TE = - \frac{(Δt)^{2}}{24} \frac{\partial^{3} u (x, t_{j + \frac{1}{2}})}{\partial t^{3}} . \end{array}

Next, we approximate the Caputo fractional derivative $_{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}})$ in the following approach.

\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) & = \frac{1}{Γ (1 - β)} \int_{0}^{t_{j + \frac{1}{2}}} (t_{j + \frac{1}{2}} - s)^{- β} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂s} ds \\ = \frac{1}{Γ (1 - β)} \int_{0}^{t_{j}} (t_{j + \frac{1}{2}} - s)^{- β} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂s} ds \\ + \frac{1}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} (t_{j + \frac{1}{2}} - s)^{- β} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂s} ds \\ = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} \int_{t_{n}}^{t_{n + 1}} (t_{j + \frac{1}{2}} - s)^{- β} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂s} ds \\ + \frac{1}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} (t_{j + \frac{1}{2}} - s)^{- β} \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂s} ds \end{array}

Using a second-order central finite difference approximation for the first-order derivative given in EquationEquation (6)(6) $\frac{\partial u (x, t_{j + \frac{1}{2}})}{∂t} = \frac{u (x, t_{j + 1}) - u (x, t_{j})}{Δt} + TE, where$ (6) , we obtain the following equation.

(7)

\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) & = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} \int_{t_{n}}^{t_{n + 1}} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt} + O ((Δt)^{2})] \\ {((j + \frac{1}{2}) Δt - s)}^{- β} ds + \frac{1}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt} \\ + O ((Δt)^{2})] {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt} + O ((Δt)^{2})] \\ \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds + \frac{1}{Γ (1 - β)} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt} \\ + O ((Δt)^{2})] \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt}] \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{O ((Δt)^{2})}{Γ (1 - β)} \sum_{n = 0}^{j - 1} \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{1}{Γ (1 - β)} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt}] \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{O ((Δt)^{2})}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{(Δt)^{- β}}{Γ (2 - β)} \sum_{n = 0}^{j - 1} [(u^{n + 1} (x) - u^{n} (x)) ({(j - n + \frac{1}{2})}^{1 - β} \\ - (j - n - \frac{1}{2})^{1 - β})] \\ + \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} [{(j - n + \frac{1}{2})}^{1 - β} - {(j - n - \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ + \frac{(Δt)^{- β}}{Γ (2 - β)} (\frac{u^{j + 1} (x) - u^{j} (x)}{2^{1 - β}}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}

(7)

Rearranging the above row of EquationEquation (7)(7) $\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) & = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} \int_{t_{n}}^{t_{n + 1}} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt} + O ((Δt)^{2})] \\ {((j + \frac{1}{2}) Δt - s)}^{- β} ds + \frac{1}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt} \\ + O ((Δt)^{2})] {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt} + O ((Δt)^{2})] \\ \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds + \frac{1}{Γ (1 - β)} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt} \\ + O ((Δt)^{2})] \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{1}{Γ (1 - β)} \sum_{n = 0}^{j - 1} [\frac{u^{n + 1} (x) - u^{n} (x)}{Δt}] \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{O ((Δt)^{2})}{Γ (1 - β)} \sum_{n = 0}^{j - 1} \int_{t_{n}}^{t_{n + 1}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{1}{Γ (1 - β)} [\frac{u^{j + 1} (x) - u^{j} (x)}{Δt}] \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ + \frac{O ((Δt)^{2})}{Γ (1 - β)} \int_{t_{j}}^{t_{j + \frac{1}{2}}} {((j + \frac{1}{2}) Δt - s)}^{- β} ds \\ = \frac{(Δt)^{- β}}{Γ (2 - β)} \sum_{n = 0}^{j - 1} [(u^{n + 1} (x) - u^{n} (x)) ({(j - n + \frac{1}{2})}^{1 - β} \\ - (j - n - \frac{1}{2})^{1 - β})] \\ + \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} [{(j - n + \frac{1}{2})}^{1 - β} - {(j - n - \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ + \frac{(Δt)^{- β}}{Γ (2 - β)} (\frac{u^{j + 1} (x) - u^{j} (x)}{2^{1 - β}}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}$ (7) , we get an equivalent equation of the form:

(8)

\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & \frac{(Δt)^{- β}}{Γ (2 - β)} (\frac{u^{j + 1} (x) - u^{j} (x)}{2^{1 - β}}) + \frac{(Δt)^{- β}}{Γ (2 - β)} \\ \times \sum_{n = 0}^{j - 1} [(u^{n + 1} (x) - u^{n} (x)) ((j - n + \frac{1}{2})^{1 - β} \\ - {(j - n - \frac{1}{2})}^{1 - β})] \frac{1}{Γ (2 - β)} \\ \times \sum_{n = 0}^{j - 1} [(j - n + \frac{1}{2})^{1 - β} - {(j - n - \frac{1}{2})}^{1 - β}] \\ \times O ((Δt)^{3 - β}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}

(8)

Let

(9)

ϑ = \frac{Δ t^{- β}}{Γ (2 - β)},

(9)

(10)

σ = \frac{Δ t^{- β}}{Γ (2 - β) 2^{(1 - β)}}, and

(10)

(11)

b_{j} = {(j + \frac{1}{2})}^{1 - β} - {(j - \frac{1}{2})}^{1 - β} .

(11)

Substituting EquationEquations (9)(9) $ϑ = \frac{Δ t^{- β}}{Γ (2 - β)},$ (9) , (Equation10(10) $σ = \frac{Δ t^{- β}}{Γ (2 - β) 2^{(1 - β)}}, and$ (10) ) and (Equation11(11) $b_{j} = {(j + \frac{1}{2})}^{1 - β} - {(j - \frac{1}{2})}^{1 - β} .$ (11) ) into EquationEquation (8)(8) $\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & \frac{(Δt)^{- β}}{Γ (2 - β)} (\frac{u^{j + 1} (x) - u^{j} (x)}{2^{1 - β}}) + \frac{(Δt)^{- β}}{Γ (2 - β)} \\ \times \sum_{n = 0}^{j - 1} [(u^{n + 1} (x) - u^{n} (x)) ((j - n + \frac{1}{2})^{1 - β} \\ - {(j - n - \frac{1}{2})}^{1 - β})] \frac{1}{Γ (2 - β)} \\ \times \sum_{n = 0}^{j - 1} [(j - n + \frac{1}{2})^{1 - β} - {(j - n - \frac{1}{2})}^{1 - β}] \\ \times O ((Δt)^{3 - β}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}$ (8) , we get

(12)

\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & σ (u^{j + 1} (x) - u^{j} (x)) \\ + ϑ \sum_{n = 0}^{j - 1} (u^{n + 1} (x) - u^{n} (x)) b_{j - n} \\ + \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} b_{j - n} O ((Δt)^{3 - β}) \\ + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}

(12)

EquationEquation (12)(12) $\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & σ (u^{j + 1} (x) - u^{j} (x)) \\ + ϑ \sum_{n = 0}^{j - 1} (u^{n + 1} (x) - u^{n} (x)) b_{j - n} \\ + \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} b_{j - n} O ((Δt)^{3 - β}) \\ + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \end{array}$ (12) can be equivalently written as

(13)

\begin{array}{l} \begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) & = σ (u^{j + 1} (x) - u^{j} (x)) + ϑ [- b_{j} u^{0} (x) \\ + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) + b_{1} u^{j} (x)] + R, \\ where \end{array} \end{array}

(13)

(14)

\begin{array}{l} \begin{array}{l} R & = \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} b_{j - n} O ((Δt)^{3 - β}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [\sum_{n = 0}^{j - 1} b_{j - n} + \frac{1}{2^{1 - β}}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(j + \frac{1}{2})}^{1 - β} - {(\frac{1}{2})}^{1 - β} + \frac{1}{2^{1 - β}}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(j + \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(\frac{t_{j}}{Δt} + \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ \leq C {(Δt)}^{2} . \end{array} \end{array}

(14)

Hence, the Crank–Nicolson approximation of the Caputo fractional derivative at $(x, t_{j + \frac{1}{2}})$ is given by

(15)

\begin{array}{l} \begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & ϑ [b_{1} u^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) \\ - b_{j} u^{0} (x)] + σ (u^{j + 1} (x) - u^{j} (x)) \\ + O ((Δt)^{2}) . \end{array} \end{array}

(15)

Thus, from EquationEquation (15)(15) $\begin{array}{l} \begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & ϑ [b_{1} u^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) \\ - b_{j} u^{0} (x)] + σ (u^{j + 1} (x) - u^{j} (x)) \\ + O ((Δt)^{2}) . \end{array} \end{array}$ (15) , we observed a second-order approximation to $_{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) .$ This means that the Crank–Nicolson scheme is a second-order convergence for the time-fractional derivative.

Substituting EquationEquation (15)(15) $\begin{array}{l} \begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) = & ϑ [b_{1} u^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) \\ - b_{j} u^{0} (x)] + σ (u^{j + 1} (x) - u^{j} (x)) \\ + O ((Δt)^{2}) . \end{array} \end{array}$ (15) by neglecting the truncation error into EquationEquation (3)(3) $\begin{array}{l} _{0} D_{t}^{β} u (x, t_{j + \frac{1}{2}}) - ϵ \frac{\partial^{2} u (x, t_{j + \frac{1}{2}})}{\partial x^{2}} + p (x) \frac{\partial u (x, t_{j + \frac{1}{2}})}{∂x} \\ + q (x, t_{j + \frac{1}{2}}) u (x, t_{j + \frac{1}{2}}) = - r (x, t_{j + \frac{1}{2}}) u (x, t_{j + \frac{1}{2} - K}) + g (x, t_{j + \frac{1}{2}}) \end{array}$ (3) , we have the following approximation formula:

(16)

\begin{array}{l} \begin{array}{l} ϑ [b_{1} {\tilde{u}}^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) {\tilde{u}}^{n} (x) - b_{j} {\tilde{u}}^{0} (x)] \\ + σ ({\tilde{u}}^{j + 1} (x) - {\tilde{u}}^{j} (x)) - ϵ \frac{\partial^{2} {\tilde{u}}^{j + \frac{1}{2}} (x)}{\partial x^{2}} \\ + p (x) \frac{\partial {\tilde{u}}^{j + \frac{1}{2}} (x)}{∂x} + q^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2}} (x) \\ = - r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) + g^{j + \frac{1}{2}} (x), \end{array} \end{array}

(16)

where ${\tilde{u}}^{j + \frac{1}{2}} (x)$ represents the approximate solution of the exact solution $u (x, t_{j + \frac{1}{2}})$ at ${(j + \frac{1}{2})}^{th}$ time level. Now, introducing the notation ${\tilde{u}}^{j + \frac{1}{2}} (x) = \frac{\tilde{u} (x, t_{j + 1}) + \tilde{u} (x, t_{j})}{2}$ into the above Equation and rearranging, we get the following equivalent form:

(17)

\begin{array}{l} \begin{array}{l} - ϵ \frac{\partial^{2} {\tilde{u}}^{j + 1} (x)}{\partial x^{2}} + p (x) \frac{\partial {\tilde{u}}^{j + 1} (x)}{∂x} + q^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + 1} (x) \\ = ϵ \frac{\partial^{2} {\tilde{u}}^{j} (x)}{\partial x^{2}} - p (x) \frac{\partial {\tilde{u}}^{j} (x)}{∂x} - q^{j + \frac{1}{2}} (x) {\tilde{u}}^{j} (x) \\ - 2 r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) + 2 g^{j + \frac{1}{2}} (x) \\ - 2 ϑ [b_{1} u^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) - b_{j} u^{0} (x)] \\ - 2 σ (u^{j + 1} (x) - u^{j} (x)), \end{array} \end{array}

(17)

where

\begin{array}{l} q^{j + \frac{1}{2}} (x) = \frac{q^{j + 1} (x) + q^{j} (x)}{2}, \end{array}

\begin{array}{l} r^{j + \frac{1}{2}} (x) = \frac{r^{j + 1} (x) + r^{j} (x)}{2}, \end{array}

\begin{array}{l} g^{j + \frac{1}{2}} (x) = \frac{g^{j + 1} (x) + g^{j} (x)}{2} . \end{array}

Rearranging EquationEquation (17)(17) $\begin{array}{l} \begin{array}{l} - ϵ \frac{\partial^{2} {\tilde{u}}^{j + 1} (x)}{\partial x^{2}} + p (x) \frac{\partial {\tilde{u}}^{j + 1} (x)}{∂x} + q^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + 1} (x) \\ = ϵ \frac{\partial^{2} {\tilde{u}}^{j} (x)}{\partial x^{2}} - p (x) \frac{\partial {\tilde{u}}^{j} (x)}{∂x} - q^{j + \frac{1}{2}} (x) {\tilde{u}}^{j} (x) \\ - 2 r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) + 2 g^{j + \frac{1}{2}} (x) \\ - 2 ϑ [b_{1} u^{j} (x) + \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) u^{n} (x) - b_{j} u^{0} (x)] \\ - 2 σ (u^{j + 1} (x) - u^{j} (x)), \end{array} \end{array}$ (17) , we get the following equivalent form.

(18)

\begin{array}{l} \begin{array}{l} - ϵ \frac{\partial^{2} {\tilde{u}}^{j + 1} (x)}{\partial x^{2}} + p (x) \frac{\partial {\tilde{u}}^{j + 1} (x)}{∂x} + (q^{j + \frac{1}{2}} (x) + 2 σ) {\tilde{u}}^{j + 1} (x) \\ = ϵ \frac{\partial^{2} {\tilde{u}}^{j} (x)}{\partial x^{2}} - p (x) \frac{\partial {\tilde{u}}^{j} (x)}{∂x} - (q^{j + \frac{1}{2}} (x) - 2 σ + 2 ϑ b_{1}) {\tilde{u}}^{j} (x) \\ - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) {\tilde{u}}^{n} (x) + 2 ϑ b_{j} {\tilde{u}}^{0} (x) \\ - 2 r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) + 2 g^{j + \frac{1}{2}} (x), x \in Ω_{x}, j = 0, 1, \dots, N - 1. \end{array} \end{array}

(18)

Using discrete operator, EquationEquation (18)(18) $\begin{array}{l} \begin{array}{l} - ϵ \frac{\partial^{2} {\tilde{u}}^{j + 1} (x)}{\partial x^{2}} + p (x) \frac{\partial {\tilde{u}}^{j + 1} (x)}{∂x} + (q^{j + \frac{1}{2}} (x) + 2 σ) {\tilde{u}}^{j + 1} (x) \\ = ϵ \frac{\partial^{2} {\tilde{u}}^{j} (x)}{\partial x^{2}} - p (x) \frac{\partial {\tilde{u}}^{j} (x)}{∂x} - (q^{j + \frac{1}{2}} (x) - 2 σ + 2 ϑ b_{1}) {\tilde{u}}^{j} (x) \\ - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) {\tilde{u}}^{n} (x) + 2 ϑ b_{j} {\tilde{u}}^{0} (x) \\ - 2 r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) + 2 g^{j + \frac{1}{2}} (x), x \in Ω_{x}, j = 0, 1, \dots, N - 1. \end{array} \end{array}$ (18) can be written as

(19)

L^{N} {\tilde{u}}^{j + 1} (x) = F (x, t_{j + 1}), x \in Ω_{x}, j = 0, 1, \dots, N - 1,

(19)

where

(20)

\begin{array}{l} \begin{array}{l} F (x, t_{j + 1}) = & ϵ \frac{\partial^{2} {\tilde{u}}^{j} (x)}{\partial x^{2}} - p (x) \frac{\partial {\tilde{u}}^{j} (x)}{∂x} \\ - (q^{j + \frac{1}{2}} (x) - 2 σ + 2 ϑ b_{1}) {\tilde{u}}^{j} (x) \\ - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) {\tilde{u}}^{n} (x) \\ + 2 ϑ b_{j} {\tilde{u}}^{0} (x) - 2 r^{j + \frac{1}{2}} (x) {\tilde{u}}^{j + \frac{1}{2} - K} (x) \\ + 2 g^{j + \frac{1}{2}} (x) \end{array} \end{array}

(20)

with the interval boundary conditions:

(21)

\begin{array}{l} \begin{array}{l} {\tilde{u}}^{j + 1} (x) & = φ_{b} (x, t_{j + 1}), x \in Ω_{x}, - (K + 1) \leq j \leq - 1, \\ {\tilde{u}}^{j + 1} (0) & = φ_{l} (t_{j + 1}), {\tilde{u}}^{j + 1} (1) = φ_{r} (t_{j + 1}) for 0 \leq j \leq N - 1. \end{array} \end{array}

(21)

The operator $L^{N}$ is defined as

(22)

\begin{array}{l} L^{N} {\tilde{u}}^{j + 1} (x) \equiv & - ϵ \frac{\partial^{2} {\tilde{u}}^{j + 1} (x)}{\partial x^{2}} + p (x) \frac{\partial {\tilde{u}}^{j + 1} (x)}{∂x} \\ + (q^{j + \frac{1}{2}} (x) + 2 σ) {\tilde{u}}^{j + 1} (x), \end{array}

(22)

satisfies the next semi-discrete maximum principle.

Lemma 2.1.

(The Semi-discrete maximum principle) Let ${\tilde{Φ}}^{j + 1} (x)$ be a smooth function on ${\overset{―}{Ω}}_{x}$ . If ${\tilde{Φ}}^{j + 1} (0) \geq 0, \break {\tilde{Φ}}^{j + 1} (1) \geq 0$ and $L^{N} {\tilde{Φ}}^{j + 1} (x) \geq 0, ∀x \in Ω_{x},$ then ${\tilde{Φ}}^{j + 1} (x) \geq 0, ∀x \in {\overset{―}{Ω}}_{x} .$

Proof.

To proof is done by contradiction. Let $x^{⋆} \in {\overset{―}{Ω}}_{x}$ such that

\begin{array}{l} {\tilde{Φ}}^{j + 1} (x^{⋆}) = min_{x \in {\overset{―}{Ω}}_{x}} {\tilde{Φ}}^{j + 1} (x) \end{array}

and assume that ${\tilde{Φ}}^{j + 1} (x^{⋆}) < 0$ . From the hypotheses, we have $x^{⋆} \neq {0, 1}$ and $\frac{\partial {\tilde{Φ}}^{j + 1} (x^{⋆})}{∂x} = 0, \frac{\partial^{2} {\tilde{Φ}}^{j + 1} (x^{⋆})}{\partial x^{2}} \geq 0$ , hence

\begin{array}{l} L^{N} {\tilde{Φ}}^{j + 1} (x^{⋆}) = - ϵ \frac{\partial^{2} {\tilde{Φ}}^{j + 1} (x^{⋆})}{\partial x^{2}} + p (x^{⋆}) \frac{\partial {\tilde{Φ}}^{j + 1} (x^{⋆})}{∂x} \\ + (q^{j + \frac{1}{2}} (x^{⋆}) + 2 σ) {\tilde{Φ}}^{j + 1} (x^{⋆}) < 0, \end{array}

which is a contradiction to our assumption of $L^{N} {\tilde{Φ}}^{j + 1} (x^{⋆}) \geq 0$ in $Ω_{x}$ . This shows that our assumption is wrong and then ${\tilde{Φ}}^{j + 1} (x^{⋆}) \geq 0.$ Therefore, we conclude that ${\tilde{Φ}}^{j + 1} (x^{⋆}) \geq 0, \break \forall x \in {\overset{―}{Ω}}_{x}$ .

Theorem 2.2.

Under the smoothness and compatibility conditions, the semi-discretized solution ${\tilde{u}}^{j + 1} (x)$ and its derivatives satisfy the following bounds:

\begin{array}{l} | \frac{d^{i} {\tilde{u}}^{j + 1} (x)}{d x^{i}} | \leq & C (1 + ϵ^{- i} \exp (- \frac{α (1 - x)}{ϵ})), \\ i = 0, 1, 2, 3, 4. \end{array}

Proof.

The proof is found in (Kellogg et al., Citation1978).

2.2. Spatial discretization

Now, we consider the spatial discretization of EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) by the non-standard finite difference method on a uniform mesh. Let M be a positive integer and divide the domain $[0, 1]$ into M equal sub-intervals with step size $Δx = 1 / M$ :

Ω_{x}^{M} = {0 = x_{0} < x_{1} < \dots < x_{M} = 1},

\begin{array}{l} x_{i} = x_{0} + i Δx, i = 0, 1, \dots, M . \end{array}

To discretize the space derivatives using the method of non-standard finite difference scheme on a uniform mesh, we will follow the procedures used in (Mickens Citation2005). Using these principles, the denominator function denoted by $ϕ_{i}^{2}$ is given by (Mbroh et al., Citation2019)

(23)

ϕ_{i}^{2} (Δx, ϵ) = \frac{ϵ Δx}{p (x_{i})} (\exp (\frac{p (x_{i}) Δx}{ϵ}) - 1) .

(23)

Let $U_{i, j}$ represent the approximation of $\tilde{u} (x_{i}, t_{j}) .$ Using this denominator function and the simple upwind finite difference approximation for the first derivative of EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) , we obtain the following fully-discrete problem.

(24)

\begin{array}{l} L^{M, N} U_{i, j + 1} \equiv & - ϵ [\frac{U_{i - 1, j + 1} - 2 U_{i, j + 1} + U_{i + 1, j + 1}}{ϕ_{i}^{2}}] \\ + p_{i} [\frac{U_{i, j + 1} - U_{i - 1, j + 1}}{Δx}] \\ + [\frac{1}{2} (q_{i, j} + q_{i, j + 1}) + 2 σ] U_{i, j + 1} \\ = & ϵ [\frac{U_{i - 1, j} - 2 U_{i, j} + U_{i + 1, j}}{ϕ_{i}^{2}}] - p_{i} [\frac{U_{i, j} - U_{i - 1, j}}{Δx}] \\ - [\frac{1}{2} (q_{i, j} + q_{i, j + 1}) - 2 σ + 2 ϑ b_{1}] U_{i, j} \\ - \frac{1}{2} [r_{i, j} + r_{i, j + 1}] [U_{i, j + 1 - K} + U_{i, j - K}] \\ + (g_{i, j} + g_{i, j + 1}) \\ - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) U_{i, n} + 2 ϑ b_{j} U_{i, 0}, \\ i & = 1, 2, \dots, M - 1, j = 0, 1, 2, \dots, N - 1, \end{array}

(24)

with the interval boundary conditions:

(25)

\begin{array}{l} U (x_{i}, t_{j + 1}) = φ_{b} (x_{i}, t_{j + 1}), i = 0, 1, \dots, \\ M, j = - (K + 1), - K + 1, \dots, - 1, \\ U (0, t_{j + 1}) = φ_{l} (t_{j + 1}), U (1, t_{j + 1}) = φ_{r} (t_{j + 1}), \\ j = 0, 1, \dots, N - 1. \end{array}

(25)

After making some rearrangement, we have the following equation.

(26)

\begin{array}{l} - (\frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}) U_{i - 1, j + 1} + (\frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) + 2 σ) \\ U_{i, j + 1} - (\frac{ϵ}{ϕ_{i}^{2}}) U_{i + 1, j + 1} = (\frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}) U_{i - 1, j} \\ - (\frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) - 2 σ + 2 ϑ b_{1}) U_{i, j} \\ + (\frac{ϵ}{ϕ_{i}^{2}}) U_{i + 1, j} - \frac{1}{2} [r_{i, j} + r_{i, j + 1}] [U_{i, j + 1 - K} + U_{i, j - K}] \\ + (g_{i, j} + g_{i, j + 1}) - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) U_{i, n} + 2 ϑ b_{j} U_{i, 0}, \\ i = 1, 2, \dots, M - 1, j = 0, 1, 2, \dots, N - 1. \end{array}

(26)

Equivalently, we write EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) in linear system form as follows:

(27)

\begin{array}{l} A_{i} U_{i - 1, j + 1} + B_{i} U_{i, j + 1} + C_{i} U_{i + 1, j + 1} = & D_{i} U_{i - 1, j} + E_{i} U_{i, j} \\ + F_{i} U_{i + 1, j} + G_{i, j + 1}, \end{array}

(27)

where

(28)

\begin{array}{l} A_{i} = & - (\frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}), \\ B_{i} = & \frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) + 2 σ, \\ C_{i} = & - \frac{ϵ}{ϕ_{i}^{2}}, D_{i} = \frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}, \\ E_{i} = & - (\frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) - 2 σ + 2 ϑ b_{1}), \\ F_{i} = & \frac{ϵ}{ϕ_{i}^{2}}, G_{i, j + 1} = - \frac{1}{2} [r_{i, j} + r_{i, j + 1}] [U_{i, j + 1 - K} + U_{i, j - K}] \\ + (g_{i, j} + g_{i, j + 1}) - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) U_{i, n} + 2 ϑ b_{j} U_{i, 0} \end{array}

(28)

Since we let $q (x, t) > 0$ in the domain $Ω$ , again let us suppose that $v_{i, j} = \frac{1}{2} (q_{i, j} + q_{i, j + 1}) + 2 σ \geq ξ > 0.$

Lemma 2.3.

(The Fully-discrete Maximum Principle). Assume that the discrete function $Ψ (x_{i}, t_{j})$ satisfying $Ψ (x_{i}, t_{j}) \geq 0 on (x_{i}, t_{j}) \in Γ^{M, N}$ and $L^{M, N} Ψ (x_{i}, t_{j}) \geq 0 \breakin (x_{i}, t_{j}) \in Ω^{M, N}$ . Then $Ψ (x_{i}, t_{j}) \geq 0,$ for all $(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}$ .

Proof.

Let $s, s^{⋆}$ be the indexes such that

\begin{array}{l} Ψ_{s, s^{⋆}} = min_{(x_{i}, t_{j}) \in {\bar{Ω}}^{M, N}} Ψ (x_{i}, t_{j}) < 0. \end{array}

Then by the given hypothesis, we have $0 < s < M \breakand 0 < s^{⋆} < N$ . It follows that $Ψ_{s + 1, s^{⋆}} - Ψ_{s, s^{⋆}} \geq 0, \break Ψ_{s, s^{⋆}} - Ψ_{s - 1, s^{⋆}} \leq 0.$ Hence,

\begin{array}{l} L^{M, N} Ψ_{s, s^{⋆}} = & - ϵ [\frac{Ψ_{s + 1, s^{⋆}} - 2 Ψ_{s, s^{⋆}} + Ψ_{s - 1, s^{⋆}}}{ϕ_{s}^{2}}] \\ + p_{s} [\frac{Ψ_{s, s^{⋆}} - Ψ_{s - 1, s^{⋆}}}{Δx}] \\ + v_{s, s^{⋆}} Ψ_{s, s^{⋆}} = - ϵ [\frac{Ψ_{s + 1, s^{⋆}} - Ψ_{s, s^{⋆}}}{ϕ_{s}^{2}} - \frac{Ψ_{s, s^{⋆}} - Ψ_{s - 1, s^{⋆}}}{ϕ_{s}^{2}}] \\ + p_{s} [\frac{Ψ_{s, s^{⋆}} - Ψ_{s - 1, s^{⋆}}}{Δx}] + v_{s, s^{⋆}} Ψ_{s, s^{⋆}} \leq 0, \\ which is a contradiction to our assumption. \end{array}

This shows that our assumption is wrong and $Ψ_{s, s^{⋆}} \geq 0.$ Therefore, we have $Ψ_{i, j} \geq 0, \forall (x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}$ .

An immediate consequence of Lemma 2.3 is the uniform stability estimate given below.

Lemma 2.4.

The solution $U_{i, j}$ of the fully-discretized problem in EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) satisfy the following bound.

\begin{array}{l} | U_{i, j} | \leq & ξ^{- 1} max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} | L^{M, N} U_{i, j} | \\ + max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} {| ψ_{b} (x_{i}, t_{j}) |, max_{t_{j} \in {\overset{―}{Ω}}_{t}^{N}} {| ψ_{l} (t) |, | ψ_{r} (t) |}} . \end{array}

Proof.

Let us denote

\begin{array}{l} κ = & ξ^{- 1} max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} | L^{M, N} U_{i, j} | \\ + max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} {| ψ_{b} (x_{i}, t_{j}) |, max_{t_{j} \in {\overset{―}{Ω}}_{t}^{N}} {| ψ_{l} (t) |, | ψ_{r} (t) |}} \end{array}

and define the barrier functions $Ψ_{i, j}^{\pm} = κ \pm U_{i, j} .$

At the boundary points

\begin{array}{l} Ψ_{0, j}^{\pm} = κ \pm U_{0, j} = κ \pm U (0, t_{j}) = κ \pm ψ_{l} (t_{j}) \geq 0, \end{array}

\begin{array}{l} Ψ_{M, j}^{\pm} = κ \pm U_{M, j} = κ \pm U (1, t_{j}) = κ \pm ψ_{r} (t_{j}) \geq 0. \end{array}

On the discretized domain $x_{i}, i = 1, 2, \dots, N - 1,$ we have

\begin{array}{l} L^{M} Ψ_{i, j}^{\pm} = & - ϵ (\frac{κ \pm U_{i + 1, j} - 2 (κ \pm U_{i, j}) + κ \pm U_{i - 1, j}}{ϕ_{i}^{2}}) \\ + p_{i} (\frac{κ \pm U_{i, j} - (κ \pm U_{i - 1, j})}{Δx}) + v_{i, j} (κ \pm U_{i, j}) \\ = & v_{i, j} κ \pm L^{M, N} U_{i, j}, = v_{i, j} (ξ^{- 1} max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} | L^{M, N} U_{i, j} | \\ + max_{(x_{i}, t_{j}) \in {\overset{―}{D}}^{M, N}} {| ψ_{b} (x_{i}, t_{j}) |, max_{t_{j} \in {\overset{―}{Ω}}_{t}^{N}} {| ψ_{l} (t) |, | ψ_{r} (t) |}}) \\ \pm L^{M} U_{i, j} \geq 0. \end{array}

By Lemma 2.3, we obtain $Ψ_{i, j}^{\pm} \geq 0$ , for all $(x_{i}, t_{j}) \in {\bar{Ω}}^{M, N}$ . This implies that

\begin{array}{l} | U_{i, j} | \leq & ξ^{- 1} max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} | L^{M, N} U_{i, j} | + max_{(x_{i}, t_{j}) \in {\overset{―}{Ω}}^{M, N}} {| ψ_{b} (x_{i}, t_{j}) |, \\ max_{t_{j} \in {\overset{―}{Ω}}_{t}^{N}} {| ψ_{l} (t) |, | ψ_{r} (t) |}} . \end{array}

In the next section, we estimate the error associated with the spatial discretization.

3. Convergence analysis

By using the information presented above, we were able to demonstrate that the fully- discrete operator denoted by $L^{M, N}$ satisfies both the discrete maximum principle and the uniform stability estimate. Let $U_{i, j}$ denotes the approximate solution in the fully-discretization of the exact solution $u (x, t)$ at $(x_{i}, t_{j}), i = 0, 1, \dots, M, j = 0, 1, \dots, N .$ Next, let us denote and define the forward, backward and central finite difference approximations of the space derivatives, respectively, as

\begin{array}{l} D_{x}^{+} u_{i, j} = \frac{u_{i + 1, j} - u_{i, j}}{Δx}, D_{x}^{-} u_{i, j} = \frac{u_{i, j} - u_{i - 1, j}}{Δx}, and \end{array}

\begin{array}{l} δ_{x}^{2} u_{i, j} = D_{x}^{+} D_{x}^{-} u_{i, j} = \frac{D^{+} u_{i, j} - D^{-} u_{i, j}}{Δx} . \end{array}

Then, the truncation error of the scheme (24)-(25) is given by

\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | = | L^{M, N} u (x_{i}, t_{j}) - L^{M} U_{i, j} | \\ = | - ϵ \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} + p_{i} \frac{∂u (x_{i}, t_{j})}{∂x} \\ + ϵ [\frac{U_{i - 1, j} - 2 U_{i, j} + U_{i + 1, j}}{ϕ_{i}^{2}}] - p_{i} [\frac{U_{i, j} - U_{i - 1, j}}{Δx}] | \\ = ϵ | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} - \frac{D_{x}^{+} D_{x}^{-} (Δx)^{2} u (x_{i}, t_{j})}{ϕ_{i}^{2}} \\ + p_{i} (\frac{∂u (x_{i}, t_{j})}{∂x} - D_{x}^{-} u (x_{i}, t_{j})) | \\ \leq ϵ | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} - D_{x}^{+} D_{x}^{-} u (x_{i}, t_{j}) | \\ + Cϵ | (\frac{(Δx)^{2}}{ϕ_{i}^{2}} - 1) D_{x}^{+} D_{x}^{-} u (x_{i}, t_{j}) | \\ + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | \end{array}

Here, we have $ϵ | \frac{(Δx)^{2}}{ϕ_{i}^{2}} - 1 | \leq C Δx .$ To show this, let us define $λ = p_{i} Δx / ϵ, λ \in (0, \in fty) .$ Using EquationEquation (23)(23) $ϕ_{i}^{2} (Δx, ϵ) = \frac{ϵ Δx}{p (x_{i})} (\exp (\frac{p (x_{i}) Δx}{ϵ}) - 1) .$ (23) , we obtain the following relation:

\begin{array}{l} ϵ | \frac{(Δx)^{2}}{ϕ_{i}^{2}} - 1 | = p_{i} Δx | \frac{1}{\exp (λ) - 1} - \frac{1}{λ} | = p_{i} ΔxR (λ), \\ where \end{array}

\begin{array}{l} R (λ) = \frac{\exp (λ) - λ - 1}{λ (\exp (λ) - 1)} . \end{array}

Because of $lim_{λ \to 0} R (λ) = 1 / 2, lim_{λ \to \in fty} R (λ) = 0,$ we have $R (λ) \leq C, ∀λ \in (0, \in fty) .$ Using this relation, we have the following:

\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | \leq Cϵ (Δx)^{2} | \frac{\partial^{4} u (x_{i}, t_{j})}{\partial x^{4}} | \\ + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | \\ \leq Cϵ (Δx)^{2} | \frac{\partial^{4} u (x_{i}, t_{j})}{\partial x^{4}} | + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | \end{array}

So, the error estimate in the spatial discretization is bounded as

(29)

\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | \leq & Cϵ (Δx)^{2} | \frac{\partial^{4} u (x_{i}, t_{j})}{\partial x^{4}} | \\ + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | \end{array}

(29)

Using the boundedness of derivatives of solution given in Theorem 2.2 and EquationEquation (29)(29) $\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | \leq & Cϵ (Δx)^{2} | \frac{\partial^{4} u (x_{i}, t_{j})}{\partial x^{4}} | \\ + C Δx | \frac{\partial^{2} u (x_{i}, t_{j})}{\partial x^{2}} | \end{array}$ (29) , we obtain

(30)

\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | & \leq Cϵ (Δx)^{2} | 1 + ϵ^{- 4} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 2} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ \leq C (Δx)^{2} | ϵ + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 2} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ \leq C (Δx)^{2} | 1 + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) |, \\ \leq C Δx | 1 + max_{i} ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | . \end{array}

(30)

Lemma 3.1.

For a fixed mesh and ϵ → 0, it holds

(31)

lim_{ϵ \to 0} max_{i} \frac{\exp (- \frac{α (1 - x_{i})}{ϵ})}{ϵ^{m}} = 0, m = 1, 2, 3, \dots

(31)

where $x_{i} = i Δx, Δx = 1 / M,$ for all $i = 1, 2, \dots, M - 1.$

Proof.

The proof is found in (Turuna et al., Citation2020).

Using the result of Lemma 3.1 into EquationEquation (30)(30) $\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | & \leq Cϵ (Δx)^{2} | 1 + ϵ^{- 4} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 2} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ \leq C (Δx)^{2} | ϵ + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 2} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ \leq C (Δx)^{2} | 1 + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | \\ + C Δx | 1 + ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) |, \\ \leq C Δx | 1 + max_{i} ϵ^{- 3} \exp (- \frac{α (1 - x_{i})}{ϵ}) | . \end{array}$ (30) , we obtain

\begin{array}{l} | L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx . \end{array}

Define the barrier functions

\begin{array}{l} Φ_{i, j}^{\pm} = CΔx \pm (u (x_{i}, t_{j}) - U_{i, j}) . \end{array}

The values $Φ_{i, j}^{\pm}$ at boundaries are

Φ_{i, 0}^{\pm} = CΔx \pm (u (x_{i}, t_{0}) - U_{i, 0}) = C Δx \geq 0,

Φ_{0, j}^{\pm} = CΔx \pm (u (x_{0}, t_{j}) - U_{0, j}) = C Δx \geq 0,

Φ_{M, j}^{\pm} = CΔx \pm (u (x_{M}, t_{j}) - U_{M, j}) = C Δx \geq 0, and

\begin{array}{l} L^{M, N} Φ_{i, j}^{\pm} & = L^{M, N} (C Δx \pm (u (x_{i}, t_{j}) - U_{i, j})), \\ = L^{M, N} (C Δx) \pm L^{M, N} (u (x_{i}, t_{j}) - U_{i, j}), \\ \geq 0, 0 < i < M, 0 < j < N . \end{array}

Using maximum principle, we have

\begin{array}{l} Φ_{i, j}^{\pm} \geq 0, 0 \leq i \leq M, 0 \leq j \leq N . \end{array}

From this it follows that

(32)

| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .

(32)

Theorem 3.2.

(ϵ-uniform convergence) If $u (x_{i}, t_{j})$ be the exact solution of the continuous problem EquationEquation (1)(1) $\begin{array}{l} \begin{array}{l} L u (x, t) : & \equiv \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + p (x) \frac{∂u (x, t)}{∂x} \\ + q (x, t) u (x, t) \\ = - r (x, t) u (x, t - η) + g (x, t), \forall (x, t) \in Ω, \end{array} \end{array}$ (1) -EquationEquation (2)(2) $\begin{array}{l} \begin{array}{l} u (x, t) & = φ_{b} (x, t), \forall (x, t) \in [0, 1] \times [- η, 0], \\ u (0, t) & = φ_{l} (t), ∀t \in (0, T), \\ u (1, t) & = φ_{r} (t), ∀t \in (0, T), \end{array} \end{array}$ (2) and $U_{i, j}$ be the numerical solution of the fully discretized scheme EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) -EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) , then the following error estimate holds

sup_{0 < ϵ << 1} ‖ U (x_{i}, t_{j}) - u (x_{i}, t_{j}) ‖_{Ω_{x}^{M} \times Ω_{t}^{N}} \leq C (Δx + (Δt)^{2}) .

Proof.

From EquationEquation (14)(14) $\begin{array}{l} \begin{array}{l} R & = \frac{1}{Γ (2 - β)} \sum_{n = 0}^{j - 1} b_{j - n} O ((Δt)^{3 - β}) + \frac{1}{Γ (2 - β) 2^{1 - β}} O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [\sum_{n = 0}^{j - 1} b_{j - n} + \frac{1}{2^{1 - β}}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(j + \frac{1}{2})}^{1 - β} - {(\frac{1}{2})}^{1 - β} + \frac{1}{2^{1 - β}}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(j + \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ = \frac{1}{Γ (2 - β)} [{(\frac{t_{j}}{Δt} + \frac{1}{2})}^{1 - β}] O ((Δt)^{3 - β}) \\ \leq C {(Δt)}^{2} . \end{array} \end{array}$ (14) and EquationEquation (32)(32) $| u (x_{i}, t_{j}) - U_{i, j}) | \leq C Δx, 0 \leq i \leq M, 0 \leq j \leq N .$ (32) , we can obtain the above error bound. To show this,

\begin{array}{l} | U (x_{i}, t_{j}) - u (x_{i}, t_{j}) | & = | U (x_{i}, t_{j}) - \tilde{u} (x_{i}, t_{j}) + \tilde{u} (x_{i}, t_{j}) - u (x_{i}, t_{j}) | \\ \leq | U (x_{i}, t_{j}) - \tilde{u} (x_{i}, t_{j}) | + | \tilde{u} (x_{i}, t_{j}) - u (x_{i}, t_{j}) | \\ \leq C (Δx + (Δt)^{2}) . \end{array}

4. Numerical examples

To show the effectiveness of the proposed scheme given in EquationEquation (26)(26) $\begin{array}{l} - (\frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}) U_{i - 1, j + 1} + (\frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) + 2 σ) \\ U_{i, j + 1} - (\frac{ϵ}{ϕ_{i}^{2}}) U_{i + 1, j + 1} = (\frac{ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx}) U_{i - 1, j} \\ - (\frac{2 ϵ}{ϕ_{i}^{2}} + \frac{p_{i}}{Δx} + \frac{1}{2} (q_{i, j} + q_{i, j + 1}) - 2 σ + 2 ϑ b_{1}) U_{i, j} \\ + (\frac{ϵ}{ϕ_{i}^{2}}) U_{i + 1, j} - \frac{1}{2} [r_{i, j} + r_{i, j + 1}] [U_{i, j + 1 - K} + U_{i, j - K}] \\ + (g_{i, j} + g_{i, j + 1}) - 2 ϑ \sum_{n = 1}^{j - 1} (b_{j - n + 1} - b_{j - n}) U_{i, n} + 2 ϑ b_{j} U_{i, 0}, \\ i = 1, 2, \dots, M - 1, j = 0, 1, 2, \dots, N - 1. \end{array}$ (26) for the given problem (1), we consider the following two examples. Since the exact solution to these problems is not known, to determine the maximum pointwise error and rate of convergence, we use the double mesh principle (Doolan et al., Citation1980). To apply this technique, the space and time domains are divided into 2 M and 2 N equal sub-intervals respectively. Then, the maximum pointwise error and the convergence rate are obtained as follows. The maximum pointwise error is

(33)

E^{M, N} = max_{0 \leq i, j \leq M, N} | U^{M, N} (x_{i}, t_{j}) - U^{2 M, 2 N} (x_{2 i}, t_{2 j}) |,

(33)

where $U^{M, N} (x_{i}, t_{j})$ denotes the numerical solution obtained on a mesh containing M +1 points in the spatial direction and N +1 points in the temporal direction, and $U^{2 M, 2 N} (x_{i}, t_{j})$ is obtained by further dividing each sub-interval in space and time direction. In addition, the rate of convergence of the problems are obtained by the formula

(34)

R^{M, N} = \log_{2} (\frac{E^{M, N}}{E^{2 M, 2 N}}) .

(34)

Example 4.1.

Consider the following time-fractional singularly perturbed parabolic partial differential equation with delay η = 1:

\begin{array}{l} \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + (2 - x^{2}) \frac{∂u (x, t)}{∂x} \\ + xu (x, t) = - u (x, t - η) + 10 t^{2} \exp (- t) x (1 - x), \\ (x, t) \in Ω = (0, 1) \times (0, 2), \end{array}

with the interval and boundary conditions $φ_{b} (x, t) = 0, \break φ_{l} (t) = 0, φ_{r} (t) = 0.$

Example 4.2.

Consider the following time-fractional singularly perturbed parabolic partial differential equation with delay η = 1:

\begin{array}{l} \frac{\partial^{β} u (x, t)}{\partial t^{β}} - ϵ \frac{\partial^{2} u (x, t)}{\partial x^{2}} + (2 - x^{2}) \frac{∂u (x, t)}{∂x} \\ + ((x + 1) (t + 1)) u (x, t) = - u (x, t - η) \\ + 10 t^{2} \exp (- t) x (1 - x), (x, t) \in Ω = (0, 1) \times (0, 2), \end{array}

with the interval and boundary conditions $φ_{b} (x, t) = 0, φ_{l} (t) = 0, φ_{r} (t) = 0.$

For different values of $ϵ, N, M$ and fixed β = 0.8, the maximum pointwise error and the corresponding rate of convergence for Example 4.1 and Example 4.2 are given in Tables and respectively. The last two rows in these tables represent the ϵ- uniform error estimate $(E^{M, N})$ and the ϵ- uniform convergence rate $(R^{M, N}) .$ Similarly, for different values of $β, N$ and M with fixed $ϵ = 10^{- 10}$ , the maximum pointwise error and the rate of convergence are given in Tables and respectively. From the Tables , one can clearly observe that the maximum pointwise error decreases monotonically as the mesh size M and N increases, which confirms that the proposed scheme is ϵ - uniformly convergent. To show the existence of boundary layer near x = 1 for small ϵ (ϵ → 0), the surface plots of solutions of Example 4.1 and Example 4.2 are given in Figures and respectively for different values of β. In addition, to show the effect of ϵ on the width of the boundary layer, the surface plots of solutions of Example 4.1 and Example 4.2 are given in Figures and respectively. Similarly, to show the effects of the fractional order β on the solution of Example 4.1 and Example 4.2, the surface plots are given in Figures and respectively. Figure shows the log-log plot of the maximum pointwise error for the solutions for Example 4.1 and Example 4.2. From this figure, one can observe the parameter uniform convergence of the proposed method. Moreover, the formulated scheme is shown to be first-order convergence globally.

Figure 1. Surface plots of the numerical solution of example 4.1 for different values of β with $ϵ = 10^{- 10}, M = 80$ and N = 64.

Figure 2. Surface plots of the numerical solution of example 4.1 for different values of ϵ with $β = 0.4, M = 80$ and N = 64.

Figure 3. Surface plots of the numerical solution of example 4.2 for different values of β with $ϵ = 10^{- 10}, M = 80$ and N = 64.

Figure 4. Surface plots of the numerical solution of example 4.2 for different values of ϵ with $β = 0.6, M = 80$ and N = 64.

Figure 5. Loglog plot of maximum pointwise error for examples 4.1 and 4.2 for β = 0.8.

Table 1. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using fractional order β = 0.8 with different values of $ϵ .$

Display Table

Table 2. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using $ϵ = 10^{- 10}$ with different values of β

Display Table

Table 3. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using fractional order β = 0.8 with different values of ϵ

Display Table

Table 4. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using $ϵ = 10^{- 10}$ with different values of β

Display Table

5. Conclusion

In this paper, a parameter uniformly convergent numerical scheme for a class of time-fractional singularly perturbed partial differential equations with a large time delay is considered. The proposed numerical scheme consists of the Crank–Nicolson method for discretizing the time-fractional derivative, which is given in Caputo sense and the non-standard finite difference method together with simple upwind scheme to discretize the spatial derivatives on a uniform mesh. The method is shown to be second-order accurate with respect to the time variable and first-order accurate with respect to the space variable. Moreover, the proposed method is uniformly convergent with respect to the perturbation parameter ϵ. To validate the theoretical predictions, numerical experimentation is performed and the numerical results support the theoretical findings.

It is very important to note that previously the method of Crank–Nicolson and NSFDM are developed for only integer order singularly perturbed delay partial differential equations. According to the best of our knowledge, this has not been done so far for the case of TF-SPDPDEs. Hence, this is the first work based on the fractional order singularly perturbed delay partial differential equations. In the future, we will extend this work for higher order numerical schemes to solve TF-SPDPDEs.

PUBLIC INTEREST STATEMENT

Nowadays, due to the non-locality of fractional differential operators, it is confirmed that many anomalous diffusion processes associated with complex systems in science and engineering are described more accurately by fractional partial differential equations than classical partial differential equations. If the small perturbation parameter multiplies the highest order derivative of the time-fractional partial differential equation and when it contains the delay parameter, then the problem is singularly perturbed time-fractional delay partial differential equation. The solution of such types of problems exhibit boundary layer phenomena. Due to the existence of boundary layer in the solution, if we use classical numerical methods, large oscillations occur in the computed solution, leading to numerical instability and inaccurate solution. To overcome this problem, we developed uniformly convergent numerical method with respect to the perturbation parameter.

Supplemental material

sn-sample-bib.tex

Download Latex File (1.1 KB)

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/27684830.2023.2293373.

Additional information

Notes on contributors

Habtamu Getachew Kumie

Awoke Andargie Tirunehi is an Associate Professor of Mathematics and member of Department of Mathematics, College of Natural Sciences at Bahir Dar University, Ethiopia. His research interest is in the areas of Applied Mathematics and Numerical Analysis, including the Theory of Perturbation. So far he has published more than 30 research articles in reputable journals.

Habtamu Getachew Kumie is a PhD candidate at Bahir Dar University, Ethiopia. His research interest is in the areas of Applied Mathematics and Numerical Analysis, including the Numerical Solution of Singularly Perturbed Fractional Differential Equations.

Getachew Adamu Derese

Getachew Adamu Derese is an Associate Professor of Mathematics and member of Department of Mathematics, College of Natural Sciences at Bahir Dar University, Ethiopia. His research interests are in the areas of Applied Mathematics and Numerical Analysis, including the lubrication of bearings. So far he has published more than 20 research articles in reputable journals.

References

Ansari, A. R., Bakr, S. A., & Shishkin, G. I. (2007). A parameter-robust finite difference method for singularly perturbed delay parabolic partial differential equations. Journal of Computational and Applied Mathematics, 205(1), 552–566. https://doi.org/10.1016/j.cam.2006.05.032
Web of Science ®Google Scholar
Ata, E., & Kıymaz, I. O.(2023). New generalized Mellin transform and applications to partial and fractional differential equations. International Journal of Mathematics and Computer in Engineering, 2023.
Google Scholar
Atangana, A. (2014). On the singular perturbations for fractional differential equation. Scientific World Journal, 2014, 1–9. https://doi.org/10.1155/2014/752371
Google Scholar
Bashier, E. B. M., & Patidar, K. C. (2011). A second-order fitted operator finite difference method for a singularly perturbed delay parabolic partial differential equation. Journal of Difference Equations and Applications, 17(5), 779–794. https://doi.org/10.1080/10236190903305450
Web of Science ®Google Scholar
Chen, C.-M., Liu, F., Turner, I., & Anh, V. (2007). A Fourier method for the fractional diffusion equation describing sub-diffusion. Journal of Computational Physics, 227(2), 886–897. https://doi.org/10.1016/j.jcp.2007.05.012
Web of Science ®Google Scholar
Choudhary, R., Singh, S., & Kumar, D. (2022). A second-order numerical scheme for the time-fractional partial differential equations with a time delay. Computational and Applied Mathematics, 41(3), 1–28. https://doi.org/10.1007/s40314-022-01810-9
Google Scholar
Dai, Z., Peng, Y., Mansy, H. A., Sandler, R. H., & Royston, T. J. (2015). A model of lung parenchyma stress relaxation using fractional viscoelasticity. Medical Engineering & Physics, 37(8), 752–758. https://doi.org/10.1016/j.medengphy.2015.05.003
PubMed Web of Science ®Google Scholar
Daraghmeh, A., Qatanani, N., & Saadeh, A. (2020). Numerical solution of fractional differential equations. Applied Mathematics, 11(11), 1100–1115. https://doi.org/10.4236/am.2020.1111074
Google Scholar
Doolan, E. P., Miller, J. J., & Schilders, W. H. (1980). Uniform numerical methods for problems with initial and boundary layers. Boole Press.
Google Scholar
Duressa, G. F., Daba, I. T., & Deressa, C. T. (2023). A systematic review on the solution methodology of singularly perturbed differential difference equations. Mathematics, 11(5), 1108. https://doi.org/10.3390/math11051108
Web of Science ®Google Scholar
Govindarao, L., & Mohapatra, J. (2018). A second order numerical method for singularly perturbed delay parabolic partial differential equation. Engineering Computations, 36(2), 420–444. https://doi.org/10.1108/EC-08-2018-0337
Web of Science ®Google Scholar
Gowrisankar, S., & Natesan, S. (2017). ε-uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations. International Journal of Computer Mathematics, 94(5), 902–921. https://doi.org/10.1080/00207160.2016.1154948
Web of Science ®Google Scholar
Karatay, I., & Bayramoglu, S. R. (2012). An efficient difference scheme for time fractional advection dispersion equations. Applied Mathematical Sciences, 6(98), 4869–4878. https://doi.org/10.1155/2012/548292
Google Scholar
Kellogg, R. B., & Tsan, A. (1978). Analysis of some difference approximations for a singular perturbation problem without turning points. Mathematics of Computation, 32(144), 1025–1039. https://doi.org/10.1090/S0025-5718-1978-0483484-9
Web of Science ®Google Scholar
Kilbas, A. A., Srivastava, H. M., & Trujillo, J. J. (2006). Theory and applications of fractional differential equations. North-Holland204.
Google Scholar
Kumar, D., & Kumari, P. (2019). A parameter-uniform numerical scheme for the parabolic singularly perturbed initial boundary value problems with large time delay. Journal of Applied Mathematics and Computing, 59(1–2), 179–206. https://doi.org/10.1007/s12190-018-1174-z
Web of Science ®Google Scholar
Kumar, D., & Kumari, P. (2020). A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag. Numerical Methods for Partial Differential Equations, 36(4), 868–886. https://doi.org/10.1002/num.22455
Web of Science ®Google Scholar
Kumar, S., Pandey, R. K., Kumar, K., Kamal, S., & Dinh, T. N. (2022). Finite difference–collocation method for the generalized fractional diffusion equation. Fractal and Fractional, 6(7), 387. https://doi.org/10.3390/fractalfract6070387
Web of Science ®Google Scholar
Kumar, K., & Vigo-Aguiar, J. (2021). Numerical solution of time-fractional singularly perturbed convection–diffusion problems with a delay in time. Mathematical Methods in the Applied Sciences, 44(4), 3080–3097. https://doi.org/10.1002/mma.6477
Web of Science ®Google Scholar
Lin, Y., & Xu, C. (2007). Finite difference/spectral approximations for the time-fractional diffusion equation. Journal of Computational Physics, 225(2), 1533–1552. https://doi.org/10.1016/j.jcp.2007.02.001
Web of Science ®Google Scholar
Liu, X., Abbas, M., Yang, H., Qin, X., & Nazir, T. (2021). Novel finite point approach for solving time-fractional convection-dominated diffusion equations. Advances in Difference Equations, 2021(1), 1–22. https://doi.org/10.1186/s13662-020-03178-8
Web of Science ®Google Scholar
Mbroh, N. A., & Munyakazi, J. B. (2019). A fitted operator finite difference method of lines for singularly perturbed parabolic convection–diffusion problems. Mathematics and Computers in Simulation, 165, 156–171. https://doi.org/10.1016/j.matcom.2019.03.007
Web of Science ®Google Scholar
Mickens, R. E. (2005). Advances in the applications of nonstandard finite difference schemes. World Scientific.
Google Scholar
Miller, J. J., O’riordan, E., & Shishkin, G. I. (1996). Fitted numerical methods for singular perturbation problems: Error estimates in the maximum norm for linear problems in one and two dimensions. World Scientific.
Google Scholar
Negero, N., & Duressa, G. (2022). An efficient numerical approach for singularly perturbed parabolic convection-diffusion problems with large time-lag. Journal of Mathematical Modeling, 10(2), 173–110.
Google Scholar
Odibat, Z. M., & Momani, S. (2006). Application of variational iteration method to nonlinear differential equations of fractional order. International Journal of Nonlinear Sciences and Numerical Simulation, 7(1), 27–34. https://doi.org/10.1515/IJNSNS.2006.7.1.27
Web of Science ®Google Scholar
Qiu, W., Chen, H., & Zheng, X. (2019). An implicit difference scheme and algorithm implementation for the one-dimensional time-fractional burgers equations. Mathematics and Computers in Simulation, 166, 298–314. https://doi.org/10.1016/j.matcom.2019.05.017
Web of Science ®Google Scholar
Rihan, F. A. (2010). Computational methods for delay parabolic and time-fractional partial differential equations. Numerical Methods for Partial Differential Equations, 26(6), 1556–1571. https://doi.org/10.1002/num.20504
Web of Science ®Google Scholar
Sun, H., Zhang, Y., Baleanu, D., Chen, W., & Chen, Y. (2018). A new collection of real world applications of fractional calculus in science and engineering. Communications in Nonlinear Science and Numerical Simulation, 64, 213–231. https://doi.org/10.1016/j.cnsns.2018.04.019
Web of Science ®Google Scholar
Tirunehi, A. A., Derese, G. A., & Ayele, M. A. (2023). Singularly perturbed reaction diffusion problem with large spatial delay via non-standard fitted operator method. Research in Mathematics, 10(1). https://doi.org/10.1080/27684830.2023.2171698
Google Scholar
Turuna, D. A., Woldaregay, M. M., & Duressa, G. F. (2020). Uniformly convergent numerical method for singularly perturbed convection-diffusion problems. Kyungpook Mathematical Journal, 60(3), 629–645.
Web of Science ®Google Scholar

Crank–Nicolson method for solving time-fractional singularly perturbed delay partial differential equations

ABSTRACT

1. Introduction

1.1. Properties of the continuous problem

2. Numerical scheme

2.1. Time discretization

2.2. Spatial discretization

3. Convergence analysis

4. Numerical examples

Table 1. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using fractional order β = 0.8 with different values of $ϵ .$

Table 2. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using $ϵ = 10^{- 10}$ with different values of β

Table 3. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using fractional order β = 0.8 with different values of ϵ

Table 4. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using $ϵ = 10^{- 10}$ with different values of β

5. Conclusion

PUBLIC INTEREST STATEMENT

sn-sample-bib.tex

Disclosure statement

Supplementary material

Notes on contributors

Habtamu Getachew Kumie

Getachew Adamu Derese

References

Information for

Open access

Opportunities

Help and information

Crank–Nicolson method for solving time-fractional singularly perturbed delay partial differential equations

ABSTRACT

1. Introduction

1.1. Properties of the continuous problem

2. Numerical scheme

2.1. Time discretization

2.2. Spatial discretization

3. Convergence analysis

4. Numerical examples

Table 1. Maximum pointwise error (EM,N) and rate of convergence (RM,N) for example 4.1 using fractional order β = 0.8 with different values of ϵ.

Table 2. Maximum pointwise error (EM,N) and rate of convergence (RM,N) for example 4.1 using ϵ=10−10 with different values of β

Table 3. Maximum pointwise error (EM,N) and rate of convergence (RM,N) for example 4.2 using fractional order β = 0.8 with different values of ϵ

Table 4. Maximum pointwise error (EM,N) and rate of convergence (RM,N) for example 4.2 using ϵ=10−10 with different values of β

5. Conclusion

PUBLIC INTEREST STATEMENT

sn-sample-bib.tex

Disclosure statement

Supplementary material

Additional information

Notes on contributors

Habtamu Getachew Kumie

Getachew Adamu Derese

References

Related research

To cite this article:

Download citation

Your download is now in progress and you may close this window

Login or register to access this feature

Information for

Open access

Opportunities

Help and information

Keep up to date

Table 1. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using fractional order β = 0.8 with different values of $ϵ .$

Table 2. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.1 using $ϵ = 10^{- 10}$ with different values of β

Table 3. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using fractional order β = 0.8 with different values of ϵ

Table 4. Maximum pointwise error $(E^{M, N})$ and rate of convergence $(R^{M, N})$ for example 4.2 using $ϵ = 10^{- 10}$ with different values of β