Layer adapted finite element method for two-parameter singularly perturbed delay parabolic problems

ABSTRACT

A layer-adapted finite element method is applied to solve a singularly perturbed delay parabolic differential equation with two perturbation parameters. The numerical method involves two steps of discretization in the decoupled approach. First, the temporal variable is discretized uniformly and approximated by a finite element piecewise linear continuous function. The time derivative of this approximation in each mesh is constant in the time direction. It is a function of the space variable. Then, the time derivative term in the equation is approximated by the equivalent Crank-Nicholson finite difference method. Next, for the semi-discretized problem, a fitted mesh (layer-adapted) finite element method is applied using the Shishkin mesh for the spatial discretization. The solution in the spatial elements is approximated by piecewise quadratic polynomial functions. Finally, the finite element solution is presented as a linear combination of the piecewise linear temporal FE approximate functions and the piecewise quadratic spatial FE approximate functions. That is, a piece wise smooth continuous function in the given domain, which empowers finding an approximation solution at any desired point of the domain. It is proven that the numerical method is convergent with the second order of convergence in the temporal direction and with almost the second order of convergence in the spatial direction. Using the MATLAB software, some numerical experiments are performed to verify the applicability of the numerical method.

KEYWORDS:

• Two-parameter singularly perturbed problem
• delay differential equation
• Shishkin Mesh
• finite element
• fitted mesh method

1. Introduction

Two-parameter singularly perturbed delay parabolic differential equations refer to a class of mathematical equations that involve partial derivatives, delay in the time variable, and underlying parameters that multiply the two higher spatial derivative terms that characterize the degree of perturbation. The solution to such problems possesses a highly steep or descent gradient, especially when the perturbation parameters are very small. Standard finite difference methods (FDM) and finite element methods (FEM) lead to unsuitable approximations in these areas unless an unacceptably large number of meshes are used, which is not practical due to the limited memory capacity of computers and their exposure to rounding errors (M. K. Kadalbajoo & Gupta, 2010). A class of these types of equations, which involve two perturbation parameters and a delay term, arise in different science and engineering fields.

Two-parameter singularly perturbed differential equations have an important role in the application of fluid flow through unsaturated porous media such as soils and industrial filters (Bhathawala & Verma, 1975), computational fluid dynamics of incomprehensible Navier-Stokes equation (Heywood & Rannacher, 1990), and oscillating drop in an electric field (Feng, 1990). Professor Robert E. O’Malley is the first researcher to study two-parameter perturbation problems. In his 1967 paper, he considers the boundary value problems of the form $\mathit{ϵ}{y}^{″}+\mathit{\mu a}(\mathit{x})y^{\prime }+\mathit{b}(\mathit{x})\mathit{y}=0$ for $0\le \mathit{x}\le 1,$ where ϵ and µ are small positive interrelated parameters (O’Malley, 1967). The behavior of the asymptotic solution to the problem is well studied for different parities of the coefficient functions a(x) and b(x) in different prescribed cases of initial and boundary conditions. It is also mentioned that the existence of the solution and appearance of boundary layers are mostly based on the condition whether ${\displaystyle \frac{ϵ}{{\mu }^2}\to 0}$ as $\mathit{\mu }\to 0$ or ${\displaystyle \frac{{\mu }^2}{\mathit{ϵ}}\to 0}$ as $\mathit{ϵ}\to 0$ (Gao et al., 2023; M. K. Kadalbajoo & Gupta, 2010; O’Malley, 1967). The findings of his studies serve as a foundation for further and extended literary works like two-dimensional partial differential equations and parabolic problems. Many studies have been done on the topic of two-parameter singularly perturbed parabolic differential equations (Bullo et al., 2019; P. Das & Mehrmann, 2016; A. Das & Natesan, 2015; Jha & Kadalbajoo, 2015; M. Kadalbajoo & Yadaw, 2012; Munyakazi, 2015; O’Riordan et al., 2006) that applied different numerical approaches to solve the problem. But all the mentioned papers did not consider the delay case.

There are mathematical models of singularly perturbed differential equations that depend not only on the instant but also on the past state. The delay term represents the effect of the past state of the solution on the current state. This delay term makes the equation more complex but realistic. Time-delay parabolic problems with a perturbation parameter are applicable in different areas of science and engineering. For example, a respiratory control system in the case of unstable patterns of ventilation, monitoring the temperature distribution in a metal sheet, population dynamics, and chemical reactions or pollution in fluids are some areas of application that characterize the delay term in the model of the differential equation. These kinds of problems are addressed by some literature (Bashier & Patidar, 2011; Mbroh et al., 2020; Salama & Al-Amery, 2017; Woldaregay et al., 2021).

According to our recent review, singularly perturbed parabolic differential equations that involve two perturbation parameters and a time delay have been studied by several authors. Some applied the fitted mesh (layer-adapted) finite difference method, and others used the fitted operator finite difference method to develop a uniformly convergent scheme independent of the parameters. Govindarao et al. (2019) solved the problem by applying implicit Euler finite difference approximation on uniform mesh in the direction of time and upwind finite difference approximation on Shishkin and Bakhvalov-Shishkin mesh as well in the space direction. They proved that their scheme is parameter-uniform with $\mathcal{O}(\mathit{h}+\mathit{k})$, where h and k are spatial and temporal mesh sizes, respectively. Kumar et al. (2020) improved the spatial order of convergence to almost second order by implementing a hybrid monotone finite difference scheme in the direction of space. Singh et al. (2022) applied the Crank-Nicolson scheme in temporal direction and the B-spline collocation method using piecewise uniform mesh for spatial direction. Ayele et al. (2023) found a numerical solution to the problem by using a combination of the fitted operator scheme of a cubic spline with θ - method on a uniform mesh grid. The other two articles, Negero et al. (2023) and Mohye et al. (2023), also applied the nonstandard finite difference method. Both articles followed the Crank-Nicolson approach for temporal discretization and the nonstandard finite difference scheme for spatial discretization on uniform mesh grids. Worku et al. (2023) considered this problem when the diffusion term is dominant, using a finite difference scheme to obtain almost a second order accuracy.

The finite element method (FEM) is more advantageous than FDM in terms of its ability to provide more precise and continuous approximations in the given domain, especially when the domain is irregular (Li et al., 2017). The finite element model of such time-dependent problems can be developed in two ways: the coupled formulation, in which the time t is treated as an additional coordinate along with the spatial coordinate x, and the decoupled formulation, where the time and spatial variations are assumed to be separable (Reddy, 2019). This study aims to develop a parameter uniform convergent scheme to solve a two-parameter singularly perturbed delay parabolic differential equation numerically by applying a layer-adapted Galerkin finite element method (GFEM). The decoupled formulation approach will be employed, which is a two-step process. To control the perturbation effect in the boundary layers, a piecewise uniform (Shishkin) mesh discretization is used in the spatial direction, followed by a uniform discretization along the temporal direction.

This article consists of the following structure: the existence and uniqueness of the analytical solution, and necessary bounds for its derivatives are discussed in Section 2. In Section 3, the discretization of the domain and the development of the numerical scheme derived. In Section 4, the convergence of the finite element solution is proved. Section 5 presents theoretical results using numerical solutions of counter-examples with tables and graphs. Finally, a brief conclusion is given in Section 6.

Notations: M and N represents the number of meshes in time and space direction, respectively, C denotes the generic positive constant independent of perturbation parameters and mesh sizes. For a given continuous functions f(x) and g(x) on $\overline{D}$, $||\mathit{f}||$ denotes the supreme(maximum) norm, ${\displaystyle \underset{(x)\in \overline{D}}{max}\{|\mathit{f}(\mathit{x})|\}}$ or ${\displaystyle \underset{1\le \mathit{i}\le \mathit{N}}{max}|\mathit{f}(x_i)|}$ if f is discrete function on the discrete domain ${\displaystyle {\overline{D}}^{\mathit{N}}}$.

The inner product operator $⟨.⟩$ is given by

$\begin{array}{l}⟨\mathit{f}(\mathit{x}),\mathit{g}(\mathit{x})⟩={\int }_{\overline{D}}\mathit{f}(\mathit{x}).\mathit{g}(\mathit{x})\mathit{dx}.\end{array}$

We consider two-parameters singularly perturbed large time-delay one-dimensional parabolic problem defined on the rectangular domain,

$\begin{array}{l}\mathcal{D}={\mathcal{D}}_x\times {\mathcal{D}}_t,\mathit{for}{\mathcal{D}}_x=(0,1),{\mathcal{D}}_t=(0,\mathit{T}]\end{array}$

given by

(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1)

with initial and boundary conditions

(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2)

where ɛ and µ are perturbation parameters such that $0<\mathit{\epsilon }\ll 1$ and $0<\mathit{\mu }\ll 1$, considering the large time-delay parameter τ, such that the terminal time, $\mathit{T}=\mathit{k\tau }$ for some positive integer k. The functions $\mathit{a}(\mathit{x},\mathit{t}),\mathit{b}(\mathit{x},\mathit{t}),\mathit{c}(\mathit{x},\mathit{t})$, and $\mathit{f}(\mathit{x},\mathit{t}),$ are assumed to be sufficiently smooth and bounded; those satisfy

$\begin{array}{l}\mathit{a}(\mathit{x},\mathit{t})\ge \mathit{\alpha }>0,\mathit{b}(\mathit{x},\mathit{t})\ge \mathit{\beta }>0\mathit{and}\\ \mathit{c}(\mathit{x},\mathit{t})\ge \mathit{\vartheta }>0\mathit{on}\overline{\mathcal{D}}=[0,1]\times [0,\mathit{T}].\end{array}$

Let $\mathrm{\Gamma }={\mathrm{\Gamma }}_b\cup {\mathrm{\Gamma }}_l\cup {\mathrm{\Gamma }}_r$ and define the differential operator, ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}$ by

(3) ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}\equiv \frac{\partial }{\mathit{\partial t}}-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}+\mathit{b}(\mathit{x},\mathit{t}).$(3)

2. Behavior of the analytical solution

For all $(\mathit{x},\mathit{t})\in [0,1]\times [0,\mathit{\tau }],$ the delay term $\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })$ is given by the initial condition. Then, the problem is treated as a non-delay two-parameter singularly perturbed parabolic problem (P. Das & Mehrmann, 2016; M. Kadalbajoo & Yadaw, 2012). The computation extends for $(\mathit{x},\mathit{t})\in [0,1]\times [\mathit{\tau },2\mathit{\tau }]$ after solving $\mathit{u}(\mathit{x},\mathit{t})$ in $[0,1]\times [0,\mathit{\tau }]$. Repeating the method step by step, completes the solution for the whole domain $[0,1]\times [0,\mathit{T}]$.

Problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) )-(2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) become a kind of reaction-diffusion when µ = 0 and convection-diffusion when µ = 1. These problems gives important clues about the properties of the solutions in our concern where $0<\mathit{\epsilon },\mathit{\mu }\ll 1$ and both $\mathit{\epsilon },\mathit{\mu }\to 0$. The analysis for two parameters problems mainly started by O’Malley, based on asymptotic expansion methods (O’Malley, 1967). His study notated two basic cases, when ${\mu }^2\le \mathit{C\epsilon }$ and ${\mu }^2\ge \mathit{C\epsilon }$. In the first case, the problem resembles properties of the case when µ = 0 and hence $\mathcal{O}(\sqrt{\mathit{\epsilon }})$ width layers appear in the neighborhood of x = 0 and x = 1. For the second case different layer widths $\mathcal{O}(\frac{\epsilon }{\mu })$ and $\mathcal{O}(\mathit{\mu })$ appear in the neighborhood of x = 0 and x = 1 respectively.

The existence and uniqueness of problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) )-(2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) is established assuming that the data is Holder’s continuous and imposing appropriate compatibility conditions at its boundary (Li et al., 2017; Woldaregay et al., 2021). Therefore, the following hold at the corner points $(0,0)$, $(1,0),(0,-\mathit{\tau })$ and $(1,-\mathit{\tau })$.

(4) $\begin{array}{ll}{\gamma }_b(0,0)& =0,\end{array}$(4)

(5) $\begin{array}{ll}{\gamma }_b(1,0)& =0,\end{array}$(5)

(6) ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}(\mathit{u}(0,0))=\mathit{c}(0,0)\mathit{u}(0,-\mathit{\tau })+\mathit{f}(0,0)\text{and}$(6)

(7) $\begin{array}{l}{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}(\mathit{u}(1,0))=\mathit{c}(1,0)\mathit{u}(1,-\mathit{\tau })+\mathit{f}(1,0).\end{array}$(7)

In addition to that, the continuity and boundedness of the solution and its derivatives have to be certain. To show this, it is necessary to first validate that the differential operator, ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}$ holds the maximum principle.

Lemma.

2.1. Continuous maximum principle:

Suppose $\mathit{\chi }(\mathit{x},\mathit{t})\in {\mathcal{C}}^{2,1}(\mathcal{D})\cap {\mathcal{C}}^{0,0}(\overline{\mathcal{D}})$, such that $\mathit{\chi }(\mathit{x},\mathit{t})\ge 0,$ $\mathrm{\forall }(\mathit{x},\mathit{t})\in \mathrm{\Gamma }$ and ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}(\mathit{\chi }(\mathit{x},\mathit{t}))\ge 0,$ $\mathrm{\forall }(\mathit{x},\mathit{t})\in \mathcal{D}$, then $\mathit{\chi }(\mathit{x},\mathit{t})\ge 0,$ $\mathrm{\forall }(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}$.

Proof.

Suppose $(x^{\prime },t^{\prime })\in \overline{\mathcal{D}}$, ${\displaystyle \mathit{\chi }(x^{\prime },t^{\prime })=\underset{(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}}{min}\{\mathit{\chi }(\mathit{x},\mathit{t})\}}$ and $\mathit{\chi }(x^{\prime },t^{\prime })<0.$

By the hypothesis, $\mathit{\chi }(\mathit{x},\mathit{t})\ge 0,$ $\mathrm{\forall }(\mathit{x},\mathit{t})\in \mathrm{\Gamma },$ which implies $(x^{\prime },t^{\prime })\notin \mathrm{\Gamma }.$ Since the minimum is at $(x^{\prime },t^{\prime }),$ ${\chi }_x(x^{\prime },t^{\prime })={\chi }_t(x^{\prime },t^{\prime })=0$ and ${\chi }_{\mathit{xx}}(x^{\prime },t^{\prime })\ge 0.$

$\begin{array}{ll}{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}(\mathit{\chi }(x^{\prime },t^{\prime }))& ={\chi }_t(x^{\prime },t^{\prime })-\mathit{\epsilon }{\chi }_{\mathit{xx}}(x^{\prime },t^{\prime })\\ -\mathit{\mu a}(x^{\prime }){\chi }_x(x^{\prime },t^{\prime })+\mathit{b}(x^{\prime },t^{\prime })\mathit{\chi }(x^{\prime },t^{\prime })<0.\end{array}$

This contradicts the given hypothesis. Therefore, $\mathit{\chi }(\mathit{x},\mathit{t})\ge 0,\mathrm{\forall }(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}$.

Lemma 2.2.

Let $\mathit{u}(\mathit{x},\mathit{t})$ be the solution of problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) ) - (2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ), then $\mathit{u}(\mathit{x},\mathit{t})$ satisfies the following bounds.

(8) $|\mathit{u}(\mathit{x},\mathit{t})-{\gamma }_b(\mathit{x},0)|\le \mathit{Ct},$(8)

(9) $\Vert \mathit{u}(\mathit{x},\mathit{t})\Vert \le \mathit{C},$(9)

for some constant C independent of ɛ and µ.

Proof.

Set two functions s and q defined as

(10) $\mathit{s}(\mathit{x},\mathit{t})=\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})-{\gamma }_b(\mathit{x},0),0\le \mathit{t}\le \mathit{T},\\ \mathit{u}(\mathit{x},\mathit{t})-{\gamma }_b(\mathit{x},\mathit{t}),-\mathit{\tau }\le \mathit{t}\le 0,\end{array}$(10)

and

(11) $\mathit{q}(\mathit{x},\mathit{t})=\{\begin{array}{ll}\mathit{Ct},& 0\le \mathit{t}\le \mathit{T},\\ 0.& -\mathit{\tau }\le \mathit{t}\le 0.\end{array}$(11)

Applying the differential operator, ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }},$

$\begin{array}{l}{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}[\mathit{s}(\mathit{x},\mathit{t})]=\{\begin{array}{l}\mathit{c}(\mathit{x},\mathit{t}){\gamma }_b(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})+\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{\gamma }_b(\mathit{x},0)\\ +\mathit{\mu a}\frac{d}{\mathit{dx}}{\gamma }_b(\mathit{x},0)-\mathit{b}{\gamma }_b(\mathit{x},0),& 0\le \mathit{t}\le \mathit{T},\\ 0,& -\mathit{\tau }\le \mathit{t}\le 0.\end{array}\\ {\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}[\mathit{q}(\mathit{x},\mathit{t})]=\{\begin{array}{ll}\mathit{C}(\mathit{bt}+1),& 0\le \mathit{t}\le \mathit{T},\\ 0.& -\mathit{\tau }\le \mathit{t}\le 0.\end{array}\end{array}$

For sufficiently large C, we have

$\begin{array}{l}|{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}[\mathit{s}(\mathit{x},\mathit{t})]|\le {\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}[\mathit{q}(\mathit{x},\mathit{t})].\end{array}$

Since ${\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}$ satisfies the maximum principle, it holds the relation

$\begin{array}{l}|\mathit{s}(\mathit{x},\mathit{t})|\le \mathit{q}(\mathit{x},\mathit{t}).\text{ That is}|\mathit{u}(\mathit{x},\mathit{t})-{\gamma }_b(\mathit{x},0)|\le \mathit{Ct}\end{array}$

The comparisons $|\mathit{u}(\mathit{x},\mathit{t})|-|{\gamma }_b(\mathit{x},0)|\le |\mathit{u}(\mathit{x},\mathit{t})-{\gamma }_b(\mathit{x},0)|\le \mathit{Ct}$ and $\mathit{t}\le \mathit{T},\mathrm{\forall t}\in [0,\mathit{T}]$ imply $|\mathit{u}(\mathit{x},\mathit{t})|\le \mathit{Ct}+|{\gamma }_b(\mathit{x},0)|\le \mathit{CT}+|{\gamma }_b(\mathit{x},0)|.$

Then, taking the supreme of all $(\mathit{x},\mathit{t})\in \overline{\mathcal{D}},$ gives $\Vert \mathit{u}(\mathit{x},\mathit{t})\Vert \le \mathit{C}.$

Lemma 2.3.

For non-negative integers i and j, such that $0\le \mathit{i}+2\mathit{j}\le 4$ the derivatives of solution of problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) )-(2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) satisfies the following bounds

(12) $\begin{array}{ll}\Vert \frac{{\mathit{\partial }}^{\mathit{i}+\mathit{j}}\mathit{u}}{\mathit{\partial }{x}^i\mathit{\partial }{t}^j}\Vert & \le \{\begin{array}{ll}\mathit{C}{\epsilon }^{-\frac{1}{2}\mathit{i}},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}(\frac{\mu }{\epsilon }{)}^i(\frac{{\mu }^2}{\mathit{\epsilon }}{)}^j,& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon },\end{array}\text{where}\\ \mathit{\varsigma }& =\underset{(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}}{min}\{\frac{\mathit{b}(\mathit{x},\mathit{t})}{\mathit{a}(\mathit{x},\mathit{t})}\}.\end{array}$(12)

Proof.

Refer (O’Riordan et al., 2006)

The bounds of the regular and singular component solutions and their derivatives are necessary for later error estimation. The decompositions of the solution $\mathit{u}(\mathit{x},\mathit{t})$ of the problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) ) - (2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) are the regular component $\mathit{v}(\mathit{x},\mathit{t})$, which characterizes the solution behavior outside the boundary layers, and a singular component $\mathit{w}(\mathit{x},\mathit{t})$, which characterizes the solution behavior inside the boundary layers. Further, we decompose the singular component into the sum of the left and right singular parts, $\mathit{w}(\mathit{x},\mathit{t})=w_L(\mathit{x},\mathit{t})+w_R(\mathit{x},\mathit{t})$, so that $\mathit{u}(\mathit{x},\mathit{t})=\mathit{v}(\mathit{x},\mathit{t})+w_L(\mathit{x},\mathit{t})+w_R(\mathit{x},\mathit{t})$. Each components satisfy

(13) $\{\begin{array}{l}\left({\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}\right)(\mathit{v})=\mathit{c}(\mathit{x},\mathit{t})\mathit{v}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t}),\\ \mathit{v}(\mathit{x},\mathit{t})=\mathit{u}(\mathit{x},\mathit{t}),\text{on }{\mathrm{\Gamma }}_b,\\ \mathit{v}(0,\mathit{t})={\phi }_{\mathit{ol}}(\mathit{t}),\text{(conveniently chosen),}\\ \mathit{v}(1,\mathit{t})={\phi }_{\mathit{or}}(\mathit{t}),\text{(conveniently chosen),}\end{array}$(13)

(14) $\{\begin{array}{l}\left({\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}\right)(w_L)=\mathit{c}(\mathit{x},\mathit{t})w_L(\mathit{x},\mathit{t}-\mathit{\tau }),\\ w_L(\mathit{x},\mathit{t})=0,\text{on }{\mathrm{\Gamma }}_b,\\ w_L(1,\mathit{t})=0,\\ w_L(0,\mathit{t})=\mathit{u}(0,\mathit{t})-\mathit{v}(0,\mathit{t}),\end{array}$(14)

(15) $\{\begin{array}{l}\left({\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}\right)(w_R)=\mathit{c}(\mathit{x},\mathit{t})w_R(\mathit{x},\mathit{t}-\mathit{\tau }),\\ w_R(\mathit{x},\mathit{t})=0,\text{on }{\mathrm{\Gamma }}_b,\\ w_R(0,\mathit{t})=0,\\ w_R(1,\mathit{t})=\mathit{u}(1,\mathit{t})-\mathit{v}(1,\mathit{t}).\end{array}$(15)

Lemma 2.4.

The bounds of singular components have given by

(16) $\begin{array}{ll}|w_L|& \le \{\begin{array}{l}\mathit{Cexp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}\mathit{x}),\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{Cexp}(-\frac{\mathit{\alpha \mu }}{\epsilon }\mathit{x}),\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon },\end{array}\text{and }\\ |w_R|& \le \{\begin{array}{l}\mathit{Cexp}(-\frac{\sqrt{\alpha \varsigma }}{2\sqrt{\mathit{\epsilon }}}(1-\mathit{x})),\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{Cexp}(-\frac{\varsigma }{2\mathit{\mu }}(1-\mathit{x})),\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}\end{array}$(16)

Proof.

Refer (Kumar et al., 2020; O’Riordan et al., 2006)

Lemma 2.5.

The derivative of regular and singular components of $\mathit{u}(\mathit{x},\mathit{t})$ satisfy the bounds.

i. If $\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }$,then

(17) $\begin{array}{l}\Vert \frac{{\mathit{\partial }}^{\mathit{i}+\mathit{j}}\mathit{v}}{\mathit{\partial }{x}^i\mathit{\partial }{t}^j}\Vert \le \mathit{C}(1+(\frac{\epsilon }{\mu }{)}^{3-\mathit{i}}(\frac{{\mu }^2}{\mathit{\epsilon }}{)}^{\mathit{j}}),\mathit{for}0\le \mathit{i}+2\mathit{j}\le 4,\end{array}$(17)

(18) $\begin{array}{ll}\Vert \frac{{\mathit{\partial }}^{\mathit{i}}{w}_L}{\mathit{\partial }{x}^i}\Vert & \le \mathit{C}{\left(\frac{\mu }{\epsilon }\right)}^{\mathit{i}},\mathit{for}0\le \mathit{i}\le 4\mathit{and}\\ \Vert \frac{{\mathit{\partial }}^2{w}_L}{\mathit{\partial }{t}^2}\Vert & \le \mathit{C}(1+{\mu }^2{\epsilon }^{-1}),\end{array}$(18)

(19) $\begin{array}{ll}\Vert \frac{{\mathit{\partial }}^{\mathit{i}}{w}_R}{\mathit{\partial }{x}^i}\Vert & \le \mathit{C}({\mu }^{-\mathit{i}}+{\mu }^{-1}{\epsilon }^{2-\mathit{i}}),\mathit{for}0\le \mathit{i}\le 4\mathit{and}\\ \Vert \frac{{\mathit{\partial }}^{\mathit{j}}{w}_R}{\mathit{\partial }{t}^j}\Vert & \le \mathit{Cforj}=1,2.\end{array}$(19)

ii. If $\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon }$, then

(20) $\begin{array}{l}\Vert \frac{{\mathit{\partial }}^{\mathit{i}+\mathit{j}}\mathit{v}}{\mathit{\partial }{x}^i\mathit{\partial }{t}^j}\Vert \le \mathit{C},\mathit{for}0\le \mathit{i}+2\mathit{j}\le 4,\end{array}$(20)

w_L and w_R satisfy the derivative bounds of Lemma 2.3.

Proof.

Refer (Kumar et al., 2020; O’Riordan et al., 2006)

3. Formulation of numerical scheme

We use two steps in the procedure to generate the finite element model. At the beginning of the construction of the semi-discrete finite element model, there is no integration by part for the time variable. Using the Galerkin finite element method in combination with the uniform discretization of time Crank-Nicolson approach, the numerical scheme is formulated on a piecewise uniform Shishkin mesh type in the spatial direction.

3.1. Time discretization

Let us discretize the temporal variable interval, $[-\mathit{\tau },\mathit{T}]$ into m + M uniform meshes, where m and M are the number of meshes in the time-lag interval $[-\mathit{\tau },0]$ and the domain time component, $[0,\mathit{T}]$ respectively.

The time step $t_{\mathit{J}-\mathit{m}}=(\mathit{J}-\mathit{m}-1)\mathrm{\Delta t},$ for $\mathit{J}=1,..,\mathit{m}+\mathit{M}+1$, where $\mathrm{\Delta t}$ is the uniform time mesh size, in which $\mathrm{\Delta t}=\mathit{\tau }/\mathit{m}=\mathit{T}/\mathit{M}$ as we impose $\mathit{T}=\mathit{k\tau }$ for some positive integer k. The discrete temporal domain is given by

$\begin{array}{l}{\mathcal{D}}_t^{\mathit{\text{M}}}=\{t_{\mathit{j}}=(\mathit{j}-1)\mathrm{\Delta t},\mathit{j}=1,\dots ,\mathit{M}+1\}.\end{array}$

The FE solutions are piece-wise smooth, they may not differentiable at the mesh points. So we need a function space of L_p space whose generalized derivatives also belongs to L_p, the so called Sobolev space. The Sobolev space, $W_2^1$ is the set of functions from L₂ whose first generalized derivative exist and belongs to $L_2.$ It is standard to use the Hilbert space $H^1=W_2^1.$ $H_0^1$ is the subspace of H¹ in which the functions vanish at the boundary (Grossmann et al., 2007).

Define linear base functions ${\displaystyle \{{\psi }_j(\mathit{t}){\}}_1^M}$ with respect to variable t, by

(21) $\{{\psi }_j(\mathit{t}){\}}_1^M=\{\begin{array}{ll}\frac{\mathit{t}-{\mathit{t}}_{\mathit{j}-1}}{\mathrm{\Delta t}},& t_{\mathit{j}-1}\le \mathit{t}\le t_j,\\ \frac{{\mathit{t}}_{\mathit{j}+1}-\mathit{t}}{\mathrm{\Delta t}},& t_j\le \mathit{t}\le t_{\mathit{j}+1},\\ 0,& \mathit{otherwise}.\end{array}$(21)

Let ${\displaystyle V_1({\overline{\mathcal{D}}}_{\mathit{t}}^{\mathit{\text{M}}})=\mathit{Span}\{{\psi }_j(\mathit{t}){\}}_1^M}$ and set a Sobolev space $\mathit{V}=H_0^1({\overline{\mathcal{D}}}_{\mathit{t}}^{\mathit{\text{M}}}).$ clearly, at each x the piecewise linear approximation $\mathit{U}(\mathit{t})\in V_1({\overline{\mathcal{D}}}_{\mathit{t}}^{\mathit{\text{M}}})\subset \mathit{V}$.

The solution $\mathit{u}(\mathit{x},\mathit{t})$ is approximated as a linear combination of $U_j(\mathit{x})$ and ${\psi }_j(\mathit{t})$.

(22) $\mathit{u}(\mathit{x},\mathit{t})\approx \mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}}{U}_j(\mathit{x}){\psi }_j(\mathit{t}),$(22)

where the function $U_j(\mathit{x})$ traces the function $\mathit{u}(\mathit{x},t_j)$ at $\mathit{t}=t_j,$ which is going to be determined in a later step (Reddy, 2019).

The first-order time derivative of the approximation in 22(22) $\mathit{u}(\mathit{x},\mathit{t})\approx \mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}}{U}_j(\mathit{x}){\psi }_j(\mathit{t}),$(22) is independent of time. That is

(23) $\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})\approx \sum _{\mathit{j}=1}^{\mathit{M}}\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}.$(23)

At a fixed time interval $[t_{\mathit{j}-1},t_j],$ for some $\mathit{j}=1,2,\dots ,\mathit{M},$ 23(23) $\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})\approx \sum _{\mathit{j}=1}^{\mathit{M}}\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}.$(23) becomes

(24) $\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})\approx \frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}.$(24)

So, we can write problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) ) - (2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) at each j^th time level as,

(25) $\begin{array}{ll}\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}& -\mathit{\epsilon }\frac{{\partial }^2{U}_j(\mathit{x})}{\mathit{\partial }{x}^2}-\mathit{\mu }{a}_j(\mathit{x})\frac{\partial U_j(\mathit{x})}{\mathit{\partial x}}\\ +b_j(\mathit{x})U_j(\mathit{x})=c_j(\mathit{x})U_{\mathit{j}-\mathit{m}}(\mathit{x})+f_j(\mathit{x}).\end{array}$(25)

In our procedure, we used the Crank-Nicolson finite difference approximation of ${\displaystyle \frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},t_{j-\frac{1}{2}})}$, at the mean point of $(\mathit{x},t_{\mathit{j}})$ and $(\mathit{x},t_{\mathit{j}-1})$.

From the Taylor’s series expansion of $\mathit{u}(\mathit{x},t_{\mathit{j}})$ and $\mathit{u}(\mathit{x},t_{\mathit{j}-1})$ about $(\mathit{x},t_{j-\frac{1}{2}})$,

$\begin{array}{l}{U}_j(\mathit{x})=U_{j-\frac{1}{2}}(\mathit{x})+\frac{\mathrm{\Delta t}}{2}\frac{\partial }{\mathit{\partial t}}{U}_{j-\frac{1}{2}}(\mathit{x})+\frac{1}{2}{\left(\frac{\mathrm{\Delta t}}{2}\right)}^2\frac{{\partial }^2}{\mathit{\partial }{t}^2}{U}_{j-\frac{1}{2}}(\mathit{x})+...\end{array}$

and

$\begin{array}{l}{U}_{\mathit{j}-1}(\mathit{x})=U_{j-\frac{1}{2}}(\mathit{x})-\frac{\mathrm{\Delta t}}{2}\frac{\partial }{\mathit{\partial t}}{U}_{j-\frac{1}{2}}(\mathit{x})+\frac{1}{2}{\left(\frac{\mathrm{\Delta t}}{2}\right)}^2\frac{{\partial }^2}{\mathit{\partial }{t}^2}{U}_{j-\frac{1}{2}}(\mathit{x})+...,\end{array}$

we can derive ,

(26) $\frac{\partial }{\mathit{\partial t}}{U}_{j-\frac{1}{2}}(\mathit{x})=\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}+\mathcal{O}((\mathrm{\Delta }\mathit{t}{)}^2)$(26)

and

(27) $U_{j-\frac{1}{2}}(\mathit{x})=\frac{U_j(\mathit{x})+U_{\mathit{j}-1}(\mathit{x})}{2}+\mathcal{O}((\mathrm{\Delta }\mathit{t}{)}^2).$(27)

The governing problem (1) at this mean point, $(\mathit{x},t_{j-\frac{1}{2}})$ is

(28) $\begin{array}{l}(-\mathit{\epsilon }\frac{{\partial }^2\mathit{u}}{\mathit{\partial }{x}^2}-\mathit{\mu a}\frac{\mathit{\partial u}}{\mathit{\partial x}}+\mathit{bu}+\frac{\mathit{\partial u}}{\mathit{\partial t}})(\mathit{x},t_{j-\frac{1}{2}})\\ =-\mathit{c}(\mathit{x},t_{j-\frac{1}{2}})\mathit{u}(\mathit{x},t_{j-\frac{1}{2}-\mathit{m}})+\mathit{f}(\mathit{x},t_{j-\frac{1}{2}}).\end{array}$(28)

Substituting 26(26) $\frac{\partial }{\mathit{\partial t}}{U}_{j-\frac{1}{2}}(\mathit{x})=\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}+\mathcal{O}((\mathrm{\Delta }\mathit{t}{)}^2)$(26) for time derivative and average values between $(\mathit{x},t_j)$ and $(\mathit{x},t_{\mathit{j}-1})$ for spatial derivatives in 28(28) $\begin{array}{l}(-\mathit{\epsilon }\frac{{\partial }^2\mathit{u}}{\mathit{\partial }{x}^2}-\mathit{\mu a}\frac{\mathit{\partial u}}{\mathit{\partial x}}+\mathit{bu}+\frac{\mathit{\partial u}}{\mathit{\partial t}})(\mathit{x},t_{j-\frac{1}{2}})\\ =-\mathit{c}(\mathit{x},t_{j-\frac{1}{2}})\mathit{u}(\mathit{x},t_{j-\frac{1}{2}-\mathit{m}})+\mathit{f}(\mathit{x},t_{j-\frac{1}{2}}).\end{array}$(28) gives

(29) $\begin{array}{l}(-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}+b_{j-\frac{1}{2}}(\mathit{x}))\left(\frac{U_j(\mathit{x})+U_{\mathit{j}-1}(\mathit{x})}{2}\right)\\ +\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}=c_{j-\frac{1}{2}}(\mathit{x})\frac{{\mathit{U}}_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x})}{2}+f_{j-\frac{1}{2}}(\mathit{x}).\end{array}$(29)

29(29) $\begin{array}{l}(-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}+b_{j-\frac{1}{2}}(\mathit{x}))\left(\frac{U_j(\mathit{x})+U_{\mathit{j}-1}(\mathit{x})}{2}\right)\\ +\frac{U_j(\mathit{x})-U_{\mathit{j}-1}(\mathit{x})}{\mathrm{\Delta t}}=c_{j-\frac{1}{2}}(\mathit{x})\frac{{\mathit{U}}_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x})}{2}+f_{j-\frac{1}{2}}(\mathit{x}).\end{array}$(29) is a second-order boundary differential equation that can be described iteratively by

(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30)

with the boundary conditions,

(31) $\{\begin{array}{l}{U}_j(0)=0,\\ U_j(1)=0\end{array}\begin{array}{l}\mathit{Giventhat}{U}_{\mathit{j}-\mathit{m}}(\mathit{x})={\gamma }_b(\mathit{x},t_{\mathit{j}-\mathit{m}}),\\ \mathit{when}0\le \mathit{j}\le \mathit{m}.\end{array}$(31)

(30) and (31) constitute the semi-discrete problem of problem (28)-(2). Once again, we define a differential operator,

(32) ${\mathcal{L}}_{\epsilon ,\mu }^M\equiv -\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x})).$(32)

Considering the semi-discretization, at each time step j, the operator ${\mathcal{L}}_{\epsilon ,\mu }^M$ satisfies the maximum principle.

Lemma 3.1.

Semi-discrete maximum principle:

Suppose $\mathit{\chi }(\mathit{x},t_{\mathit{j}})\in {\mathcal{C}}^2(0,1)$, such that $\mathit{\chi }(0,t_{\mathit{j}}),\mathit{\chi }(1,t_{\mathit{j}})\ge 0,$ and ${\mathcal{L}}_{\epsilon ,\mu }^M(\mathit{\chi }(\mathit{x},t_{\mathit{j}}))\ge 0,\mathrm{\forall x}\in (0,1)$, then $\mathit{\chi }(\mathit{x},t_{\mathit{j}})\ge 0,\mathrm{\forall x}\in [0,1]$.

Proof.

Ref (Worku et al., 2023)

Lemma 3.2.

Semi-disretization Uniform Stability

The solution $\mathit{u}(\mathit{x},t_{\mathit{j}})$ of (30) - (31) satisfies ${\displaystyle \Vert \mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \frac{1}{\beta }\Vert \mathit{g}\Vert ,}$ where

$\begin{array}{ll}\mathit{g}(\mathit{x},t_{\mathit{j}})& =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})\\ +(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}).\end{array}$

Proof.

Define barrier functions

$\begin{array}{l}{\theta }^{\pm }(\mathit{x},t_{\mathit{j}})=\frac{1}{\beta }\Vert \mathit{g}\Vert \pm \mathit{u}(\mathit{x},t_{\mathit{j}})\end{array}$

$\begin{array}{l}{\theta }^{\pm }(0,t_{\mathit{j}})=\frac{1}{\beta }\Vert \mathit{g}\Vert \pm 0\ge 0,{\theta }^{\pm }(1,t_{\mathit{j}})=\frac{1}{\beta }\Vert \mathit{g}\Vert \pm 0\ge 0\text{ and}\end{array}$

$\begin{array}{ll}{\mathcal{L}}_{\epsilon ,\mu }^M({\theta }^{\pm }(\mathit{x},t_{\mathit{j}}))& =\pm {\mathcal{L}}_{\epsilon ,\mu }^M\mathit{u}(\mathit{x},t_{\mathit{j}})+(\frac{2}{\mathrm{\Delta t}}+\mathit{b}).(\frac{1}{\beta }\Vert \mathit{g}\Vert )\\ =\pm \mathit{g}(\mathit{x},t_{\mathit{j}})+\frac{b}{\beta }\Vert \mathit{g}\Vert +\frac{2}{\Delta \mathit{t\beta }}\Vert \mathit{g}\Vert \\ \ge 0.\end{array}$

${\displaystyle \frac{1}{\beta }\Vert \mathit{g}\Vert \pm \mathit{u}(\mathit{x},t_{\mathit{j}})\ge 0}$ by the semi-discrete maximum principle of Lemma  3.1.

Then

$\begin{array}{l}\Vert \mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \frac{1}{\beta }\Vert \mathit{g}\Vert .\end{array}$

Let $U_{\mathit{j}}(\mathit{x})$ be the numerical solution of the semi-discrete problem (30(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30) ), whose exact solution is $\mathit{u}(\mathit{x},t_{\mathit{j}})$ at a fixed $\mathit{t}=t_{\mathit{j}}.$ $U_{\mathit{j}}(\mathit{x})$ is generated in an iterative sequence from the prior numerical solution, $U_{\mathit{j}-1}(\mathit{x}).$ Assume $U_j^{\ast }(\mathit{x})$ is generated from the exact prior solution $\mathit{u}(\mathit{x},t_{\mathit{j}-1}),$ then the following holds.

(33) $\begin{array}{l}\frac{\mathrm{\Delta t}}{2}\{-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}+b_{j-\frac{1}{2}}(\mathit{x})\}{U}_j^{\ast }(\mathit{x})+U_j^{\ast }(\mathit{x})\\ =\frac{\mathrm{\Delta t}}{2}\{\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}-b_{j-\frac{1}{2}}(\mathit{x})\}\mathit{u}(\mathit{x},t_{\mathit{j}-1})+\mathit{u}(\mathit{x},t_{\mathit{j}-1})\\ +\mathrm{\Delta t}{c}_{j-\frac{1}{2}}(\mathit{x})\left(\frac{{\mathit{U}}_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x})}{2}\right)+\mathrm{\Delta t}{f}_{j-\frac{1}{2}}(\mathit{x}).\end{array}$(33)

Then, $e_j=|U_{\mathit{j}}(\mathit{x})-U_j^{\ast }(\mathit{x})|$ is the absolute local error, the error caused by one step of iteration. The bound of this error is measured in the $(\mathit{j}{)}^{\mathit{th}}$ local truncation error Te_j given by

$\begin{array}{ll}\mathit{T}{e}_{\mathit{j}}& ={\mathcal{L}}_{\epsilon ,\mu }^M(U_{\mathit{j}}(\mathit{x})-U_j^{\ast }(\mathit{x}))\\ =\frac{\mathrm{\Delta t}}{2}\{-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}+b_{j-\frac{1}{2}}(\mathit{x})\}\mathit{u}(\mathit{x},t_{\mathit{j}})+\mathit{u}(\mathit{x},t_{\mathit{j}})\\ -\frac{\mathrm{\Delta t}}{2}\{\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}-b_{j-\frac{1}{2}}(\mathit{x})\}\mathit{u}(\mathit{x},t_{\mathit{j}-1})-\mathit{u}(\mathit{x},t_{\mathit{j}-1})\\ -\Delta \mathit{tc}(\mathit{x},t_{j-\frac{1}{2}})\left(\frac{\mathit{u}(\mathit{x},{\mathit{t}}_{\mathit{j}-1-\mathit{m}})+\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})}{2}\right)-\Delta \mathit{tf}(\mathit{x},t_{j-\frac{1}{2}})\end{array}$

$\begin{array}{ll}|\mathit{T}{e}_{\mathit{j}}|& \le \mathrm{\Delta t}|(-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu a}(\mathit{x},t_{j-\frac{1}{2}})\frac{d}{\mathit{dx}}+\mathit{b}(\mathit{x},t_{j-\frac{1}{2}}))\\ \times \left(\frac{\mathit{u}(\mathit{x},{\mathit{t}}_{\mathit{j}})+\mathit{u}(\mathit{x},t_{\mathit{j}-1})}{2}\right)+\left(\frac{\mathit{u}(\mathit{x},{\mathit{t}}_{\mathit{j}})-\mathit{u}(\mathit{x},t_{\mathit{j}-1})}{\mathrm{\Delta t}}\right)\\ -(\mathit{c}(\mathit{x},t_{j-\frac{1}{2}})\mathit{u}(\mathit{x},t_{j-\frac{1}{2}-\mathit{m}})+\mathit{f}(\mathit{x},t_{j-\frac{1}{2}}))|\\ \le \mathrm{\Delta t}|(-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}-\mathit{\mu a}(\mathit{x},t_{j-\frac{1}{2}})\frac{d}{\mathit{dx}}+\mathit{b}(\mathit{x},t_{j-\frac{1}{2}}))(\mathit{u}(\mathit{x},t_{j-\frac{1}{2}}))\\ +\left(\frac{\mathit{u}(\mathit{x},{\mathit{t}}_{\mathit{j}})-\mathit{u}(\mathit{x},t_{\mathit{j}-1})}{\mathrm{\Delta t}}\right)-{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}\mathit{u}(\mathit{x},t_{j-\frac{1}{2}})|\\ \le \mathrm{\Delta t}|\frac{\mathit{u}(\mathit{x},{\mathit{t}}_{\mathit{j}})-\mathit{u}(\mathit{x},t_{\mathit{j}-1})}{\mathrm{\Delta t}}-\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},t_{j-\frac{1}{2}})|\\ \le \mathrm{\Delta t}|\frac{(\mathrm{\Delta }t{)}^2}{24}\frac{{\partial }^3}{\mathit{\partial }{t}^3}\mathit{u}(\mathit{x},t^{\ast })|,\begin{array}{l}\text{for some }{t}^{\ast }\text{such that }\\ t_{\mathit{j}-1}\le t^{\ast }\le t_j\end{array}\\ \le \mathit{C}(\mathrm{\Delta t}{)}^3\end{array}$

This argument with Lemma 3.1 implies the error estimat of local error,

(34) $\Vert e_j\Vert \le \mathit{C}(\mathrm{\Delta t}{)}^3.$(34)

Lemma 3.3.

The global error estimate, E of the semi-discrete problem (30(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30) ) and (31(31) $\{\begin{array}{l}{U}_j(0)=0,\\ U_j(1)=0\end{array}\begin{array}{l}\mathit{Giventhat}{U}_{\mathit{j}-\mathit{m}}(\mathit{x})={\gamma }_b(\mathit{x},t_{\mathit{j}-\mathit{m}}),\\ \mathit{when}0\le \mathit{j}\le \mathit{m}.\end{array}$(31) ) is given by

(35) $\Vert \mathit{E}\Vert \le \mathit{C}(\mathrm{\Delta t}{)}^2$(35)

Proof.

From (34(34) $\Vert e_j\Vert \le \mathit{C}(\mathrm{\Delta t}{)}^3.$(34) ) we have $\Vert e_j\Vert \le \mathit{C}(\mathrm{\Delta t}{)}^3$.

$\Vert e_j\Vert$ is the contribution of each time step error to the global error, then the global error, E will be

$\begin{array}{ll}{\displaystyle \Vert \mathit{\text{E}}\Vert }& =\Vert \sum _{\mathit{j}=1}^{\mathit{M}}{e}_j\Vert \le \sum _{\mathit{j}=1}^{\mathit{M}}\Vert e_j\Vert \\ \le \sum _{\mathit{j}=1}^{\mathit{M}}\mathit{C}(\mathrm{\Delta t}{)}^3=\mathit{CM}(\mathrm{\Delta t}{)}^3=\mathit{CT}(\mathrm{\Delta t}{)}^2=\mathit{C}(\mathrm{\Delta t}{)}^2.\end{array}$

Then, $\Vert \mathit{\text{E}}\Vert \le \mathit{C}(\mathrm{\Delta t}{)}^2$

The boundary layer properties are determined according to the characteristic equation of 30(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30) which is

(36) $-\mathit{\epsilon }(\mathit{\lambda }(\mathit{x}){)}^2-\mathit{\mu a\lambda }(\mathit{x})+(\mathit{b}+\frac{2}{\mathrm{\Delta t}})=0,$(36)

whose roots are

$\begin{array}{ll}{\lambda }_1(\mathit{x})& =-\frac{\mathit{\mu a}}{2\mathit{\epsilon }}[1+\sqrt{1+\frac{4\mathit{\epsilon }}{{\mu }^2{a}^2}(\mathit{b}+\frac{2}{\mathrm{\Delta t}})}]\\ {\lambda }_2(\mathit{x})& =-\frac{\mathit{\mu a}}{2\mathit{\epsilon }}[1-\sqrt{1+\frac{4\mathit{\epsilon }}{{\mu }^2{a}^2}(\mathit{b}+\frac{2}{\mathrm{\Delta t}})}].\end{array}$

let

(37) ${\mathrm{\Lambda }}_1=-\underset{x}{max}\{{\lambda }_1(\mathit{x})\}\mathit{and}{\mathrm{\Lambda }}_2=\underset{x}{min}\{{\lambda }_2(\mathit{x})\}$(37)

Lemma 3.4.

The i^th derivative of the solution, $\mathit{u}(\mathit{x},t_{\mathit{j}})$ of (30) - (31) satisfies ${\displaystyle |\frac{d^i}{\mathit{d}{x}^i}\mathit{u}(\mathit{x},t_{\mathit{j}})|\le \mathit{C}[1+{\mathrm{\Lambda }}_1^i\mathit{exp}(-\mathit{\rho }{\mathrm{\Lambda }}_1\mathit{x})+{\mathrm{\Lambda }}_2^i\mathit{exp}(-\mathit{\rho }{\mathrm{\Lambda }}_2(1-\mathit{x}))]}$, for $0\le \mathit{i}\le \mathit{q}$ up to a certain prescribed order q, and any ρ,$0<\mathit{\rho }<1$.

Proof.

Refer (M. Kadalbajoo & Yadaw, 2012; Rajan & Reddy, 2022)

3.2. Spatial discretization

The boundary layers’ widths depend on the values of perturbation parameters, while their position is at the left and right neighbourhoods of x = 0 and x = 1. Since we are implementing the fitted mesh method, the spatial discretization is piecewise uniform Shishkin mesh type. The spatial domain component, $[0,1]$ is split into three subintervals: $[0,{\sigma }_1],[{\sigma }_1,1-{\sigma }_2]$ and $[1-{\sigma }_2,1]$, assuming the left transition point is σ₁ and the right is σ₂. The boundary layer subintervals $[0,{\sigma }_1]$ and $[1-{\sigma }_2,1]$ will have $\frac{N}{4}$ fine elements each and $\frac{N}{2}$coarse elements will distribute along the outer layer region, $[{\sigma }_1,1-{\sigma }_2]$ uniformly, where a 4^th multiple N is number of elements in spatial direction (Miller et al., 1998). The transition points σ₁ and σ₂ are determined by

(38) $\begin{array}{ll}{\sigma }_1& =\{\begin{array}{ll}min\{\frac{1}{4},\frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\varsigma \alpha }}\mathit{lnN}\},& {\mu }^2\le \frac{\mathit{\varsigma \epsilon }}{\alpha },\\ min\{\frac{1}{4},\frac{4\mathit{\epsilon }}{\mathit{\mu \alpha }}\mathit{lnN}\},& {\mu }^2\ge \frac{\mathit{\varsigma \epsilon }}{\alpha }\end{array}\text{and }\\ {\sigma }_2& =\{\begin{array}{ll}min\{\frac{1}{4},\frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\varsigma \alpha }}\mathit{lnN}\},& {\mu }^2\le \frac{\mathit{\varsigma \epsilon }}{\alpha },\\ min\{\frac{1}{4},\frac{4\mathit{\mu }}{\varsigma }\mathit{lnN}\},& {\mu }^2\ge \frac{\mathit{\varsigma \epsilon }}{\alpha }.\end{array}\end{array}$(38)

The sizes of piecewise uniform meshes are provided by

(39) $h_i=\{\begin{array}{ll}\frac{4{\sigma }_1}{\mathit{N}},& \mathit{i}=1,2,\dots ,\frac{N}{4}+1,\\ \frac{2(1-{\sigma }_1-{\sigma }_2)}{\mathit{N}},& \mathit{i}=\frac{N}{4}+1,\dots ,\frac{3\mathit{N}}{4}+1,\\ \frac{4{\sigma }_2}{\mathit{N}},& \mathit{i}=\frac{3\mathit{N}}{4}+1,\dots ,\mathit{N}+1.\end{array}$(39)

Accordingly the node points are given as

(40) $x_i=\{\begin{array}{ll}\frac{4}{N}{\sigma }_1(\mathit{i}-1),& \mathit{i}=1,2,\dots ,\frac{N}{4}+1,\\ {\sigma }_1+\frac{2}{N}(1-{\sigma }_1-{\sigma }_2)(\mathit{i}-1-\frac{N}{4}),& \mathit{i}=\frac{N}{4}+1,\dots ,\frac{3\mathit{N}}{4}+1,\\ 1-{\sigma }_2+\frac{4}{N}{\sigma }_2(\mathit{i}-1-\frac{3\mathit{N}}{4}),& \mathit{i}=\frac{3\mathit{N}}{4}+1,\dots ,\mathit{N}+1.\end{array}$(40)

So, the spatial discretization , ${\mathcal{D}}_x^N=\{x_i=x_{\mathit{i}-1}+h_i,x_1=0\text{and }\mathit{i}=2,\dots ,\mathit{N}+1\}.$ The following inequality needs to be considered for error analysis:

(41) $\begin{array}{l}1<\mathit{i}<\mathit{N},h_i\le \frac{2}{N},\end{array}$(41)

(42) $\begin{array}{l}1<\mathit{i}<\frac{N}{4},h_i\le \frac{4{\sigma }_1}{\mathit{N}}\le \frac{1}{N},\end{array}$(42)

(43) $\begin{array}{l}\frac{3\mathit{N}}{4}<\mathit{i}<\mathit{N},h_i\le \frac{4{\sigma }_2}{\mathit{N}}\le \frac{1}{N},\end{array}$(43)

(44) $\begin{array}{l}\frac{N}{4}<\mathit{i}<\frac{3\mathit{N}}{4},h_i\ge \frac{1}{N}.\end{array}$(44)

Across the spatial direction, we consider the quadratic base functions. For each individual element there are three local nodes, which implies $2\mathit{N}+1$ global nodes through the spatial domain. Let ${\displaystyle \{{\varphi }_i(\mathit{x}){\}}_1^{2N}}$ be the set of base function and ${\displaystyle V_2({\overline{\mathcal{D}}}_{\mathit{x}}^{2\mathit{N}})}$ be the space of second degree polynomials those condense at the boundaries x = 0 and x = 1. ${\displaystyle V_2({\overline{\mathcal{D}}}_{\mathit{x}}^{\mathit{\text{2N}}})=\mathit{Span}\{{\varphi }_j(\mathit{x}){\}}_1^{2N}}$ and set a Sobolev space $\mathit{V}=H_0^1({\overline{\mathcal{D}}}_{\mathit{x}}^{\mathit{\text{2N}}}).$ Clearly, at each t, the piecewise quadratic approximation $U_j(\mathit{x})\in V_2({\overline{\mathcal{D}}}_{\mathit{x}}^{\mathit{\text{2M}}})\subset \mathit{V}$. The finite element approximation of $\mathit{u}(\mathit{x},t_j)$ is given as linear combination of these base functions.

(45) $\mathit{u}(\mathit{x},t_j)\approx U_j(\mathit{x})=\sum _{\mathit{i}=1}^{2\mathit{N}}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}),$(45)

where $U_{\mathit{i},\mathit{j}}$ is the functional value of $\mathit{u}(x_i,t_j)$ which is going to be determined later. The required two variable finite element solution is obtained by substituting (45(45) $\mathit{u}(\mathit{x},t_j)\approx U_j(\mathit{x})=\sum _{\mathit{i}=1}^{2\mathit{N}}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}),$(45) ) into the Rothe’s function (22(22) $\mathit{u}(\mathit{x},\mathit{t})\approx \mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}}{U}_j(\mathit{x}){\psi }_j(\mathit{t}),$(22) ).

(46) $\mathit{u}(\mathit{x},\mathit{t})\approx \mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}}\sum _{\mathit{i}=1}^{2\mathit{N}}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}){\psi }_j(\mathit{t}).$(46)

Now, to determine the unknown coefficients $U_{\mathit{i},\mathit{j}}$, write the weak or variation form of (30(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30) ) by multiplying a piecewise smooth weight function w(x) which vanishes at the boundary ${\mathrm{\Gamma }}_l\cup {\mathrm{\Gamma }}_r$ and integrate over $[0,1]$. Since the Dirichlet boundary conditions are homogeneous, integration by part implies

(47) $\begin{array}{l}\mathit{\epsilon }⟨\frac{d{U}_j}{\mathit{dx}}.\frac{\mathit{dw}}{\mathit{dx}}⟩-\mathit{\mu }\overline{a}⟨\frac{d{U}_j}{\mathit{dx}}.\mathit{w}⟩+(\frac{2}{\mathrm{\Delta t}}+\overline{b})⟨U_j,\mathit{w}⟩=\\ -\mathit{\epsilon }⟨\frac{\mathit{d}{\mathit{U}}_{\mathit{j}-1}}{\mathit{dx}}.\frac{\mathit{dw}}{\mathit{dx}}⟩+\mathit{\mu }\overline{a}⟨\frac{\mathit{d}{\mathit{U}}_{\mathit{j}-1}}{\mathit{dx}}.\mathit{w}⟩+(\frac{2}{\mathrm{\Delta t}}-\overline{b})⟨U_{\mathit{j}-1},\mathit{w}⟩\\ +\overline{c}⟨U_{\mathit{j}-\mathit{m}}+U_{\mathit{j}-1-\mathit{m}},\mathit{w}⟩+2⟨f_{j-\frac{1}{2}}(\mathit{x}),\mathit{w}⟩.\end{array}$(47)

In the case of Galerkin finite element method, the weight function w(x) is chosen to be equal to the approximation function $\mathit{\varphi }(\mathit{x})$. So, the algebraic equation of 47(47) $\begin{array}{l}\mathit{\epsilon }⟨\frac{d{U}_j}{\mathit{dx}}.\frac{\mathit{dw}}{\mathit{dx}}⟩-\mathit{\mu }\overline{a}⟨\frac{d{U}_j}{\mathit{dx}}.\mathit{w}⟩+(\frac{2}{\mathrm{\Delta t}}+\overline{b})⟨U_j,\mathit{w}⟩=\\ -\mathit{\epsilon }⟨\frac{\mathit{d}{\mathit{U}}_{\mathit{j}-1}}{\mathit{dx}}.\frac{\mathit{dw}}{\mathit{dx}}⟩+\mathit{\mu }\overline{a}⟨\frac{\mathit{d}{\mathit{U}}_{\mathit{j}-1}}{\mathit{dx}}.\mathit{w}⟩+(\frac{2}{\mathrm{\Delta t}}-\overline{b})⟨U_{\mathit{j}-1},\mathit{w}⟩\\ +\overline{c}⟨U_{\mathit{j}-\mathit{m}}+U_{\mathit{j}-1-\mathit{m}},\mathit{w}⟩+2⟨f_{j-\frac{1}{2}}(\mathit{x}),\mathit{w}⟩.\end{array}$(47) will be

(48) $\begin{array}{l}\sum _{\mathit{i}=1}^{2\mathit{N}}\{\mathit{\epsilon }⟨\frac{d{\varphi }_i}{\mathit{dx}},\frac{d{\varphi }_l}{\mathit{dx}}⟩-\mathit{\mu }\overline{a}⟨\frac{d{\varphi }_i}{\mathit{dx}},{\varphi }_l⟩\\ +(\frac{2}{\mathrm{\Delta t}}+\overline{b})⟨{\varphi }_i,{\varphi }_l⟩\}{U}_{\mathit{i},\mathit{j}}\\ =\sum _{\mathit{i}=1}^{2\mathit{N}}\{-\mathit{\epsilon }⟨\frac{d{\varphi }_i}{\mathit{dx}},\frac{d{\varphi }_l}{\mathit{dx}}⟩+\mathit{\mu }\overline{a}⟨\frac{d{\varphi }_i}{\mathit{dx}},{\varphi }_l⟩\\ +(\frac{2}{\mathrm{\Delta t}}-\overline{b})⟨{\varphi }_i,{\varphi }_l⟩\}{U}_{\mathit{i},\mathit{j}-1}\\ +\overline{c}\sum _{\mathit{i}=1}^{2\mathit{N}}(U_{\mathit{i},\mathit{j}-1-\mathit{m}}+U_{\mathit{i},\mathit{j}-\mathit{m}})⟨{\varphi }_i,{\varphi }_l⟩+2⟨f_j(\mathit{x}),{\varphi }_l⟩,\end{array}$(48)

for each $\mathit{j}=1,2,\dots ,\mathit{M}+1,$ for all $\mathit{i},\mathit{\text{l}}=1,2,\dots ,2\mathit{N}$ and the bar functions $\overline{a}=a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\approx \mathit{a}(x_i+\frac{1}{2}{h}_i,t_j+\frac{1}{2}\mathrm{\Delta t})$ and $\overline{b}=b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\approx \mathit{b}(x_i+\frac{1}{2}{h}_i,t_j+\frac{1}{2}\mathrm{\Delta t})$.

48 can be expressed in a matrix notation of the form

(49) $\mathit{A}{U}^{\mathit{j}}=\mathit{B}{U}^{\mathit{j}-1}+\mathit{C}(U^{\mathit{j}-1-\mathit{m}}+U^{\mathit{j}-\mathit{m}})+\mathbf{\text{textit{f}}},$(49)

where the matrices, $\mathit{A}=(a_{\mathit{i},\mathit{\text{l}}})$, $\mathit{B}=(b_{\mathit{i},\mathit{\text{l}}})$, $\mathit{C}=(c_{\mathit{i},\mathit{\text{l}}})$ and the vector $\mathbf{\text{textit{f}}}=(f_{\mathit{\text{l}}})$, whose entries are given bellow.

$\begin{array}{ll}{\displaystyle a_{\mathit{i},\mathit{\text{l}}}}& =\mathit{\epsilon }⟨\frac{d{\varphi }_i}{\mathit{dx}},\frac{d{\varphi }_l}{\mathit{dx}}⟩-\mathit{\mu }\overline{a}⟨\frac{d{\varphi }_i}{\mathit{dx}},{\varphi }_l⟩+(\frac{2}{\mathrm{\Delta t}}+\overline{b})⟨{\varphi }_i,{\varphi }_l⟩,\\ b_{\mathit{i},\mathit{\text{l}}}& =-\mathit{\epsilon }⟨\frac{d{\varphi }_i}{\mathit{dx}},\frac{d{\varphi }_l}{\mathit{dx}}⟩+\mathit{\mu }\overline{a}⟨\frac{d{\varphi }_i}{\mathit{dx}},{\varphi }_l⟩+(\frac{2}{\mathrm{\Delta t}}-\overline{b})⟨{\varphi }_i,{\varphi }_l⟩,\\ c_{\mathit{i},\mathit{\text{l}}}& =\overline{c}⟨{\varphi }_i,{\varphi }_l⟩\text{and }\\ f_{\mathit{\text{l}}}& ={\displaystyle 2⟨f_j(\mathit{x}),{\varphi }_l⟩.}\end{array}$

3.2.1. Formulation of element equation

Consider an arbitrary space element $\mathit{e}=[x_{\mathit{i}-1},x_{\mathit{i}}]$ and a quadratic approximation for this particular element, ${\displaystyle U_j^e(\mathit{x})=U_{\mathit{i}-1,\mathit{j}}{\varphi }_1^e(\mathit{x})+U_{i-\frac{1}{2},\mathit{j}}{\varphi }_2^e(\mathit{x})+U_{\mathit{i},\mathit{j}}{\varphi }_3^e(\mathit{x})}$. We need three local nodes to determine $U_j^e(\mathit{x})$ in terms of the Lagrange quadratic interpolation functions and the nodal values. Take the two endpoints ($x_{\mathit{i}-1}=x_1^e,x_{\mathit{i}}=x_3^e$) and the midpoint ($x_{i-\frac{1}{2}}=x_2^e$) for convenience.

The parity and unity properties of Lagrange interpolation functions ${\varphi }_1^e(\mathit{x}),{\varphi }_2^e(\mathit{x}),{\varphi }_3^e(\mathit{x})$,

(50) ${\varphi }_p^e(x_q^e)=\{\begin{array}{l}0,\text{if }\mathit{p}\ne \mathit{q},\\ 1,\text{if }\mathit{p}=\mathit{q},\end{array}\mathit{p},\mathit{q}=1,2,3\text{and }\sum _{\mathit{p}=1}^3{\varphi }_p^e(\mathit{x})=1$(50)

help us to derive the Lagrange quadratic interpolation polynomials as

(51) $\begin{array}{ll}{\varphi }_1^e(\mathit{x})& =\frac{1}{D}[x_{i-\frac{1}{2}}(x_i{)}^2-x_i(x_{i-\frac{1}{2}}{)}^2+((x_{i-\frac{1}{2}}{)}^2-(x_i{)}^2)\mathit{x}\\ +((x_i)-x_{i-\frac{1}{2}})x^2],\end{array}$(51)

(52) $\begin{array}{ll}{\varphi }_2^e(\mathit{x})& =\frac{1}{D}[x_{\mathit{i}}(x_{\mathit{i}-1}{)}^2-x_{\mathit{i}-1}(x_{\mathit{i}}{)}^2+((x_{\mathit{i}}{)}^2-(x_{\mathit{i}-1}{)}^2)\mathit{x}\\ +(x_{\mathit{i}-1}-x_{\mathit{i}})x^2],\end{array}$(52)

(53) $\begin{array}{ll}{\varphi }_3^e(\mathit{x})& =\frac{1}{D}[x_{\mathit{i}-1}(x_{i-\frac{1}{2}}{)}^2-x_{i-\frac{1}{2}}(x_{\mathit{i}-1}{)}^2+((x_{\mathit{i}-1}{)}^2-(x_{i-\frac{1}{2}}{)}^2)\mathit{x}\\ +(x_{i-\frac{1}{2}}-x_{\mathit{i}-1})x^2],\end{array}$(53)

where $\mathit{D}=x_{i-\frac{1}{2}}(x_i{)}^2-x_i(x_{i-\frac{1}{2}}{)}^2+x_i(x_{\mathit{i}-1}{)}^2-x_{\mathit{i}-1}(x_i{)}^2+x_{\mathit{i}-1}(x_{i-\frac{1}{2}}{)}^2-x_{i-\frac{1}{2}}(x_{\mathit{i}-1}{)}^2.$

The inner product terms of the weak form with in this particular element are computed as follows.

(54) $\begin{array}{ll}{\displaystyle ⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},\frac{d{\varphi }_q}{\mathit{dx}}⟩}& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}}.\frac{d{\varphi }_q}{\mathit{dx}}\mathit{dx}=\frac{1}{3h_i}\left[\begin{array}{lll}7& -8& 1\\ -8& 16& -8\\ 1& -8& 7\end{array}\right],\end{array}$(54)

(55) $\begin{array}{ll}{\displaystyle ⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},{\varphi }_q⟩}& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}}.{\varphi }_q\mathit{dx}=\frac{1}{6}\left[\begin{array}{lll}-3& -4& 1\\ 4& 0& -4\\ -1& 4& 3\end{array}\right],\end{array}$(55)

(56) $\begin{array}{ll}{\displaystyle ⟨{\varphi }_{\mathit{p}},{\varphi }_q⟩}& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}{\varphi }_{\mathit{p}}.{\varphi }_q\mathit{dx}=\frac{h_i}{30}\left[\begin{array}{lll}4& 2& -1\\ 2& 16& 2\\ -1& 2& 4\end{array}\right],\end{array}$(56)

(57) $\begin{array}{ll}⟨f_j(\mathit{x}),{\varphi }_q⟩& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}{f}_j(\mathit{x}).{\varphi }_q\mathit{dx}=\left[\begin{array}{l}{\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_1\mathit{dx}\\ {\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_2\mathit{dx}\\ {\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_3\mathit{dx}\end{array}\right].\end{array}$(57)

The element matrices A^e and B^e are given by

$\begin{array}{ll}{A}^e& =\mathit{\epsilon }⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},\frac{d{\varphi }_q}{\mathit{dx}}⟩-\mathit{\mu }{a}_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},{\varphi }_q⟩\\ +(\frac{2}{\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}})⟨{\varphi }_{\mathit{p}},{\varphi }_q⟩,\end{array}$

$\begin{array}{ll}{B}^e& =-\mathit{\epsilon }⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},\frac{d{\varphi }_q}{\mathit{dx}}⟩+\mathit{\mu }{a}_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},{\varphi }_q⟩\\ +(\frac{2}{\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}})⟨{\varphi }_{\mathit{p}},{\varphi }_q⟩.\end{array}$

The results through (54(54) $\begin{array}{ll}{\displaystyle ⟨\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}},\frac{d{\varphi }_q}{\mathit{dx}}⟩}& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\frac{\mathit{d}{\mathit{\varphi }}_{\mathit{p}}}{\mathit{dx}}.\frac{d{\varphi }_q}{\mathit{dx}}\mathit{dx}=\frac{1}{3h_i}\left[\begin{array}{lll}7& -8& 1\\ -8& 16& -8\\ 1& -8& 7\end{array}\right],\end{array}$(54) ) - (57(57) $\begin{array}{ll}⟨f_j(\mathit{x}),{\varphi }_q⟩& ={\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}{f}_j(\mathit{x}).{\varphi }_q\mathit{dx}=\left[\begin{array}{l}{\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_1\mathit{dx}\\ {\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_2\mathit{dx}\\ {\int }_{{\mathit{x}}_{\mathit{i}-1}}^{x_i}\mathit{f}(\mathit{x},t_j).{\varphi }_3\mathit{dx}\end{array}\right].\end{array}$(57) ) imply

(58) $A^e=\left[\begin{array}{lll}{a}_{1,1}^e& a_{1,2}^e& a_{1,3}^e\\ a_{2,1}^e& a_{2,2}^e& a_{2,3}^e\\ a_{3,1}^e& a_{3,2}^e& a_{3,3}^e\end{array}\right]\text{and }{B}^e=\left[\begin{array}{lll}{b}_{1,1}^e& b_{1,2}^e& b_{1,3}^e\\ b_{2,1}^e& b_{2,2}^e& b_{2,3}^e\\ b_{3,1}^e& b_{3,2}^e& b_{3,3}^e\end{array}\right]$(58)

where the entries are

$\begin{array}{ll}{a}_{1,1}^e& =\frac{7\mathit{\epsilon }}{3h_i}+\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{2}+\frac{4h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{2h_i}{15},\\ b_{1,1}^e& =-\frac{7\mathit{\epsilon }}{3h_i}-\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{2}+\frac{4h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{2h_i}{15}\\ a_{1,2}^e& =-\frac{8\mathit{\epsilon }}{3h_i}+\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15},\\ b_{1,2}^e& =\frac{8\mathit{\epsilon }}{3h_i}-\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15}\\ a_{1,3}^e& =\frac{\epsilon }{3h_i}-\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{6}-\frac{h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{30},\\ b_{1,3}^e& =-\frac{\epsilon }{3h_i}+\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{6}-\frac{h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{30}\\ a_{2,1}^e& =-\frac{8\mathit{\epsilon }}{3h_i}-\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15},\\ b_{2,1}^e& =\frac{8\mathit{\epsilon }}{3h_i}+\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15}\\ a_{2,2}^e& =\frac{16\mathit{\epsilon }}{3h_i}+0+\frac{16h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{8h_i}{15},\\ b_{2,2}^e& =-\frac{16\mathit{\epsilon }}{3h_i}+0+\frac{16h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{8h_i}{15}\\ a_{2,3}^e& =-\frac{8\mathit{\epsilon }}{3h_i}+\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15},\\ b_{2,3}^e& =\frac{8\mathit{\epsilon }}{3h_i}-\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15}\\ a_{3,1}^e& =\frac{\epsilon }{3h_i}+\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{6}-\frac{h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{30},\\ b_{3,1}^e& =-\frac{\epsilon }{3h_i}+\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{6}-\frac{h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{30},\\ a_{3,2}^e& =-\frac{8\mathit{\epsilon }}{3h_i}-\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15},\\ b_{3,2}^e& =\frac{8\mathit{\epsilon }}{3h_i}+\frac{2\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{3}+\frac{2h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{h_i}{15}\\ a_{3,3}^e& =\frac{7\mathit{\epsilon }}{3h_i}-\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{2}+\frac{4h_i}{15\mathrm{\Delta t}}+b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{2h_i}{15},\\ b_{3,3}^e& =-\frac{7\mathit{\epsilon }}{3h_i}+\frac{\mu a_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{2}+\frac{4h_i}{15\mathrm{\Delta t}}-b_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}\frac{2h_i}{15}.\end{array}$

The coefficient matrix of delay term,

(59) $C^e=\frac{h_i{c}_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{30}\left[\begin{array}{lll}4& 2& -1\\ 2& 16& 2\\ -1& 2& 4\end{array}\right].$(59)

The element source vector,

(60) ${\mathit{f}}^{\mathit{e}}=\frac{h_i{f}_{i-\frac{1}{2},\mathit{j}-\frac{1}{2}}}{6}\left[\begin{array}{l}1\\ 4\\ 1\end{array}\right].$(60)

Then, an element equation of the arbitrary space element $\mathit{e}=[x_{\mathit{i}-1},x_{\mathit{i}}]$ recurrently given by

(61) $A^e{U}_j^e=B^e{U}_{j-1}^e+C^e(U_{j-1-m}^e+U_{j-m}^e)+{\mathit{f}}^e.$(61)

3.2.2. Assembling element equations

If the linear elements are connected head to tail, then the assembled equation will be a $2\mathit{N}-1\times 2\mathit{N}-1$ system of equation of the form

(62) $\mathit{A}{U}^j=\mathit{B}.U^{\mathit{j}-1}+\mathit{C}(U^{\mathit{j}-1-\mathit{m}}+U^{\mathit{j}-\mathit{m}})+\mathit{f},$(62)

for each j^th time step with the discrete initial and boundary conditions,

(63) $\begin{array}{ll}{U}_{\mathit{i},\mathit{j}}& ={\gamma }_b(x_i,t_j)\mathrm{\forall ij}=-\mathit{m},-\mathit{m}+1,\cdots ,0,\\ U_{0,\mathit{j}}& =0\mathrm{\forall j},\\ U_{1,\mathit{j}}& =0\mathrm{\forall j}.\end{array}$(63)

The stiffness matrix A is a penta-diagonal given by

(64) $\mathit{A}=\left[\begin{array}{lllllll}{a}_{1,1}^1& a_{1,2}^1& a_{1,3}^1& 0& 0& 0& 0\\ a_{2,1}^1& a_{2,2}^1& a_{2,3}^1& 0& 0& 0& 0\\ a_{3,1}^1& a_{3,2}^1& a_{3,3}^1+a_{1,1}^2& a_{1,2}^2& a_{1,3}^2& 0& 0\\ 0& 0& a_{2,1}^2& a_{2,2}^2& a_{2,3}^2& 0& 0\\ 0& 0& a_{3,1}^2& a_{3,2}^2& a_{3,3}^2+a_{1,1}^3& a_{1,2}^3& a_{1,3}^2\\ 0& 0& 0& 0& a_{2,1}^3& a_{2,2}^3& a_{2,3}^3\\ 0& 0& 0& 0& a_{3,1}^2& a_{3,2}^2& a_{3,3}^2+a_{1,1}^3\\ \ddots \\ a_{3,3}^{N+1}\end{array}\right].$(64)

The matrices B and C will be established in the same manner and the source vector function,

(65) $\mathbf{\text{textit{f}}}=\frac{1}{6}\left[\begin{array}{l}{h}_1{f}_{\frac{1}{2},\mathit{j}-\frac{1}{2}}\\ 4h_1{f}_{\frac{1}{2},\mathit{j}-\frac{1}{2}}\\ h_1{f}_{\frac{1}{2},\mathit{j}-\frac{1}{2}}+h_2{f}_{\frac{3}{2},\mathit{j}-\frac{1}{2}}\\ 4h_2{f}_{\frac{3}{2},\mathit{j}-\frac{1}{2}}\\ h_2{f}_{\frac{3}{2},\mathit{j}-\frac{1}{2}}+h_3{f}_{\frac{5}{2},\mathit{j}-\frac{1}{2}}\\ ⋮\\ h_N{f}_{\frac{2\mathit{N}-1}{2},\mathit{j}-\frac{1}{2}}\end{array}\right].$(65)

Based on the assumption of the perturbation parameter given in the governing problem ($0<\mathit{\epsilon },\mathit{\mu }\ll 1$), the stiffness matrix is non-singular. The algebraic system of 62(62) $\mathit{A}{U}^j=\mathit{B}.U^{\mathit{j}-1}+\mathit{C}(U^{\mathit{j}-1-\mathit{m}}+U^{\mathit{j}-\mathit{m}})+\mathit{f},$(62) implies the unknown vector values $U_{\mathit{i},\mathit{j}}$ and the finite element approximate solution, $\mathit{U}(\mathit{x},\mathit{t})$ is determined by

(66) $\mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}+1}\sum _{\mathit{i}=1}^{2\mathit{N}+1}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}){\psi }_j(\mathit{t}),$(66)

$\mathrm{\forall }(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}$, where ${\varphi }_i(\mathit{x})$ is space base function and ${\psi }_j(\mathit{t})$ is time base functions defined on (51(51) $\begin{array}{ll}{\varphi }_1^e(\mathit{x})& =\frac{1}{D}[x_{i-\frac{1}{2}}(x_i{)}^2-x_i(x_{i-\frac{1}{2}}{)}^2+((x_{i-\frac{1}{2}}{)}^2-(x_i{)}^2)\mathit{x}\\ +((x_i)-x_{i-\frac{1}{2}})x^2],\end{array}$(51) ) and (21(21) $\{{\psi }_j(\mathit{t}){\}}_1^M=\{\begin{array}{ll}\frac{\mathit{t}-{\mathit{t}}_{\mathit{j}-1}}{\mathrm{\Delta t}},& t_{\mathit{j}-1}\le \mathit{t}\le t_j,\\ \frac{{\mathit{t}}_{\mathit{j}+1}-\mathit{t}}{\mathrm{\Delta t}},& t_j\le \mathit{t}\le t_{\mathit{j}+1},\\ 0,& \mathit{otherwise}.\end{array}$(21) ).

4. Convergence analysis

The convergence of a finite element solution depends on the degree of approximate functions, the order of the derivative of $\mathit{u}(\mathit{x},\mathit{t})$ in the weak form, and the size of the mesh elements (Reddy, 2019). We used piecewise linear polynomials in time and quadratic Lagrange interpolation functions in space direction to approximate $\mathit{u}(\mathit{x},\mathit{t})$ throughout the domain. Hence, we are treating the h-convergence analysis, that is, the convergence by refining the size of the elements. Not increasing the degree of the approximate functions.

Errors introduced into finite element solutions are due to domain approximation, arithmetic, and approximation solutions. In our case, the domain is rectangular. Hence, the domain approximation error is not included (Surana & Reddy, 2016). The round-off error is reasonably small or zero in the desired decimal point accuracy. Then, the only inherent error introduced into the finite element solution is because of the approximation of the solution. The error estimate of the semi-discrete problem is discussed under Section 3 in maximum norm. Again, we apply the maximum norm to evaluate the error estimate in the spatial direction.

It is possible to decompose the numerical solution into $\mathit{U}=\mathit{V}+W_L+W_R,$ just as the continuous solution did in Section 3, decomposed into regular and singular components. The error, $\mathit{U}-\mathit{u}=(\mathit{V}-\mathit{v})+(W_L-w_L)+(W_R-w_R),$ and then triangle inequality gives its bound with

$\begin{array}{l}\Vert \mathit{U}-\mathit{u}\Vert \le \Vert \mathit{V}-\mathit{v}\Vert +\Vert W_L-w_L\Vert +\Vert W_R-w_R\Vert .\end{array}$

Now, take an arbitrary spatial element ${\overline{\mathcal{D}}}_{\mathit{i}}=[x_{\mathit{i}-1},x_i]$ for a fixed time level j.

The finite element approximation at $\mathit{t}=t_j$, $\mathit{u}(\mathit{x},t_j)$ with in this interval is

$\begin{array}{ll}{U}_j(\mathit{x})& =U_{\mathit{i}-1,\mathit{j}}{\varphi }_{\mathit{i}-1}(\mathit{x})+U_{\mathit{i},\mathit{j}}{\varphi }_i(\mathit{x})\\ |U_j(\mathit{x})|& \le \underset{{\mathcal{D}}_i}{max}(\mathit{u}(\mathit{x},t_j))({\varphi }_{\mathit{i}-1}(\mathit{x})+{\varphi }_i(\mathit{x})),\\ \text{since }{U}_{\mathit{i}-1,\mathit{j}},U_{\mathit{i},\mathit{j}}& \le \underset{{\mathcal{D}}_i}{max}(\mathit{u}(\mathit{x},t_j)),\mathrm{\forall x}\in [x_{\mathit{i}-1},x_i].\\ \text{Then }|U_j(\mathit{x})|& \le \underset{{\mathcal{D}}_i}{max}(|\mathit{u}(\mathit{x},t_j)|),\text{by the property }\\ {\varphi }_{\mathit{i}-1}(\mathit{x})+{\varphi }_i(\mathit{x})=1.\end{array}$

Theorem 4.1

Given a function $\mathit{u}(\mathit{x})\in C^2[0,1]$ and mesh points $x_o,x_1,x_2,\dots ,x_N$, the continuous piecewise quadratic function U(x) has an error estimate

(67) $\Vert \mathit{u}-\mathit{U}\Vert \le \mathit{C}{h}^2\Vert u^{″}\Vert .$(67)

Proof.

The error, $\mathit{e}=\mathit{u}(\mathit{x})-\mathit{U}(\mathit{x})$. At the mesh points, $\mathit{u}(x_i)=\mathit{U}(x_i)$ and $\mathit{u}(x_{\mathit{i}-1})=\mathit{U}(x_{\mathit{i}-1})$. Hence, $\mathit{e}(x_i)=\mathit{e}(x_{\mathit{i}-1})=0$. The Rolle’s theorem of calculus implies that there is $z_i\in (x_{\mathit{i}-1},x_i)$ such that $e^{\prime }(z_i)=0$. So, for any $\mathit{x}\in (x_{\mathit{i}-1},x_i),$ the Taylor expansion of $\mathit{e}(\mathit{x})=\mathit{e}(z_i+\mathit{x}-z_i)$ is

$\begin{array}{ll}\mathit{e}(\mathit{x})& =\mathit{e}(z_i)+e^{\prime }(z_i)(\mathit{x}-z_i)+\frac{e^{″}(z_i)}{2}(\mathit{x}-z_i{)}^2\\ +\frac{e^{‴}(z_i)}{3!}(\mathit{x}-z_i{)}^3+\dots \\ =\mathit{e}(z_i)+e^{\prime }(z_i)(\mathit{x}-z_i)+\frac{e^{″}(\mathit{\xi })}{2}(\mathit{x}-z_i{)}^2,\\ \text{for some }\mathit{\xi },x_{\mathit{i}-1}\le \mathit{\xi }\le x_i.\end{array}$

At $\mathit{x}=x_i$, ${\displaystyle \mathit{e}(x_i)=\mathit{e}(z_i)+\frac{1}{2}{e}^{″}(\mathit{\xi })(x_i-z_i{)}^2,}$ which implies $\mathit{e}(z_i)=-\frac{1}{2}{e}^{″}(\mathit{\xi })(x_i-z_i{)}^2.$

If ${\displaystyle x_{\mathit{i}-1}+\frac{h_i}{2}\le z_i\le x_i}$, then

$\begin{array}{ll}|\mathit{e}(z_i)|& \le \frac{1}{2}{\left(\frac{h_i}{2}\right)}^2|e^{″}(\mathit{\xi })|\\ \le \frac{1}{2}{\left(\frac{h_i}{2}\right)}^2|u^{″}(\mathit{\xi })-U^{″}(\mathit{\xi })|\\ =\frac{1}{2}{\left(\frac{h_i}{2}\right)}^2|u^{″}(\mathit{\xi })-\mathit{C}|\\ =\mathit{C}{h}_i^2\Vert u^{″}\Vert .\end{array}$

${\displaystyle \Vert \mathit{u}-\mathit{U}\Vert =\underset{1\le \mathit{i}\le \mathit{N}}{max}|\mathit{e}(z_i)|\le \mathit{C}{h}^2\Vert u^{″}\Vert .}$

If $x_{\mathit{i}-1}\le z_i\le x_i-\frac{h_i}{2},$ we arrive at the same result considering $\mathit{x}=x_{\mathit{i}-1}.$

In our case, Theorem 4.1 holds for $U_j(\mathit{x})$ and $\mathit{u}(\mathit{x},t_j)$ at any fixed time t_j. Using the second spatial derivative bounds of the continuous solution in Lemma 2.3, we can obtain the following error estimate.

(68) $|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\le \{\begin{array}{ll}\mathit{C}{h}_i^2{\epsilon }^{-1},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2},& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}$(68)

Applying the component solutions of $U_j(\mathit{x})$ and $\mathit{u}(\mathit{x},t_j)$,

$\begin{array}{ll}|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|& \le |V_j(\mathit{x})-\mathit{v}(\mathit{x},t_j)|+|W_{\mathit{Lj}}(\mathit{x})-w_L(\mathit{x},t_j)|\\ +|W_{\mathit{Rj}}(\mathit{x})-w_R(\mathit{x},t_j)|\\ \le |V_j(\mathit{x})-\mathit{v}(\mathit{x},t_j)|+|W_{\mathit{Lj}}(\mathit{x})|+|w_L(\mathit{x},t_j)|\\ +|W_{\mathit{Rj}}(\mathit{x})|+|w_R(\mathit{x},t_j)|\\ \le \mathit{C}{h}_i^2\underset{{\mathcal{D}}_i}{max}|v^{″}(\mathit{x},t_j)|+2\underset{{\mathcal{D}}_i}{max}|w_L(\mathit{x},t_j)|\\ +2\underset{{\mathcal{D}}_i}{max}|w_R(\mathit{x},t_j)|.\end{array}$

The bounds on Lemma 2.4 and Lemma 2.5 implies

(69) $\begin{array}{l}|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\\ \le \{\begin{array}{l}\mathit{C}(h_i^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}(h_i^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\varsigma }{\mu }(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}\end{array}$(69)

The error estimates in (68) and (69) depend on many factors. These error bounds are determined according to parameter values, whether $\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon }$ or $\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }$, the values of transition points ${\sigma }_1,{\sigma }_2\le 0.25$ or ${\sigma }_1,{\sigma }_2=0.25$, the position of the element, whether in the boundary regions or not and the propagation of the error in the delay term, when $\mathit{j}\le \mathit{m}$ or j > m, the way that error propagation varies. Before we get into the analysis of each case, let us state the following theorem.

Theorem 4.2

Let $U_j(\mathit{x})$ be the finite element solution of the semi-discrete problem 30(30) $\begin{array}{l}-\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_j(\mathit{x})-\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_j(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}+b_{j-\frac{1}{2}}(\mathit{x}))U_j(\mathit{x})\\ =\mathit{\epsilon }\frac{d^2}{\mathit{d}{x}^2}{U}_{\mathit{j}-1}(\mathit{x})+\mathit{\mu }{a}_{j-\frac{1}{2}}(\mathit{x})\frac{d}{\mathit{dx}}{U}_{\mathit{j}-1}(\mathit{x})+(\frac{2}{\mathrm{\Delta t}}-b_{j-\frac{1}{2}}(\mathit{x}))U_{\mathit{j}-1}(\mathit{x})\\ +c_{j-\frac{1}{2}}(\mathit{x})(U_{\mathit{j}-\mathit{m}}(\mathit{x})+U_{\mathit{j}-1-\mathit{m}}(\mathit{x}))+2f_{j-\frac{1}{2}}(\mathit{x}),\end{array}$(30) - 31(31) $\{\begin{array}{l}{U}_j(0)=0,\\ U_j(1)=0\end{array}\begin{array}{l}\mathit{Giventhat}{U}_{\mathit{j}-\mathit{m}}(\mathit{x})={\gamma }_b(\mathit{x},t_{\mathit{j}-\mathit{m}}),\\ \mathit{when}0\le \mathit{j}\le \mathit{m}.\end{array}$(31) on a fitted mesh ${\mathcal{D}}_x^N$, to the exact solution $\mathit{u}(\mathit{x},t_j)$, then

(70) $\Vert U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)\Vert \le \mathit{C}{N}^{-2}(\mathit{ln}\mathit{N}{)}^2.$(70)

Proof.

Assume $\Vert U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)\Vert \le \Vert U_{\mathit{j}-\mathit{m}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})\Vert +\mathit{k}(\mathit{N})$, for some mathematical expression k(N).

If $\mathit{j}\le \mathit{m},$ both $U_{\mathit{j}-\mathit{m}}(\mathit{x})=\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})={\gamma }_b(\mathit{x},t_{\mathit{j}-\mathit{m}})$ and $U_{\mathit{j}-\mathit{m}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})=0$. No error is encountered due to the delay term. Then $\Vert U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)\Vert \le 0+\mathit{k}(\mathit{N})$

if $\mathit{m}\le \mathit{j}\le 2\mathit{m}$ an additional error contributed due to the delay term $\Vert U_{\mathit{j}-\mathit{m}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})\Vert \le \mathit{k}(\mathit{N})$, then $\Vert U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)\Vert \le 2\mathit{k}(\mathit{N}).$

This can be proceed up to $(\mathit{k}-1)\mathit{m}\le \mathit{j}\le \mathit{km}=\mathit{T}.$ The current error bound is a constant multiple of the prior one. that is

(71) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}\Vert U_{\mathit{j}-\mathit{m}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}-\mathit{m}})\Vert .$(71)

Hence, it is sufficient to show the error estimate for $\mathit{j}\le \mathit{m}$ or for the time interval $[0,\mathit{\tau }]$. We focus the proof on the two major cases $\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon }$ and $\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }$,

Case I: $\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon }$

If ${\sigma }_1=\frac{1}{4}$, then there is uniform mesh length, $h_i=\frac{1}{N}$ and $\frac{1}{4}\le \frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\alpha \varsigma }}\mathit{ln}\mathit{N}$, which implies $\frac{1}{\epsilon }\le \mathit{C}{(\mathit{ln}\mathit{N})}^2.$ Substituting for h_i and $\frac{1}{\epsilon },$ in the 68(68) $|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\le \{\begin{array}{ll}\mathit{C}{h}_i^2{\epsilon }^{-1},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2},& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}$(68) , $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{h}_i^2{\epsilon }^{-1}$ gives

(72) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}{(\mathit{ln}\mathit{N})}^2.$(72)

If ${\sigma }_1<\frac{1}{4}$, then ${\sigma }_1={\sigma }_2=\frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\alpha \varsigma }}\mathit{ln}\mathit{N}.$ In the boundary regions , $1<\mathit{i}<\frac{N}{4}$ and $\frac{3\mathit{N}}{4}<\mathit{i}<\mathit{N}$, the mesh size $h_i=\frac{16\sqrt{\mathit{\epsilon }}}{\sqrt{\alpha \varsigma }}.\frac{\mathit{ln}\mathit{N}}{N}$ which implies $\frac{h_i}{\sqrt{\mathit{\epsilon }}}\le \frac{16}{\sqrt{\alpha \varsigma }}.\frac{\mathit{ln}\mathit{N}}{N}$. then again 68(68) $|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\le \{\begin{array}{ll}\mathit{C}{h}_i^2{\epsilon }^{-1},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2},& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}$(68) will be,

(73) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}{(\mathit{ln}\mathit{N})}^2.$(73)

In the outer boundary region, $\frac{N}{4}<\mathit{i}<\frac{3\mathit{N}}{4}$ for any i, ${\sigma }_1<x_i<{\sigma }_2$ and by 41(41) $\begin{array}{l}1<\mathit{i}<\mathit{N},h_i\le \frac{2}{N},\end{array}$(41) , $h_i\le \frac{2}{N}$ . We can have the following result from 69(69) $\begin{array}{l}|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\\ \le \{\begin{array}{l}\mathit{C}(h_i^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}(h_i^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\varsigma }{\mu }(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}\end{array}$(69) .

$\begin{array}{l}\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \\ \le \mathit{C}(h_i^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(1-x_{\mathit{i}})))\\ \le \mathit{C}((\frac{2}{N}{)}^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{\sigma }_1)+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}({\sigma }_2)))\\ \le \mathit{C}((\frac{2}{N}{)}^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}\frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\alpha \varsigma }}\mathit{ln}\mathit{N})\\ +\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(\frac{4\sqrt{\mathit{\epsilon }}}{\sqrt{\alpha \varsigma }}\mathit{ln}\mathit{N})))\\ \le \mathit{C}(N^{-2}+N^{-4}+N^{-4}).\end{array}$

Thus

(74) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}.$(74)

Case II: $\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }$

If ${\sigma }_1=\frac{1}{4}$, then $\frac{1}{4}<\frac{4\mathit{\epsilon }}{\mathit{\alpha \mu }}\mathit{ln}\mathit{N}$, which is $\frac{\mu }{\epsilon }\le \frac{16}{\alpha }\mathit{ln}\mathit{N}$. For any case of ${\sigma }_2=\frac{1}{4}$ or ${\sigma }_2<\frac{1}{4}$, the mesh size $h_i\le \frac{2}{N}.$ These relations with the second case of 68(68) $|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\le \{\begin{array}{ll}\mathit{C}{h}_i^2{\epsilon }^{-1},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2},& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}$(68) , $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2}$, gives

(75) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}{(\mathit{ln}\mathit{N})}^2.$(75)

If ${\sigma }_1<\frac{1}{4}$, then ${\sigma }_1=\frac{4\mathit{\epsilon }}{\mathit{\alpha \mu }}\mathit{ln}\mathit{N}.$ The second transition point ${\sigma }_2\le \frac{1}{4}.$

In the left boundary layer, $h_i=\frac{16\mathit{\epsilon }}{\mathit{\alpha \mu }}\frac{\mathit{ln}\mathit{N}}{N}.$ Inserting this value in 68(68) $|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\le \{\begin{array}{ll}\mathit{C}{h}_i^2{\epsilon }^{-1},& \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}{h}_i^2{\mu }^2{\epsilon }^{-2},& \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}$(68) , we get

(76) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}{(\mathit{ln}\mathit{N})}^2.$(76)

In the outer boundary layer region, ${\sigma }_1<x_i<{\sigma }_2,$ by 43(43) $\begin{array}{l}\frac{3\mathit{N}}{4}<\mathit{i}<\mathit{N},h_i\le \frac{4{\sigma }_2}{\mathit{N}}\le \frac{1}{N},\end{array}$(43) , $h_i\le \frac{2}{N}$ for all i, $\frac{n}{4}<\mathit{i}<\frac{3\mathit{N}}{4}.$ Applying the inequality $\frac{\epsilon }{\mu }\le \frac{16}{\varsigma }\mathit{ln}\mathit{N},$ (since $\frac{1}{\mu }\le \frac{16}{\varsigma }\mathit{ln}\mathit{N}$) the second case of 69(69) $\begin{array}{l}|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\\ \le \{\begin{array}{l}\mathit{C}(h_i^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}(h_i^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\varsigma }{\mu }(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}\end{array}$(69) , becomes

$\begin{array}{l}\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \\ \le \mathit{C}((\frac{2}{N}{)}^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{\sigma }_1)+\mathit{exp}(-\frac{\varsigma }{\mu }{\sigma }_2))\\ \le \mathit{C}(N^{-2}(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }\frac{4\mathit{\epsilon }}{\mathit{\alpha \mu }}\mathit{ln}\mathit{N})\\ +\mathit{exp}(-\frac{\varsigma }{\mu }\frac{4\mathit{\mu }}{\varsigma }\mathit{ln}\mathit{N}))\\ \le \mathit{C}(N^{-2}(1+\mathit{ln}\mathit{N})+\mathit{exp}(-4\mathit{ln}\mathit{N})+\mathit{exp}(-4\mathit{ln}\mathit{N}))\\ \le \mathit{C}(N^{-2}+N^{-2}\mathit{ln}\mathit{N}+N^{-4}+N^{-4}).\end{array}$

Thus,

(77) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}\mathit{ln}\mathit{N}.$(77)

In the right boundary layer region, ${\sigma }_1<{\sigma }_2<x_i<1$ and $h_i\le \frac{2}{N}$ for all i, $\frac{3\mathit{N}}{4}<\mathit{i}<\frac{N}{4}$. Taking the same 69(69) $\begin{array}{l}|U_j(\mathit{x})-\mathit{u}(\mathit{x},t_j)|\\ \le \{\begin{array}{l}\mathit{C}(h_i^2+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\sqrt{\alpha \varsigma }}{\sqrt{\epsilon }}(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon },\\ \mathit{C}(h_i^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{x}_{\mathit{i}-1})+\mathit{exp}(-\frac{\varsigma }{\mu }(1-x_{\mathit{i}}))),\\ \mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }.\end{array}\end{array}$(69) , we have

$\begin{array}{l}\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \\ \le \mathit{C}((\frac{2}{N}{)}^2(1+\mathit{\epsilon }{\mu }^{-1})+\mathit{exp}(-\frac{\mathit{\alpha \mu }}{\epsilon }{\sigma }_1)+\mathit{exp}(-\frac{\varsigma }{\mu }{\sigma }_2))\\ \le \mathit{C}(N^{-2}(1+\mathit{ln}\mathit{N})+\mathit{exp}(-4\mathit{ln}\mathit{N}))\\ \le \mathit{C}(N^{-2}+N^{-2}\mathit{ln}\mathit{N}+N^{-4}).\end{array}$

Then,

(78) $\Vert U_{\mathit{j}}(\mathit{x})-\mathit{u}(\mathit{x},t_{\mathit{j}})\Vert \le \mathit{C}{N}^{-2}\mathit{ln}\mathit{N}.$(78)

The next theorem completes our analysis with the temporal global error bond of lemma 3.3.

Theorem 4.3

Let $\mathit{u}(\mathit{x},\mathit{t})$ be the solution of problem (1(1) $\begin{array}{l}\frac{\partial }{\mathit{\partial t}}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\epsilon }\frac{{\partial }^2}{\mathit{\partial }{x}^2}\mathit{u}(\mathit{x},\mathit{t})-\mathit{\mu a}(\mathit{x},\mathit{t})\frac{\partial }{\mathit{\partial x}}\mathit{u}(\mathit{x},\mathit{t})+\mathit{b}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t})\\ =\mathit{c}(\mathit{x},\mathit{t})\mathit{u}(\mathit{x},\mathit{t}-\mathit{\tau })+\mathit{f}(\mathit{x},\mathit{t})\end{array}$(1) ) - (2(2) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})={\gamma }_b(\mathit{x},\mathit{t}),(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_b=[0,1]\times [-\mathit{\tau },0],\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_l=\{(0,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\\ \mathit{u}(\mathit{x},\mathit{t})=0,(\mathit{x},\mathit{t})\in {\mathrm{\Gamma }}_r=\{(1,\mathit{t}):-\mathit{\tau }\le \mathit{t}\le \mathit{T}\},\end{array}$(2) ) and $\mathit{U}(\mathit{x},\mathit{t})$ be the finite element approximate solution obtained by (46(46) $\mathit{u}(\mathit{x},\mathit{t})\approx \mathit{U}(\mathit{x},\mathit{t})=\sum _{\mathit{j}=1}^{\mathit{M}}\sum _{\mathit{i}=1}^{2\mathit{N}}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}){\psi }_j(\mathit{t}).$(46) ), then the following $\mathit{\epsilon },\mathit{\mu }-$ uniform error estimate holds,

(79) $\Vert \mathit{u}-\mathit{U}\Vert \le \mathit{C}(N^{-2}(\mathit{lnN}{)}^2+M^{-2})$(79)

where C is independent of ɛ and µ.

Proof.

The temporal and spatial error bounds in Lemma 3.3 and Theorem 4.2 implies the result.

5. Results and discussion

Two test problems are used for the numerical experiment. Extensive experiments are conducted on their numerical solutions. Some graphs are illustrated, and tables are presented for particular perturbation parameters and mesh numbers. The graphs and tables are used to show the maximum error estimates, convergence rates, and comparisons with previous study results.

A fine-grid solution is used to estimate errors when the exact solution to a problem is unknown. Thus, the double mesh principle is used to calculate the error estimate and rate of convergence of the proposed method. Let $\mathit{U}(\mathit{x},\mathit{t})$ on a grid mesh N × M and $\mathit{V}(\mathit{x},\mathit{t})$ on a doubled grid mesh $2\mathit{N}\times 2\mathit{M}$ be the FE solution for each perturbation parameter ɛ and µ. An element equation of a rectangular element $[x_{\mathit{i}},x_{\mathit{i}+1}]\times [t_{\mathit{j}},t_{\mathit{j}+1}]$ is given by

$\begin{array}{ll}{U}^e(\mathit{x},\mathit{t})& =\sum _{\mathit{p}=\mathit{j}}^{\mathit{j}+1}\sum _{\mathit{q}=2\mathit{i}-1}^{2\mathit{i}+1}{U}_{\mathit{p},\mathit{q}}.{\varphi }_p(\mathit{x}){\psi }_q(\mathit{t})\\ =U_{2\mathit{i}-1,\mathit{j}}.{\varphi }_1^e(\mathit{x}){\psi }_1^e(\mathit{t})+U_{2\mathit{i},\mathit{j}}.{\varphi }_2^e(\mathit{x}){\psi }_1^e(\mathit{t})\\ +U_{2\mathit{i}+1,\mathit{j}}.{\varphi }_3^e(\mathit{x}){\psi }_1^e(\mathit{t})+U_{2\mathit{i}-1,\mathit{j}+1}.{\varphi }_1^e(\mathit{x}){\psi }_2^e(\mathit{t})\\ +U_{2\mathit{i},\mathit{j}+1}.{\varphi }_2^e(\mathit{x}){\psi }_2^e(\mathit{t})+U_{2\mathit{i}+1,\mathit{j}+1}.{\varphi }_3^e(\mathit{x}){\psi }_2^e(\mathit{t})\\ =(U_{2\mathit{i}-1,\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i}-1,\mathit{j}+1}.{\psi }_2^e(\mathit{t})){\varphi }_1^e(\mathit{x})\\ +(U_{2\mathit{i},\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i},\mathit{j}+1}.{\psi }_2^e(\mathit{t})){\varphi }_2^e(\mathit{x})\\ +(U_{2\mathit{i}+1,\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i}+1,\mathit{j}+1}.{\psi }_2^e(\mathit{t})){\varphi }_3^e(\mathit{x})\end{array}$

The hat functions, ψ_j and Lagrange quadratic interpolation functions, ϕ_i referred by (21(21) $\{{\psi }_j(\mathit{t}){\}}_1^M=\{\begin{array}{ll}\frac{\mathit{t}-{\mathit{t}}_{\mathit{j}-1}}{\mathrm{\Delta t}},& t_{\mathit{j}-1}\le \mathit{t}\le t_j,\\ \frac{{\mathit{t}}_{\mathit{j}+1}-\mathit{t}}{\mathrm{\Delta t}},& t_j\le \mathit{t}\le t_{\mathit{j}+1},\\ 0,& \mathit{otherwise}.\end{array}$(21) ) and (51), in this specified element for $\mathit{j}=1,2$ and $\mathit{i}=1,2,3$ are determined by

(80) $\begin{array}{l}\{\begin{array}{l}{\psi }_1^e=\frac{1}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})\\ {\psi }_2^e=\frac{1}{\mathrm{\Delta t}}(\mathit{t}-t_{\mathit{j}})\end{array}\text{and }\\ \{\begin{array}{l}{\varphi }_1^e=\frac{1}{2}{x}_i^2{h}_i+\frac{3}{4}{x}_i{h}_i^2+\frac{1}{4}{h}_i^3-(x_i{h}_i+\frac{3}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2\\ {\varphi }_2^e=-x_i^2{h}_i-x_i{h}_i^2+(2x_i{h}_i+h_i^2)\mathit{x}-h_i{x}^2\\ {\varphi }_3^e=\frac{1}{2}{x}_i^2{h}_i+\frac{1}{4}{x}_i{h}_i^2-(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2\end{array}\end{array}$(80)

$\begin{array}{ll}{U}^e(\mathit{x},\mathit{t})& =(U_{2\mathit{i}-1,\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i}-1,\mathit{j}+1}.{\psi }_2^e(\mathit{t}))\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{3}{4}{x}_i{h}_i^2+\frac{1}{4}{h}_i^3-(x_i{h}_i+\frac{3}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2)\\ +(U_{2\mathit{i},\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i},\mathit{j}+1}.{\psi }_2^e(\mathit{t}))\\ \times (-x_i^2{h}_i-x_i{h}_i^2+(2x_i{h}_i+h_i^2)\mathit{x}-h_i{x}^2)\\ +(U_{2\mathit{i}+1,\mathit{j}}.{\psi }_1^e(\mathit{t})+U_{2\mathit{i}+1,\mathit{j}+1}.{\psi }_2^e(\mathit{t}))\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{1}{4}{x}_i{h}_i^2-(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2)\\ =U_{2\mathit{i}-1,\mathit{j}}.\frac{1}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{3}{4}{x}_i{h}_i^2+\frac{1}{4}{h}_i^3-(x_i{h}_i+\frac{3}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2)\\ +U_{2\mathit{i},\mathit{j}}.\frac{1}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})\\ \times (-x_i^2{h}_i-x_i{h}_i^2+(2x_i{h}_i+h_i^2)\mathit{x}-h_i{x}^2)\\ +U_{2\mathit{i}+1,\mathit{j}}.\frac{1}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{1}{4}{x}_i{h}_i^2-(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2)\\ +U_{2\mathit{i}-1,\mathit{j}+1}.\frac{1}{\mathrm{\Delta t}}(\mathit{t}-t_{\mathit{j}})\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{3}{4}{x}_i{h}_i^2+\frac{1}{4}{h}_i^3-(x_i{h}_i+\frac{3}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2)\\ +U_{2\mathit{i},\mathit{j}+1}.\frac{1}{\mathrm{\Delta t}}(\mathit{t}-t_{\mathit{j}})\\ \times (-x_i^2{h}_i-x_i{h}_i^2+(2x_i{h}_i+h_i^2)\mathit{x}-h_i{x}^2)\\ +U_{2\mathit{i}+1,\mathit{j}+1}.\frac{1}{\mathrm{\Delta t}}(\mathit{t}-t_{\mathit{j}})\\ \times (\frac{1}{2}{x}_i^2{h}_i+\frac{1}{4}{x}_i{h}_i^2-(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}+\frac{1}{2}{h}_i{x}^2).\end{array}$

Now, dividing each mesh in to two equal length (refer Figures 1 and 2), similar numerical solutions, ${\displaystyle V_{11}^e(\mathit{x},\mathit{t}),V_{12}^e(\mathit{x},\mathit{t}),V_{21}^e(\mathit{x},\mathit{t})}$ and ${\displaystyle V_{22}^e(\mathit{x},\mathit{t})}$ are derived for the four quarter regions of the element. For each corresponding region, ${\displaystyle [x_i,x_{i+\frac{1}{2}}]\times [t_j,t_{j+\frac{1}{2}}]}$, ${\displaystyle [x_{i+\frac{1}{2}},x_{\mathit{i}+1}]\times [t_j,t_{j+\frac{1}{2}}]}$, ${\displaystyle [x_i,x_{i+\frac{1}{2}}]\times [t_{j+\frac{1}{2}},t_{\mathit{j}+1}]}$, and ${\displaystyle [x_{i+\frac{1}{2}},x_{\mathit{i}+1}]\times [t_{j\frac{1}{2}},t_{\mathit{j}+1}]}$, the numerical solution is indicated below.

Figure 1. Global numerical solution indices refers to local nodes of arbitrary finite element $[x_{\mathit{i}},x_{\mathit{i}+1}]\times [t_{\mathit{j}},t_{\mathit{j}+1}]$.

Figure 2. Global numerical solution indices refers to local nodes of arbitrary finite element $[x_{\mathit{i}},x_{\mathit{i}+1}]\times [t_{\mathit{j}},t_{\mathit{j}+1}]$ after double split of the spatial and temporal finite elements.

If $(\mathit{x},\mathit{t})\in [x_i,x_{i+\frac{1}{2}}]\times [t_j,t_{j+\frac{1}{2}}],$ then

$\begin{array}{ll}{V}_{11}^e(\mathit{x},\mathit{t})& =V_{4\mathit{i}-3,2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{3}{16}{x}_i{h}_i^2\\ +\frac{1}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{3}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-2,2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})\\ \times (-\frac{1}{2}{x}_i^2{h}_i-\frac{1}{4}{x}_i{h}_i^2+(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})\\ \times (\frac{1}{4}{x}_i^2{h}_i+\frac{1}{16}{x}_i{h}_i^2-(\frac{1}{2}{x}_i{h}_i+\frac{1}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-3,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)(\frac{1}{4}{x}_i^2{h}_i+\frac{3}{16}{x}_i{h}_i^2\\ +\frac{1}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{3}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-2,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)\\ \times (-\frac{1}{2}{x}_i^2{h}_i-\frac{1}{4}{x}_i{h}_i^2+(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)\\ \times (\frac{1}{4}{x}_i^2{h}_i+\frac{1}{16}{x}_i{h}_i^2-(\frac{1}{2}{x}_i{h}_i+\frac{1}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2).\end{array}$

If $(\mathit{x},\mathit{t})\in [x_{i+\frac{1}{2}},x_{\mathit{i}+1}]\times [t_j,t_{j+\frac{1}{2}}],$ then

$\begin{array}{ll}{V}_{12}^e(\mathit{x},\mathit{t})& =V_{4\mathit{i}-1,2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{7}{16}{x}_i{h}_i^2\\ +\frac{3}{16}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{7}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i},2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})(-\frac{1}{2}{x}_i^2{h}_i-\frac{3}{4}{x}_i{h}_i^2-\frac{1}{4}{h}_i^3\\ +(x_i{h}_i+\frac{5}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}+1,2\mathit{j}-1}.\frac{2}{\mathrm{\Delta t}}(t_{j+\frac{1}{2}}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{5}{16}{x}_i{h}_i^2\\ +\frac{3}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{5}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)(\frac{1}{4}{x}_i^2{h}_i+\frac{7}{16}{x}_i{h}_i^2\\ +\frac{3}{16}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{7}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i},2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)(-\frac{1}{2}{x}_i^2{h}_i-\frac{3}{4}{x}_i{h}_i^2-\frac{1}{4}{h}_i^3\\ +(x_i{h}_i+\frac{5}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}+1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_j)(\frac{1}{4}{x}_i^2{h}_i+\frac{5}{16}{x}_i{h}_i^2\\ +\frac{3}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{5}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2).\end{array}$

If $(\mathit{x},\mathit{t})\in [x_i,x_{i+\frac{1}{2}}]\times [t_{j+\frac{1}{2}},t_{\mathit{j}+1}],$ then

$\begin{array}{ll}{V}_{21}^e(\mathit{x},\mathit{t})& =V_{4\mathit{i}-3,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i\\ +\frac{3}{16}{x}_i{h}_i^2+\frac{1}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{3}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-2,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(-\frac{1}{2}{x}_i^2{h}_i\\ -\frac{1}{4}{x}_i{h}_i^2+(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{1}{16}{x}_i{h}_i^2\\ -(\frac{1}{2}{x}_i{h}_i+\frac{1}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-3,2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(\frac{1}{4}{x}_i^2{h}_i+\frac{3}{16}{x}_i{h}_i^2\\ +\frac{1}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{3}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-2,2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(-\frac{1}{2}{x}_i^2{h}_i\\ -\frac{1}{4}{x}_i{h}_i^2+(x_i{h}_i+\frac{1}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(\frac{1}{4}{x}_i^2{h}_i\\ +\frac{1}{16}{x}_i{h}_i^2-(\frac{1}{2}{x}_i{h}_i+\frac{1}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2).\end{array}$

If $(\mathit{x},\mathit{t})\in [x_{i+\frac{1}{2}},x_{\mathit{i}+1}]\times [t_{j+\frac{1}{2}},t_{\mathit{j}+1}],$ then

$\begin{array}{ll}{V}_{22}^e(\mathit{x},\mathit{t})& =V_{4\mathit{i}-1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{7}{16}{x}_i{h}_i^2\\ +\frac{3}{16}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{7}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i},2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(-\frac{1}{2}{x}_i^2{h}_i-\frac{3}{4}{x}_i{h}_i^2\\ -\frac{1}{4}{h}_i^3+(x_i{h}_i+\frac{5}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}+1,2\mathit{j}}.\frac{2}{\mathrm{\Delta t}}(t_{\mathit{j}+1}-\mathit{t})(\frac{1}{4}{x}_i^2{h}_i+\frac{5}{16}{x}_i{h}_i^2\\ +\frac{3}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{5}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i}-1,2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(\frac{1}{4}{x}_i^2{h}_i+\frac{7}{16}{x}_i{h}_i^2\\ +\frac{3}{16}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{7}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2)\\ +V_{4\mathit{i},2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(-\frac{1}{2}{x}_i^2{h}_i-\frac{3}{4}{x}_i{h}_i^2\\ -\frac{1}{4}{h}_i^3+(x_i{h}_i+\frac{5}{4}{h}_i^2)\mathit{x}-\frac{1}{2}{h}_i{x}^2)\\ +V_{4\mathit{i}+1,2\mathit{j}+1}.\frac{2}{\mathrm{\Delta t}}(\mathit{t}-t_{j+\frac{1}{2}})(\frac{1}{4}{x}_i^2{h}_i+\frac{5}{16}{x}_i{h}_i^2\\ +\frac{3}{32}{h}_i^3-(\frac{1}{2}{x}_i{h}_i+\frac{5}{16}{h}_i^2)\mathit{x}+\frac{1}{4}{h}_i{x}^2).\end{array}$

Let an element numerical solution using double mesh be

(81) $V^e(\mathit{x},\mathit{t})=V_{11}^e(\mathit{x},\mathit{t})+V_{12}^e(\mathit{x},\mathit{t})+V_{21}^e(\mathit{x},\mathit{t})+V_{22}^e(\mathit{x},\mathit{t}),$(81)

then, the error committed at any $(\mathit{x},\mathit{t})\in [x_{\mathit{i}},x_{\mathit{i}+1}]\times [t_{\mathit{j}},t_{\mathit{j}+1}]$ will be

(82) $\mathit{e}(\mathit{x},\mathit{t})=U^e(\mathit{x},\mathit{t})-V^e(\mathit{x},\mathit{t}).$(82)

The error $\mathit{e}(\mathit{x},\mathit{t})$ is linear in t and quadratic in x in each quarter region. ${\displaystyle \frac{\partial }{\mathit{\partial t}}\mathit{e}(\mathit{x},\mathit{t})}$ depends on x only, which means it is constant at a fixed value of $\mathit{x}.$ Hence, the maximum value of $|\mathit{e}(\mathit{x},\mathit{t})|$ in the finite element attain at $\mathit{t}=t_j,t_{j+\frac{1}{2}}\text{or }{t}_{\mathit{j}+1}.$ Let ${\displaystyle \frac{\partial }{\mathit{\partial x}}\mathit{e}(\mathit{x},\mathit{t})=0}$ at ${\overline{x}}_1\in (x_i,x_{i+\frac{1}{2}})$ and ${\overline{x}}_2\in (x_{i+\frac{1}{2}},x_{\mathit{i}+1}).$ Therefore, the maximum absolute error in the local finite element exists at either of the 15 critical points $(x_{\mathit{i}},t_{\mathit{j}})$, $({\overline{x}}_1,t_{\mathit{j}})$, $(x_{i+\frac{1}{2}},t_{\mathit{j}})$, $({\overline{x}}_2,t_{\mathit{j}})$, $(x_{\mathit{i}+1},t_{\mathit{j}})$, $(x_{\mathit{i}},t_{j+\frac{1}{2}})$, $({\overline{x}}_1,t_{j+\frac{1}{2}})$, $(x_{i+\frac{1}{2}},t_{j+\frac{1}{2}})$, $({\overline{x}}_2,t_{j+\frac{1}{2}})$, $(x_{\mathit{i}+1},t_{j+\frac{1}{2}})$, $(x_{\mathit{i}},t_{\mathit{j}+1})$, $({\overline{x}}_1,t_{\mathit{j}+1})$, $(x_{i+\frac{1}{2}},t_{\mathit{j}+1})$, or $({\overline{x}}_2,t_{\mathit{j}+1}),(x_{\mathit{i}+1},t_{\mathit{j}+1}).$

So, the element-wise maximum absolute error $e_{\epsilon .\mu }^{N,M}$ is given by

$\begin{array}{ll}{e}_{\epsilon .\mu }^{N,M}& =max(|\mathit{e}(x_{\mathit{i}},t_{\mathit{j}})|,|\mathit{e}({\overline{x}}_1,t_{\mathit{j}})|,|\mathit{e}(x_{i+\frac{1}{2}},t_{\mathit{j}})|,|\mathit{e}({\overline{x}}_2,t_{\mathit{j}})|,\\ |\mathit{e}(x_{\mathit{i}+1},t_{\mathit{j}})|,|\mathit{e}(x_{\mathit{i}},t_{j+\frac{1}{2}})|,|\mathit{e}({\overline{x}}_1,t_{j+\frac{1}{2}})|,|\mathit{e}(x_{i+\frac{1}{2}},t_{j+\frac{1}{2}})|,\\ |\mathit{e}({\overline{x}}_2,t_{j+\frac{1}{2}})|,|\mathit{e}(x_{\mathit{i}+1},t_{j+\frac{1}{2}})|,|\mathit{e}(x_{\mathit{i}},t_{\mathit{j}+1})|,|\mathit{e}({\overline{x}}_1,t_{\mathit{j}+1})|,\\ |\mathit{e}(x_{i+\frac{1}{2}},t_{\mathit{j}+1})|,|\mathit{e}({\overline{x}}_2,t_{\mathit{j}+1})|,|\mathit{e}(x_{\mathit{i}+1},t_{\mathit{j}+1})|).\end{array}$

Because of the continuity and balance properties of the finite element solution at junction edges and error at the initial and boundary of the domain is zero, we can consider only eight of the critical points to compute the global maximum absolute error during assembly. Then, the global maximum absolute error for each ɛ and µ $(E_{\epsilon ,\mu }^{N,M})$ of the numerical solution is

$\begin{array}{ll}{E}_{\epsilon ,\mu }^{N,M}& =\underset{(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}}{sup}|\mathit{U}(\mathit{x},\mathit{t})-\mathit{V}(\mathit{x},\mathit{t})|\\ =\underset{(\mathit{x},\mathit{t})\in \overline{\mathcal{D}}}{sup}|\sum _{\mathit{j}=1}^{\mathit{M}}\sum _{\mathit{i}=1}^{2\mathit{N}}{U}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}){\psi }_j(\mathit{t})\\ -\sum _{\mathit{j}=1}^{2\mathit{M}}\sum _{\mathit{i}=1}^{4\mathit{N}}{V}_{\mathit{i},\mathit{j}}.{\varphi }_i(\mathit{x}){\psi }_j(\mathit{t})|\\ =\underset{\begin{array}{l}1\le \mathit{i}\le \mathit{N}\\ 1\le \mathit{j}\le \mathit{M}\end{array}}{max}\{|\mathit{e}(x_{\mathit{i}},t_{j+\frac{1}{2}})|,|\mathit{e}({\overline{x}}_1,t_{j+\frac{1}{2}})|,|\mathit{e}(x_{i+\frac{1}{2}},t_{j+\frac{1}{2}})|,\\ |\mathit{e}({\overline{x}}_2,t_{j+\frac{1}{2}})|,|\mathit{e}(x_{\mathit{i}},t_{\mathit{j}+1})|,|\mathit{e}({\overline{x}}_1,t_{\mathit{j}+1})|,\\ |\mathit{e}(x_{i+\frac{1}{2}},t_{\mathit{j}+1})|,|\mathit{e}({\overline{x}}_2,t_{\mathit{j}+1})|\}.\end{array}$

The corresponding rate of convergence,$R_{\epsilon ,\mu }^{N,M}$ is given by

(83) $R_{\epsilon ,\mu }^{N,M}={\log }_2\left(E_{\epsilon ,\mu }^{N,M}\right)-{\log }_2\left(E_{\epsilon ,\mu }^{2N,2M}\right).$(83)

The maximum absolute error, $E^{\mathit{N},\mathit{M}}$ and the corresponding rate of convergence, $R^{\mathit{N},\mathit{M}}$ for any perturbation parameters ɛ and µ mutually near to zero is determined by

(84) $E^{\mathit{N},\mathit{M}}=\underset{0<\mathit{\epsilon },\mathit{\mu }\ll 1}{max}\left(E_{\epsilon ,\mu }^{N,M}\right)\text{and }$(84)

(85) $R^{\mathit{N},\mathit{M}}={\log }_2\left(E^{\mathit{N},\mathit{M}}\right)-{\log }_2\left(E^{2\mathit{N},2\mathit{M}}\right)$(85)

The following two test examples are adopted from Negero and Duressa (2021) for numerical experiment.

Example 5.1.

(86) $\begin{array}{l}-\mathit{\epsilon }\frac{{\partial }^2\mathit{u}}{\mathit{\partial }{x}^2}(\mathit{x},\mathit{t})-\mathit{\mu }(\mathit{x}+1)\frac{\mathit{\partial u}}{\mathit{\partial x}}(\mathit{x},\mathit{t})+\mathit{u}(\mathit{x},\mathit{t})+\frac{\mathit{\partial u}}{\mathit{\partial t}}\\ =\mathit{u}(\mathit{x},\mathit{t}-1)-16x^2(1-\mathit{x}{)}^2,\end{array}$(86)

$(\mathit{x},\mathit{t})\in (0,1)\times (0,2],$ subjected to the initial and boundary conditions

(87) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})=0,\mathrm{\forall }(\mathit{x},\mathit{t})\in [0,1]\times [-1,0]\\ \mathit{u}(0,\mathit{t})=0,\mathrm{\forall t}\in [0,2],\\ \mathit{u}(1,\mathit{t})=0,\mathrm{\forall t}\in [0,2].\end{array}$(87)

Example 5.2.

$\begin{array}{l}-\mathit{\epsilon }\frac{{\partial }^2\mathit{u}}{\mathit{\partial }{x}^2}(\mathit{x},\mathit{t})-\mathit{\mu }(1+\mathit{x}(1-\mathit{x})+t^2)\frac{\mathit{\partial u}}{\mathit{\partial x}}(\mathit{x},\mathit{t})+(1+5\mathit{xt})\mathit{u}(\mathit{x},\mathit{t})+\frac{\mathit{\partial u}}{\mathit{\partial t}}\\ =\mathit{u}(\mathit{x},\mathit{t}-1)-\mathit{x}(1-\mathit{x})(e^x-1),(\mathit{x},\mathit{t})\in (0,1)\times (0,2]\end{array}$

subjected to the initial and boundary conditions

(88) $\{\begin{array}{l}\mathit{u}(\mathit{x},\mathit{t})=0,\mathrm{\forall }(\mathit{x},\mathit{t})\in [0,1]\times [-1,0]\\ \mathit{u}(0,\mathit{t})=0,\mathrm{\forall t}\in [0,2],\\ \mathit{u}(1,\mathit{t})=0,\mathrm{\forall }\mathit{t}\in [0,2].\end{array}$(88)

The governing equation’s assumptions and the necessary compatibility conditions in 4(4) $\begin{array}{ll}{\gamma }_b(0,0)& =0,\end{array}$(4) -(7(7) $\begin{array}{l}{\mathcal{L}}_{\mathit{\epsilon },\mathit{\mu }}(\mathit{u}(1,0))=\mathit{c}(1,0)\mathit{u}(1,-\mathit{\tau })+\mathit{f}(1,0).\end{array}$(7) ) are met. We are confident in the existence and uniqueness of the solutions for Example 5.1 and Example 5.2. However, their exact solutions are unknown. Thus, the double mesh concept is used to calculate the element-wise error estimate and rate of convergence of the proposed method. The numerical experiment is performed based on the convenience perturbation parameters and the number of finite elements so that $0<\mathit{\epsilon },\mathit{\mu }\ll 1$ and according to the performance of our personal computer. Figure 3 is the finite element solution of Example 5.1 for a $32\times 32$ rectangular elements grid and its double mesh grid. For this case, element wise absolute maximum error is illustrated in Figure 4. The $64\times 64$ rectangular elements grid and its double mesh grid numerical solutions with their corresponding error graphs are shown in Figures 5 and 6. Similar investigation is done for Example 5.2 as shown in Figures 7 to 10. Through numerical experimentation, it has been found that the point-wise (as well as the maximum) absolute error reduces as the number of finite elements increases.

Figure 3. Graphs of the FE solution of example 5.1 for $32\times 32$ mesh and double mesh grid with parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}.$.

Figure 4. The graph of element-wise absolute maximum error of the FE solution in example 5.1 for $32\times 32$ mesh when the parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}.$.

Figure 5. The graph of element-wise absolute maximum error of the FE solution in example 5.1 for $64\times 64$ mesh when the parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}.$.

Figure 6. The graph of element-wise absolute maximum error of the FE solution in example 5.1 for $64\times 64$ mesh when the parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}.$.

Figure 7. Graphs of the FE solution of example 5.2 for $32\times 32$ mesh and double mesh grid with parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}.$.

Figure 8. The graph of element-wise absolute maximum error of the FE solution in example 5.2 for $32\times 32$ mesh when the parameters $\mathit{\epsilon }={10}^{-10}$ and $\mathit{\mu }={10}^{-20}$.

Figure 9. Graphs of the FE solution of example 5.2 for 64 × 64 mesh and double mesh grid with parameters ε = 10 − 10 and μ = 10 − 20.

Figure 10. The graph of element-wise absolute maximum error of the FE solution in example 5.2 for 64 × 64 mesh when the parameters ε = 10 − 10 and μ = 10 − 20.

The error estimate and rate of convergence of the finite element solution of the first test problem in Example 5.1 are tabulated in Tables (1) and (2). That is performed by fixing one of the perturbation parameters while reducing the other so that addressing the two cases $\mathit{\alpha }{\mu }^2\le \mathit{\varsigma \epsilon }$ and $\mathit{\alpha }{\mu }^2\ge \mathit{\varsigma \epsilon }$. Table (1) shows the proposed method is uniformly convergent independent of the perturbation parameter µ, and Table (2) shows it is again uniformly convergent independent of the perturbation parameter ɛ. A similar presentation is given in Table (3) and Table (4) for the second problem of Example 5.2. The maximum absolute error decreases uniformly as the number of meshes increases, irrespective of perturbation parameters, with almost a second order of convergence. So, this experiment verified that the method is $\mathit{\epsilon },\mathit{\mu }-$ uniform.

Table 1. Error estimate and convergence rate of example 5.1, for $\mathit{\epsilon }={10}^{-10}$

$\mathit{\mu }\downarrow$		N = 16	32	64	128	256	512
		M = 16	32	64	128	256	512
10^−6	$E_{\epsilon ,\mu }^{N,M}$	6.6767e-03	1.7327e-03	4.4312e-04	1.1309e-04	2.8711e-05	7.2680e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9461	1.9673	1.9702	1.9778	1.9819
10^−8	$E_{\epsilon ,\mu }^{N,M}$	6.6769e-03	1.7328e-03	4.4313e-04	1.1310e-04	2.8712e-05	7.2681e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9461	1.9673	1.9701	1.9779	1.9820
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
10^−18	$E_{\epsilon ,\mu }^{N,M}$	6.6769e-03	1.7328e-03	4.4313e-04	1.1310e-04	2.8712e-05	7.2681e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9461	1.9673	1.9701	1.9779	1.9820
10^−20	$E_{\epsilon ,\mu }^{N,M}$	6.6769e-03	1.7328e-03	4.4313e-04	1.1310e-04	2.8712e-05	7.2681e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9461	1.9673	1.9701	1.9779	1.9820
	$E^{\mathit{N},\mathit{M}}$	6.6769e-03	1.7328e-03	4.4313e-04	1.1310e-04	2.8712e-05	7.2681e-06
	$R^{\mathit{N},\mathit{M}}$	1.9461	1.9673	1.9701	1.9779	1.9820

Table 2. Error estimate and convergence rate of example 5.1, for $\mathit{\mu }={10}^{-10}$

$\mathit{\epsilon }\downarrow$		N = 16	32	64	128	256	512
		M = 16	32	64	128	256	512
10^−18	$E_{\epsilon ,\mu }^{N,M}$	6.6831e-03	1.7347e-03	4.4377e-04	1.1335e-04	2.8789e-05	7.2902e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9458	1.9668	1.9690	1.9772	1.9815
10^−20	$E_{\epsilon ,\mu }^{N,M}$	6.6831e-03	1.7347e-03	4.4377e-04	1.1335e-04	2.8789e-05	7.2902e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9458	1.9668	1.9690	1.9772	1.9815
⋮		⋮	⋮	⋮	⋮	⋮	⋮
10^−28	$E_{\epsilon ,\mu }^{N,M}$	6.6833e-03	1.7347e-03	4.4377e-04	1.1335e-04	2.8789e-05	7.2902e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9459	1.9668	1.9690	1.9772	1.9815
10^−30	$E_{\epsilon ,\mu }^{N,M}$	6.6833e-03	1.7347e-03	4.4377e-04	1.1335e-04	2.8789e-05	7.2902e-06
	$R_{\epsilon ,\mu }^{N,M}$	1.9459	1.9668	1.9690	1.9772	1.9815
	$E^{\mathit{N},\mathit{M}}$	6.6833e-03	1.7347e-03	4.4377e-04	1.1335e-04	2.8789e-05	7.2902e-06
	$R^{\mathit{N},\mathit{M}}$	1.9459	1.9668	1.9690	1.9772	1.9815

Table 3. Error estimate and convergence rate of example 5.2, for $\mathit{\epsilon }={10}^{-10}$

$\mathit{\mu }\downarrow$		N = 16	32	64	128	256	512
		M = 16	32	64	128	256	512
10^−6	$E_{\epsilon ,\mu }^{N,M}$	5.3629e-04	1.3641e-04	3.4453e-05	8.6810e-06	2.1805e-06	5.4584e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9751	1.9852	1.9887	1.9932	1.9981
10^−8	$E_{\epsilon ,\mu }^{N,M}$	5.3630e-04	1.3642e-04	3.4454e-05	8.6812e-06	2.1806e-06	5.4586e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9750	1.9853	1.9887	1.9932	1.9981
⋮		⋮	⋮	⋮	⋮	⋮	⋮
10^−18	$R_{\epsilon ,\mu }^{N,M}$	5.3630e-04	1.3642e-04	3.4454e-05	8.6812e-06	2.1806e-06	5.4586e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9750	1.9853	1.9887	1.9932	1.9981
10^−20	$E_{\epsilon ,\mu }^{N,M}$	5.3630e-04	1.3642e-04	3.4454e-05	8.6812e-06	2.1806e-06	5.4586e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9750	1.9853	1.9887	1.9932	1.9981
	$E^{\mathit{N},\mathit{M}}$	5.3630e-04	1.3642e-04	3.4454e-05	8.6812e-06	2.1806e-06	5.4586e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9750	1.9853	1.9887	1.9932	1.9981

Table 4. Error estimate and convergence rate of example 5.2, for $\mathit{\mu }={10}^{-10}$

$\mathit{\epsilon }\downarrow$		N = 16	32	64	128	256	512
		M = 16	32	64	128	256	512
10^−18	$E_{\epsilon ,\mu }^{N,M}$	5.3713e-04	1.3668e-04	3.4621e-05	8.6793e-06	2.1722e-06	5.4316e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9745	1.9811	1.9960	1.9984	1.9997
10^−20	$E_{\epsilon ,\mu }^{N,M}$	5.3713e-04	1.3668e-04	3.4621e-05	8.6793e-06	2.1722e-06	5.4316e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9745	1.9811	1.9960	1.9984	1.9997
⋮		⋮	⋮	⋮	⋮	⋮	⋮
10^−28	$E_{\epsilon ,\mu }^{N,M}$	5.3713e-04	1.3668e-04	3.4621e-05	8.6793e-06	2.1722e-06	5.4316e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9745	1.9811	1.9960	1.9984	1.9997
10^−30	$E_{\epsilon ,\mu }^{N,M}$	5.3713e-04	1.3668e-04	3.4621e-05	8.6793e-06	2.1722e-06	5.4316e-07
	$R_{\epsilon ,\mu }^{N,M}$	1.9745	1.9811	1.9960	1.9984	1.9997
	$E^{\mathit{N},\mathit{M}}$	5.3713e-04	1.3668e-04	3.4621e-05	8.6793e-06	2.1722e-06	5.4316e-07
	$R^{\mathit{N},\mathit{M}}$	1.9745	1.9811	1.9960	1.9984	1.9997

The logarithmic scaled graph of maximum absolute error of the finite element solutions of the two test problems for some particular parameters are given in Figure 11. The graphs are almost linear with slope of − 2 which strengthen the predicted convergence rate.

Figure 11. logarithmic scaled graphs of maximum error versus number of meshes:(a)Example 5.1, and (b) example 5.2.

Assuming $\mathit{\epsilon }={10}^{-4}$ and $\mathit{\mu }={10}^{-9}$, Ayele et al. (2023) implements the $\mathit{\theta }-$ approach to solve Problem 5.1. For both θ = 1 and $\mathit{\theta }=\frac{1}{2}$, the convergence rate and error are compared with the numerical solution provided in Kumar, et al. (2020). Our numerical results and the findings from earlier studies are included in Table 5. As a result, our method offers a smaller error and a better rate of convergence for the same perturbation parameters and domain discretization.

Table 5. Comparison of the error estimate $(E^{\mathit{N},\mathit{M}})$ and rate of convergence $(R^{\mathit{N},\mathit{M}})$ of example 5.1, for $\mathit{\epsilon }={10}^{-4}$ and $\mathit{\mu }={10}^{-9}$

	$\mathit{N}=$	32	64	128	256	512
Method	$\mathit{M}=$	8	16	32	64	128
(Kumar, et al., 2020)	$E^{\mathit{N},\mathit{M}}$	1.1100e-02	5.0838e-03	2.4640e-03	1.2162e-03	6.0457e-04
	$R^{\mathit{N},\mathit{M}}$	1.1265	1.0449	1.0185	1.0084
(Ayele et al., 2023),θ = 1	$E^{\mathit{N},\mathit{M}}$	5.4545e-03	2.8706e-03	1.4692e-03	7.4328e-04	3.7384e-04
	$R^{\mathit{N},\mathit{M}}$	0.9261	0.9663	0.9830	0.9915
(Ayele et al., 2023),$\mathit{\theta }=\frac{1}{2}$	$E^{\mathit{N},\mathit{M}}$	1.3176e-03	3.2988e-04	8.2475e-05	2.0620e-05	5.1550e-06
	$R^{\mathit{N},\mathit{M}}$	1.9979	1.9999	1.9999	2.0000
The Present	$E^{\mathit{N},\mathit{M}}$	5.4397e-04	1.1935e-04	2.8367e-05	7.0915e-06	1.7729e-06
	$R^{\mathit{N},\mathit{M}}$	2.1884	2.07291	2.0001	2.0000

6. Conclusion

A finite element method is developed to solve a large time delay, two-parameter singularly perturbed one-dimensional parabolic problem. In the decoupled technique, discretization is done in two phases. First, the time derivative is substituted by the Crank-Nicolson finite difference approximation after the temporal variable is uniformly discretized and approximated by a finite element piecewise linear continuous function. Next, a fitted mesh (layer-adapted) finite element method is applied to the semi-discretized problem in order to perform spatial discretization. Nodal approximations are obtained by deriving and assembling the space element equation. Ultimately, the entire rectangular domain’s finite element solution is provided as a polynomial of t and x, which is linear in t and quadratic in x. We theoretically proved that the method is uniformly convergent, independent of perturbation parameters. We also validate numerical solutions to the governing equation by solving two test problems. The maximum element-wise error decreases uniformly independent of perturbation parameters with nearly second-order convergence as the number of meshes doubles, which is consistent with the theoretical analysis. The logarithmic-scaled graphs demonstrate the same convergence rate. When compared to previous research’s numerical approaches, the current FEM provides a higher order of accuracy.

Disclosure statement

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

Additional information

Notes on contributors

Solomon Woldu Worku

Solomon Woldu Worku received his M.Sc. from Bahir Dar University, Ethiopia. Since 2015, he has been a mathematics lecturer at Dire Dawa University, Ethiopia. He is currently a Ph.D. candidate at Arba Mich University, Ethiopia. He has published one article on the topic of singularly perturbed parabolic equations, and his research interests are numerical analysis, singularly perturbed differential equations, and computational mathematics. Department of Mathematics, Arba Minch University, Arba Minch, Ethiopia.

Gemechis File Duressa

Gemechis File Duressa received his M.Sc. from Addis Ababa University, Ethiopia. Ph.D. from the National Institute of Technology in Warangal, India. He has been the dean of the College of Natural Science at Jimma University for six years. He has been a full professor of mathematics and director of the School of Graduate Studies at Jimma University since 2022. He has published over 120 articles in the fields of Education, numerical methods, numerical analysis, and, more specifically, singularly perturbed differential equations. Currently, he is the administration and development vice president of Jimma University. He also advises thirteen Ph.D. scholars towards completion. Department of Mathematics, College of Natural Sciences, Jimma University, Jimma, Ethiopia.