Streamline Diffusion Finite Element Method for Singularly Perturbed 1D-Parabolic Convection Diffusion Differential Equations with Line Discontinuous Source

Abstract

This article presents a study on singularly perturbed 1D parabolic Dirichlet’s type differential equations with discontinuous source terms on an interior line. The time derivative is discretized using the Euler backward method, followed by the application of the streamline–diffusion finite element method (SDFEM) to solve locally one-dimensional stationary problems on a Shishkin mesh. Our proposed method is shown to achieve first-order convergence in time and second-order convergence in space. Our proposed method offers several advantages over existing techniques, including more accurate approximations of the solution on the boundary layer region, better efficiency, and robustness in dealing with discontinuous line source terms. The numerical examples presented in this paper demonstrate the effectiveness and efficiency of our method, which has practical applications in various fields, such as engineering and applied mathematics. Overall, our proposed method provides an effective and efficient solution to the challenging problem of solving singularly perturbed parabolic differential equations with discontinuous line source terms, making it a valuable tool for researchers and practitioners in various domains.

1. Introduction

In the literature, there are several articles available that deal with the numerical solution of singularly perturbed 1D parabolic differential equations with sufficiently smooth data functions, see [1,2,3,4,5]. Such problems, but with non-smooth data functions, can be seen in [6,7,8]. In [9], Clavero considered a numerical scheme with two small parameters in both the convection and diffusion terms. In [10], Gracia and O’Riordan considered a singularly perturbed reaction–diffusion parabolic problem with an initial condition that was not smooth. In [11], Clavero and Jorge considered 1D singularly perturbed parabolic convection diffusion systems and used a splitting uniformly convergent method. In [12], Yao Cheng, Yanjie Mei and H G Roos considered the local discontinuous Glerkin method for time dependent singularly perturbed convection diffusion problems on layer adapted meshes.

While there have been several studies on solving parabolic differential equations with various boundary conditions, such as Dirichlet and Neumann conditions, the problem of nonlinear parabolic stochastic differential equations with nonlinear Robin conditions remains an active area of research. For instance, the recent article, [13], discusses Well-Posedness for Nonlinear Parabolic Stochastic Differential Equations with Nonlinear Robin Conditions which provides a rigorous analysis of the mathematical properties of this problem. In this study, we focus on the related, but distinct, problem of numerically solving 1D singularly perturbed parabolic differential equations with discontinuous source terms on an interior line.

This type of problem has a regular boundary layer at $x=1$ (boundary point) as the parameter $\epsilon$ tends to zero [14]. As a layer near the boundary exists, conventional discretization techniques, such as Finite Difference Methods (FDMs) or Finite Element Methods (FEMs), cannot yield an accurate solution, unless the mesh is highly refined [14]. Therefore, it is essential that any method proposed must employ a layer-adapted mesh to achieve uniform accuracy. Mukherjee and Natesan [15] developed a hybrid finite difference approach that converges uniformly for singularly perturbed 1D parabolic initial-boundary value problems (IBVPs) on the piecewise uniform Shishkin mesh. Similarly, Das et al. [16] proposed a numerical technique on the Bakhvalov-Shishkin mesh to solve 2D delay parabolic IBVPs. Hughes and Brooks [17] introduced the Streamline–diffusion finite element method (SDFEM), which is widely recognized as an effective technique for obtaining the numerical solutions of convection-dominated flow problems. Later, Roos and Zarin [18] applied SDFEM on the Shishkin mesh to solve a singularly perturbed two-point boundary value problem with a non-smooth source function.

The main focus of this paper was to investigate the numerical treatment of 1D singularly perturbed parabolic Dirichlet’s differential equations with discontinuous source terms on an interior line. This problem contains an interior layer at $x=z$ due to the presence of line discontinuities. The authors propose a method that first uses the backward Euler method to discretize the time derivative, followed by applying SDFEM on the Shishkin mesh to solve the locally one-dimensional stationary problem. In [19], Ghiocel Groza and Nicolae Pop considered a numerical scheme for the locally one-dimensional stationary boundary value problem. Various numerical examples were used to validate the suggested method, both theoretically and numerically, and it demonstrated uniform convergence in both space and time.

The paper is organized as follows: Section 2 presents the statement of the problem, the temporal discretization, derivative estimates and stability findings of locally 1D problems. In Section 3 and Section 4, the weak formulation and the numerical scheme for solving our problem are described. In Section 5, the error estimate for the SDFEM method is provided, while Section 6 offers numerical validation through various test examples. Finally, Section 7 provides some concluding remarks.

2. Continuous Problem and Stability Analysis

2.1. Statement of Continuous Problem

Inspired by the work of [20], the following singularly perturbed 1D parabolic differential equation is investigated in this paper.

Find a function u such that

$$\begin{array}{cc}\hfill \mathfrak{M}u:=& u_t-\epsilon u_{xx}+a\left(x\right)u_x+b\left(x\right)u=\delta (x-z)g\left(t\right)+f(x,t),(x,t)\in {\Omega }^{*}\times (0,T],\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & u(x,0)=u_0\left(x\right),x\in \Omega ,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & u(0,t)=0=u(1,t),t\in [0,T],\hfill \end{array}$$

(3)

where $0<\epsilon \ll 1$ is a very small positive parameter, ${\Omega }^{*}={\Omega }_{-}\cup {\Omega }_{+},{\Omega }_{-}=(0,z),$${\Omega }_{+}=(z,1),\overline{\Omega }=[0,1]$, the functions $a\left(x\right)\ge \alpha >0,b\left(x\right)>\beta >0$, $g\left(t\right)$, $f(x,t)$ and $u_0\left(x\right)$ are sufficiently differentiable and bounded in their respective domains, $\delta (x-z)$ is the delta function, and T is some fixed positive time.

In [18], Roos and Zarin used the SDFEM method to solve a two-point boundary value problem with a point source function, which exhibits singular perturbation and requires a layer adaptive mesh like Shishkin mesh. In their model, there was only one interior layer at a single point. However, in our problem, which involves time evolution, we have an interior layer along a line source. Thus, we need to accurately capture the numerical solution along this line source.

In the paper, the positive constant C is used to denote a generic constant that is not dependent on the perturbation parameter $\epsilon$, or the discretization parameters, such as N or M. For practical purposes, the accepted convention is to assume $\epsilon \le C{N}^{-1}$ for the convection coefficient problem. Additionally, the authors use the supremum norm ${\left|\psi \right|}_D={sup}_{\mathbf{x}\in D}\left|\psi \left(\mathbf{x}\right)\right|$ to measure the error and derivative bounds.

2.2. Time Domain Discretization and Locally 1D Problems

We introduce equidistant meshes in the time domain $[0,T]$ with time step $\Delta t$, such that ${\Omega }_t^M={\{t_i=i\Delta t\}}_{i=0}^M$, $\Delta t=\frac{T}{M}$, where M represents the number of mesh elements in the time direction. Now, discretizing the time derivative by means of the Euler fractional method on uniform mesh, we obtain the following ordinary differential equation for every time step $t_n$ in the set ${\Omega }_t^M$, where n ranges from 1 to M:

$$\begin{array}{cc}\hfill (I+\Delta tL){\tilde{u}}^n=& {\tilde{u}}^{n-1}+\Delta t[f(x,t_n)+\delta (x-z)g\left(t_n\right)],x\in {\Omega }^{*}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & {\tilde{u}}^n\left(0\right)=0={\tilde{u}}^n\left(1\right),\hfill \end{array}$$

(5)

where the differential operator $L=-\epsilon \frac{d^2}{d{x}^2}+a\left(x\right)\frac{d}{dx}+b\left(x\right)$ and ${\tilde{u}}^n$ is the solution of (4) and (5). Note that ${\tilde{u}}^0\left(x\right)=0=u_0\left(x\right),x\in \overline{\Omega }.$

To achieve the stability of the scheme given by (4) and (5), one can easily prove that the operator $\mathcal{Q}=(I+\Delta tL)$ satisfies the maximum principle:

$$\parallel {\mathcal{Q}}^{-1}{\parallel }_{\infty }\le \frac{1}{1+\beta \Delta t}.$$

(6)

We can rewrite the semi-discretized problem (4) and (5) as:

$$\left\{\begin{array}{c}{L}^{*}{\tilde{u}}^n:=-\epsilon {\displaystyle \frac{d^2{\tilde{u}}^n\left(x\right)}{d^{2}x}}+a\left(x\right){\displaystyle \frac{d{\tilde{u}}^n\left(x\right)}{dx}}+c\left(x\right){\tilde{u}}^n\left(x\right)=g^n\left(x\right),x\in {\Omega }^{*},\hfill \\ {\tilde{u}}^n\left(0\right)=0={\tilde{u}}^n\left(1\right)\hfill \end{array}\right.$$

(7)

where $c\left(x\right)=b\left(x\right)+\frac{1}{\Delta t}$ and $g^n\left(x\right)=\frac{1}{\Delta t}{\tilde{u}}^{n-1}+\{f(x,t_n)+\delta (x-z)g\left(t_n\right)\}.$ The above scheme (7) is an ordinary differential equation in space variable x for each time step $t_n$.

2.3. Maximum Principle and Derivative Estimates

Lemma 1

(Maximum Principle). Suppose that there exists a function ξ, belonging to the set $C^0\left(\overline{\Omega }\right)\cap C^2\left({\Omega }^{*}\right)$, that satisfies the following conditions: $\xi \left(x\right)\ge 0$ for $x=0,1,$ $L^{*}\xi \left(x\right)\ge 0$ for all $x\in {\Omega }^{*}$ and ${\xi }^{\prime }\left(z^{-}\right)-{\xi }^{\prime }\left(z^{+}\right)\ge 0$ Then, it can be concluded that $\xi \left(x\right)\ge 0$ for all $x\in \overline{\Omega }$.

Proof. 

The method of proof for the lemma is comparable to the one used in (Lemma 2 in [21]). □

Lemma 2.

Suppose ${\tilde{u}}^n$ is the solution of problem (7) and is decomposed as ${\tilde{u}}^n=r+s$. Then, the derivatives of the regular components satisfy

$$\begin{array}{cc}\hfill & {∥{\displaystyle \frac{d^{l}r\left(x\right)}{d{x}^l}}∥}_{{\Omega }^{*}}\le C(1+{\epsilon }^{-l+2}),l=0,1,2,3,\hfill \end{array}$$

(8)

and, the derivatives of the singular components satisfy

$$\begin{array}{cc}\hfill & \left|{\displaystyle \frac{d^{l}s\left(x\right)}{d{x}^l}}\right|\le C{\epsilon }^{-l}\left\{\begin{array}{c}exp\left({\displaystyle \frac{\alpha (x-z)}{\epsilon }}\right),x\in {\Omega }_{-},l=0,1,2,3,\hfill \\ exp\left({\displaystyle \frac{\alpha (x-1)}{\epsilon }}\right),x\in {\Omega }_{+}.\hfill \end{array}\right.\hfill \end{array}$$

(9)

Proof. 

A similar strategy to the one used in (Theorem 9.1 in [22]) is employed for the proof. □

2.4. Truncation Error

Lemma 3.

Assume that, $|\frac{{\partial }^{i}u}{\partial t^i}|\le C,i=0,1,2,3,$ then $\parallel e_n{\parallel }_{\infty }\le C\Delta t^2$ and $\parallel E_n{\parallel }_{\infty }\le C\Delta t$, where $e_n=u(x,t_n)-{\overline{u}}^n$ and $E_n={\sum }_{j=1}^n{e}_j$.

Proof. 

As ${\tilde{u}}^n$ is the solution of (4) and (5), we have

$$\begin{array}{c}\hfill (I+\Delta tL){\tilde{u}}^n-\Delta t\{f(x,t_n)+\delta (x-z)g\left(t_n\right)\}=u^{n-1}.\end{array}$$

(10)

By means of Taylor expansion, we have

$$\begin{array}{c}\hfill u(x,t_{n-1})=u(x,t_n)-\Delta t{\displaystyle \frac{\partial u}{\partial t}}(x,t_n)+O(\Delta t^2).\end{array}$$

(11)

From Equation (1) and if the solution of (1) is smooth enough, we have

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial u}{\partial t}}(x,t_n)=\{f(x,t_n)+\delta (x-z)g\left(t_n\right)\}-L{u}^n.\end{array}$$

(12)

From (10)–(12), we have

$$\begin{array}{c}\hfill (I+\Delta tL)e_n=O(\Delta t^2),e_n\left(0\right)=e_n\left(1\right)=0.\end{array}$$

From this, we have $\parallel e_n{\parallel \le C(\Delta t)}^2$ by using the stability result in (6).

The stability result and consistency property of (4) and (5) together implies $\parallel E_n\parallel \le C\Delta t.$ □

3. Weak Formulation

The standard weak formulation of problem (7) for a fixed n is given as follows. Find ${\tilde{u}}^n\in V=H_0^1\left(\overline{\Omega }\right)$, such that

$$B({\tilde{u}}^n,v)=g^n\left(v\right),\forall v\in V,$$

where

$$B(y,v)=\epsilon (y^{\prime },v^{\prime })+(a{y}^{\prime }+cy,v)$$

$$g^n\left(v\right)=(g^n,v).$$

Here, (. , .) represents the inner product in $L_2\left({\Omega }^{*}\right).$

Consider the mesh in space as ${\overline{\Omega }}_x^N=\{x_0,x_1,\dots ,x_N\}$, where N is some positive integer. We define a mesh that includes the point z as one of its nodes. Let ${\varphi }_i$ be the basis for the finite-dimensional subspace $V_h$ of piece-wise linear polynomials. The basis is given by

$$\begin{array}{c}\hfill {\varphi }_i\left(x\right)=\left\{\begin{array}{cc}\frac{x-x_{i-1}}{h_i}\hfill & x_{i-1}\le x\le x_i,\hfill \\ \frac{x_{i+1}-x}{h_{i+1}}\hfill & x_i\le x\le x_{i+1},\hfill \\ 0\hfill & otherwise,\hfill \end{array}\right.\end{array}$$

where $h_i=x_i-x_{i-1}$ and $i\in \{1,2,\dots ,N\}$. We use this basis to ensure that z is one of the mesh points.

The standard Galerkin method is given as follows: find ${\tilde{u}}_h^n\in V_h\subset V=H_0^1\left(\overline{\Omega }\right)$ such that

$$B_h({\tilde{u}}_h^n,v_h)=g_h^n\left(v_h\right),\forall v_h\in V_h.$$

4. Streamline Diffusion Finite Element Formulation

The streamline diffusion weak formulation for the problem (7) is to find ${\tilde{u}}_h^n\in V_h\subset V=H_0^1\left(\overline{\Omega }\right)$, such that

$$\begin{array}{c}\hfill B_h({\tilde{u}}_h^n,v_h)=g_h^n\left(v_h\right),\forall v\in V_h,\end{array}$$

(13)

where

$$\begin{array}{c}\hfill B_h(y,w)=\epsilon (y^{\prime },w^{\prime })+(a{y}^{\prime }+cy,w)+\sum _{i=1}^N{\int }_{x_{i-1}}^{x_i}{\delta }_i(-\epsilon y^{\prime \prime }+a{y}^{\prime }+cy)a{w}^{\prime }dx,\end{array}$$

$$\begin{array}{c}\hfill g_h^n\left(w\right)=(g^n,w)+\sum _{i=1}^N{\int }_{x_{i-1}}^{x_i}{\delta }_i{g}^n\left(x\right)a\left(x\right)w^{\prime }dx.\end{array}$$

The parameter, known as the Streamline–diffusion (SD) parameter, denoted as ${\delta }_i$, is decided later.

The SDFEM’s relevant difference scheme is presented below.

$$\begin{array}{c}\hfill L_{SD}^N:=-\epsilon [D^{+}{u}_i-D^{-}{u}_i]+{\alpha }_i{D}^{-}{u}_i+{\beta }_i{D}^{+}{u}_i+{\gamma }_i{u}_i=g_h\left({\varphi }_i\right),\end{array}$$

(14)

where the symbols $D^{+}$ and $D^{-}$ are given by the following:

$$D^{+}{u}_i=\frac{u_{i+1}-u_i}{h_{i+1}},D^{-}{u}_i=\frac{u_i-u_{i-1}}{h_i},$$

and

$$\begin{array}{cc}\hfill {\alpha }_i=& {\int }_{x_{i-1}}^{x_i}a\left(x\right){\varphi }_{i}dx-h_i{\int }_{x_{i-1}}^{x_i}c\left(x\right){\varphi }_{i-1}{\varphi }_{i}dx+\frac{{\delta }_i}{h_i}{\int }_{x_{i-1}}^{x_i}a{\left(x\right)}^{2}dx\hfill \\ \hfill & -{\delta }_i{\int }_{x_{i-1}}^{x_i}a\left(x\right)c\left(x\right){\varphi }_{i-1}dx,\hfill \\ \hfill {\beta }_i=& {\int }_{x_i}^{x_{i+1}}a\left(x\right){\varphi }_{i}dx+h_{i+1}{\int }_{x_i}^{x_{i+1}}c\left(x\right){\varphi }_{i+1}{\varphi }_{i}dx-\frac{{\delta }_{i+1}}{h_{i+i}}{\int }_{x_i}^{x_{i+1}}a{\left(x\right)}^{2}dx\hfill \\ \hfill & -{\delta }_{i+1}{\int }_{x_i}^{x_{i+1}}a\left(x\right)c\left(x\right){\varphi }_{i+1}dx,\hfill \\ \hfill {\gamma }_i=& {\int }_{x_{i-1}}^{x_{i+1}}c\left(x\right){\varphi }_{i}dx+\frac{{\delta }_i}{h_i}{\int }_{x_{i-1}}^{x_i}a\left(x\right)c\left(x\right)dx-\frac{{\delta }_{i+1}}{h_{i+1}}{\int }_{x_i}^{x_{i+1}}a\left(x\right)c\left(x\right)dx.\hfill \end{array}$$

The standard Galerkin method is effective for step lengths that are sufficiently small. If this condition is not met, the method can be stabilized by utilizing the characteristics of an M-matrix. To satisfy this condition, the streamline–diffusion parameter ${\delta }_i$ can be determined such that the matrix resulting from the associated difference scheme (14) transforms into an M-matrix, as outlined below:

$${[-{\alpha }_i]}_{{\delta }_i=0}\le \epsilon ,$$

This can be written as ${\overline{C}}_i{h}_i\le {\overline{A}}_i+\epsilon$, where ${\overline{A}}_i={\int }_{x_{i-1}}^{x_i}a\left(x\right){\varphi }_{i}dx$ and ${\overline{C}}_i={\int }_{x_{i-1}}^{x_i}c\left(x\right){\varphi }_{i-1}{\varphi }_{i}dx.$ If this condition is not met, ${\alpha }_i$ is set to zero for $i\in \{1,2,\dots ,N\}$ to achieve an M-matrix, which results in the following:

$${\delta }_i=\frac{(h_i{\overline{C}}_i-{\overline{A}}_i)h_i}{{\overline{A^2}}_i-h_i{\overline{AC}}_i}$$

where, ${\overline{A^2}}_i={\int }_{x_{i-1}}^{x_i}{\left(a\left(x\right)\right)}^{2}dx$ and ${\overline{AC}}_i={\int }_{x_{i-1}}^{x_i}a\left(x\right)c\left(x\right){\varphi }_{i-1}dx.$ If we summarize the above conditions, then we obtain

$$\begin{array}{c}\hfill {\delta }_i=\left\{\begin{array}{c}0,if{\overline{C}}_i{h}_i\le {\overline{A}}_i+\epsilon \hfill \\ \frac{(h_i{\overline{C}}_i-{\overline{A}}_i)h_i}{{\overline{A^2}}_i-h_i{\overline{AC}}_i},otherwise.\hfill \end{array}\right.\end{array}$$

(15)

As a result, it is evident that ${\delta }_i=O\left(h_i\right)$. By choosing $v_h$ as ${\varphi }_i$ as well as ${\tilde{u}}_h=({\tilde{u}}_h^n\left(x_1\right),$${\tilde{u}}_h^n\left(x_2\right),$$\dots ,{\tilde{u}}_h^n\left(x_N\right))\in R^N$, we obtain the linear system of algebraic equations associated with the SDFEM scheme

$$K{\tilde{u}}_h=g_h,$$

where K represents a tridiagonal matrix $K=\left(k_{ij}\right)$ and $g_h=g_h^n\left({\varphi }_i\right)$. Upon computing the coefficients,

$$\begin{array}{cc}\hfill k_{i,i-1}=& -\epsilon h_i^{-1}-h_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right){\varphi }_{i}dx+{\int }_{x_{i-1}}^{x_i}c\left(x\right){\varphi }_{i-1}{\varphi }_{i}dx-{\delta }_i{h}_i^{-2}{\int }_{x_{i-1}}^{x_i}a{\left(x\right)}^{2}dx\hfill \\ \hfill & +{\delta }_i{h}_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right)c\left(x\right){\varphi }_{i-1}dx,\hfill \\ \hfill k_{i,i}=& \epsilon h_i^{-1}+\epsilon h_{i+1}^{-1}+h_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right){\varphi }_{i}dx-h_{i+1}^{-1}{\int }_{x_i}^{x_{i+1}}c\left(x\right){\varphi }_{i}dx+{\int }_{x_{i-1}}^{x_i}c\left(x\right){\varphi }_i^{2}dx\hfill \\ \hfill & +{\int }_{x_i}^{x_{i+1}}c\left(x\right){\varphi }_i^{2}dx+{\delta }_i{h}_i^{-2}{\int }_{x_{i-1}}^{x_i}a{\left(x\right)}^{2}dx\hfill \\ \hfill & +{\delta }_{i+1}{h}_{i+1}^{-2}{\int }_{x_i}^{x_{i+1}}a{\left(x\right)}^{2}dx+{\delta }_i{h}_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right)c\left(x\right){\varphi }_{i}dx-{\delta }_{i+1}{h}_{i+1}^{-1}{\int }_{x_i}^{x_{i+1}}a\left(x\right)c\left(x\right){\varphi }_{i}dx\hfill \\ \hfill k_{i,i+1}=& -\epsilon h_{i+1}^{-1}+h_{i+1}^{-1}{\int }_{x_i}^{x_{i+1}}a\left(x\right){\varphi }_{i}dx+{\int }_{x_i}^{x_{i+1}}c\left(x\right){\varphi }_{i+1}{\varphi }_{i}dx\hfill \\ \hfill & -{\delta }_{i+1}{h}_{i+1}^{-2}{\int }_{x_i}^{x_{i+1}}a{\left(x\right)}^{2}dx-{\delta }_{i+1}{h}_{i+1}^{-1}{\int }_{x_i}^{x_{i+1}}a\left(x\right)c\left(x\right){\varphi }_{i+1}dx.\hfill \end{array}$$

Remark 1.

We defined the SD parameter using the procedure given in [23]. Realistically, the conventional Galerkin method satisfies the criterion that defines an M-matrix for almost all 1D problems. In such a situation, the Galerkin method and the SDFEM both yield almost identical outcomes.

Remark 2.

Generally, the conventional Galerkin method is analogous to the central finite difference approach. However, in regions where convection is dominant, applying a central finite difference approximation to the convective term may lead to oscillations [14]. In such instances, the SDFEM is a superior alternative.

4.1. Discrete Green’s Function and Stability

To establish the $(l_{\infty },w_{-1,\infty })$ stability of the SDFEM, we introduce the discrete Green’s function. The i-th discrete Green’s function, denoted as ${\lambda }^i\in V_h$, is determined by solving the following problem:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{B}_h({\varphi }_j,{\lambda }^i)={\delta }_{ij},\hfill \\ {\lambda }^i\left(0\right)={\lambda }^i\left(1\right)=0,\hfill \end{array}\right.\end{array}$$

(16)

where ${\delta }_{ij}$ denotes the Kronecker delta. Additional details on these and related concepts can be found in references such as [24,25,26].

Lemma 4

([25]). If ${\lambda }^i={\sum }_{j=1}^N{\lambda }_j^i{\varphi }_j,$, then the discrete Green’s function exhibits the following conditions:

(1) 

${\lambda }_j^i\ge 0,i,j\in \{1,2,\dots ,N\},$

(2) 

$0\le {\lambda }_1^i<\dots <{\lambda }_i^i>{\lambda }_{i+1}^i\dots >{\lambda }_N^i.$

Lemma 5.

If we denote $c_0={min}_i\left\{c_i\right\},i=1,2,\dots ,N$ where $c_i=k_{i+1,i}-k_{i,i+1}$ and ${\lambda }^i={\sum }_{j=1}^N{\lambda }_j^i{\varphi }_j,$ then ${\lambda }_i^i\le c_0^{-1}.$

Proof. 

From (16), for any $v_h\in V_h,$ we have

$$B_h(v_h,{\lambda }^i)=v_h\left(x_i\right).$$

Considering $v_h={\sum }_{j=1}^i{\varphi }_j,$ one can derive

$$\sum _{j=1}^N\sum _{k=1}^j{k}_{j,k}{\lambda }_j^i=1.$$

Expanding, we obtain

$$\begin{array}{c}(k_{1,1}+k_{1,2}){\lambda }_1^i+\sum _{j=2}^i(k_{j,j-1}+k_{j,j}+k_{j,j+1}){\lambda }_j^i+\\ k_{i+1,i}({\lambda }_{i+1}^i-{\lambda }_i^i)+{\lambda }_i^i(k_{i+1,i}-k_{i,i+1})=1.\end{array}$$

K is an M-matrix implying $(k_{i,i-1}+k_{i,i}+k_{i,i+1})\ge 0,$ for $i\in \{2,\dots ,N-1\}$. Now, by means of Lemma 4 and $k_{1,1}>k_{1,2}$ (property), we have

$${\lambda }_i^i(k_{i+1,i}-k_{i,i+1})\le 1.$$

This implies that

$${\lambda }_i^i\le \frac{1}{(k_{i+1,i}-k_{i,i+1})}=c_i^{-1}.$$

So, we conclude that ${\lambda }_i^i\le c_0^{-1}.$ □

Lemma 6

([25]). The SDFEM, equipped with the streamline–diffusion parameter specified in Equation (15), is uniformly stable in the $(l_{\infty },w_{-1,\infty })$ norm, as demonstrated by the following result:

$$\parallel v_h{\parallel }_{\infty }\le \frac{2}{c_0}\underset{j}{max}\left|\sum _{k=j}^N{\left(K{v}_h\right)}_k\right|,\forall v_h\in V_h,j=1,\dots ,N.$$

4.2. Shishkin Type Mesh

We used a general type mesh introduced in [27] with adapted layers at the points $x=z$ and $x=1$. Let N be a positive even integer greater than 4 and ${\sigma }_1=min\{\frac{z}{2},\frac{{\tau }_0}{\alpha }\epsilon lnN\}$, ${\sigma }_2=min\{\frac{1-z}{2},\frac{{\tau }_0}{\alpha }\epsilon lnN\}$, ${\tau }_0\ge 2$. Here, we considered ${\sigma }_1={\sigma }_2=\frac{{\tau }_0}{\alpha }\epsilon lnN$. Let ${\Omega }_s=(0,z-{\sigma }_1)\cup (z,1-{\sigma }_2)$ and ${\Omega }_0=(z-{\sigma }_1,z)\cup (1-{\sigma }_2,1)$. On ${\Omega }_s$ and ${\Omega }_o$, the mesh is equidistant and graded respectively. The transition points are chosen to be

$$x_{\frac{N}{4}}=z-{\sigma }_1,x_{\frac{N}{2}}=z,x_{\frac{3N}{4}}=1-{\sigma }_2.$$

We chose two mesh generating functions ${\phi }_1$ and ${\phi }_2$ for the particular layers:

$${\phi }_1(\frac{1}{4})=lnN,{\phi }_1(\frac{1}{2})=0,$$

$${\phi }_2(\frac{3}{4})=lnN,{\phi }_2(1)=0.$$

The nodes of the Shishkin mesh (S-mesh) are as follows:

$$\begin{array}{c}\hfill x_i=\left\{\begin{array}{cc}\frac{4i}{N}(z-{\sigma }_1)\hfill & i=0,\dots ,\frac{N}{4},\hfill \\ z-\frac{{\tau }_0}{\alpha }\epsilon {\phi }_1\left(t_i\right)\hfill & i=\frac{N}{4}+1,\dots ,\frac{N}{2},\hfill \\ z+\frac{4(1-z-{\sigma }_2)(i-\frac{N}{2})}{N}\hfill & i=\frac{N}{2}+1,\dots ,\frac{3N}{4},\hfill \\ -\frac{{\tau }_0}{\alpha }\epsilon {\phi }_2\left(t_i\right)+1\hfill & i=\frac{3N}{4}+1,\dots ,N,\hfill \end{array}\right.\end{array}$$

(17)

where $t_i=\frac{i}{N},{\phi }_i=-ln{\psi }_i$ for $i=1,2,{\psi }_1\left(t\right)=exp(-(2-4t)lnN)$ and ${\psi }_2\left(t\right)=exp(-(4-4t)lnN)$ on Shishkin mesh (S-mesh).

${max}_i\left|{\psi }_i^{\prime }\right|$, i = 1, 2 plays a part in error analysis as we see in the upcoming section. we have $max|{\psi }^{\prime }|\le ClnN$ in the S-mesh. Assuming that the mesh generating function ${\phi }_i$, for $i=1,2,$ obeys the following:

$$\begin{array}{c}\hfill max\left\{\right|{\phi }_i^{\prime }\left|\right\}\le N{C}_0,fori=1,2.\end{array}$$

(18)

We observe that on the coarse part ${\Omega }_s$, the following holds for the S-mesh:

$$\begin{array}{c}\hfill h_{i}N\le C.\end{array}$$

(19)

On the layer part it is true that on the S-mesh,

$$\begin{array}{c}\hfill h_i\le C\epsilon \frac{lnN}{N}\end{array}$$

(20)

and, by the assumption (18), we obtain

$$\begin{array}{c}\hfill \frac{h_i}{\epsilon }\le C{N}^{-1}max\left|{\phi }^{\prime }\right|\le C.\end{array}$$

(21)

5. Error Estimations

The error at each time level $t_n$, where $1\le n\le M-1$, can be expressed as the difference between ${\tilde{u}}^n\left(x_i\right)$ and ${\tilde{u}}_h^n\left(x_i\right)$ for all $x_i\in {\Omega }^{*}$. Another way to represent this error is by utilizing the linear interpolant ${\tilde{u}}_I^n$ of ${\tilde{u}}^n$ as follows:

$$\begin{array}{cc}\hfill e\left(x_i\right)=({\tilde{u}}^n-{\tilde{u}}_h^n)\left(x_i\right)& =({\tilde{u}}^n-{\tilde{u}}_I^n)\left(x_i\right)+({\tilde{u}}_I^n-{\tilde{u}}_h^n)\left(x_i\right)\hfill \\ \hfill & =e^I+e^D.\hfill \end{array}$$

The errors of interpolation and discretization are denoted by $e^I$ and $e^D$, respectively. In this section, we begin by deriving an estimate for the upper bound of the discretization error, based on the error caused by interpolation. Following that, we estimate the upper bound of the interpolation error.

The discretization error is given by

$$\begin{array}{c}\hfill e^D=({\tilde{u}}_I^n-{\tilde{u}}_h^n)\left(x\right)=\sum _{i=0}^N{e}_i{\varphi }_i.\end{array}$$

(22)

By means of the orthogonality property, we have

$$\begin{array}{c}\hfill B_h({\tilde{u}}^n-{\tilde{u}}_h^n,v_h)=0,\forall v_h\in V_h.\end{array}$$

Using the above property, we obtain

$$\begin{array}{cc}\hfill B_h(e^D,{\varphi }_i)& =B_h({\tilde{u}}_I^n-{\tilde{u}}_h^n,{\varphi }_i)\hfill \\ \hfill & =B_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_i),i\in \{1,2,\dots ,N-1\}.\hfill \end{array}$$

Therefore, the equation for error is

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill & B_h(e^D,{\varphi }_i)=B_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_i),i\in \{1,2,\dots ,N-1\},\hfill \\ \hfill & e_0^D=0,e_N^D=0.\hfill \end{array}\right.\end{array}$$

(23)

The following lemma gives the explicit expression for $B_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_i)$.

Lemma 7.

The error Equation (23) involves the bilinear form $B_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_i)$, which can be explicitly expressed as follows:

$$B_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_i)=O_i+(P_i+P_{i+1})+(Q_i-Q_{i+1})+(R_i-R_{i+1})+(S_i-S_{i+1}),$$

where

$$\begin{array}{cc}\hfill O_i=& {\int }_{x_{i-1}}^{x_{i+1}}{a}^{\prime }\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n){\varphi }_{i}dx,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill P_i=& {\int }_{x_{i-1}}^{x_i}c\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n){\varphi }_{i}dx,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill Q_i=& \frac{{\delta }_i}{h_i}{\int }_{x_{i-1}}^{x_i}(a\left(x\right)c\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n)+a{\left(x\right)}^2{({\tilde{u}}_I^n-{\tilde{u}}^n)}^{\prime })dx,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill R_i=& -\epsilon {\delta }_i{h}_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right){({\tilde{u}}_I^n-{\tilde{u}}^n)}^{\prime \prime }dx,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill S_i=& h_i^{-1}{\int }_{x_{i-1}}^{x_i}a\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n)dx.\hfill \end{array}$$

(28)

Proof. 

One can prove the required results by a simple calculation. □

Lemma 8.

Assume that $a\left(x\right)\in W^{2,1},$ then the following estimates are true:

$$\begin{array}{cc}\hfill \left(i\right)& \sum _{i=1}^N|O_i|\le c_1\parallel {\tilde{u}}^n-{\tilde{u}}_I^n{\parallel }_{\infty },\left(ii\right)\sum _{i=1}^N|P_i|\le c_2{\parallel {\tilde{u}}^n-{\tilde{u}}_I^n\parallel }_{\infty },\hfill \\ \hfill \left(iii\right)& \parallel Q_i{\parallel }_{\infty }\le c_3\parallel {\tilde{u}}^n-{\tilde{u}}_I^n{\parallel }_{\infty },\left(iv\right)\parallel R_i{\parallel }_{\infty }\le c_4{\parallel {\tilde{u}}^n-{\tilde{u}}_I^n\parallel }_{\infty },\hfill \\ \hfill \left(v\right)& \parallel S_i{\parallel }_{\infty }\le c_5{\parallel {\tilde{u}}^n-{\tilde{u}}_I^n\parallel }_{\infty }.\hfill \end{array}$$

Proof. 

We prove this Lemma sequentially:

(i).

From (24), we know that

$$O_i={\int }_{x_{i-1}}^{x_{i+1}}{a}^{\prime }\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n){\varphi }_{i}dx.$$

Now, taking all the intervals and summing up, we have

$$\begin{array}{cc}\hfill \sum _{i=1}^N\left|O_i\right|& \le \parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }\sum _{i=1}^N{\int }_{x_{i-1}}^{x_{i+1}}{a}^{\prime }\left(x\right)dx\hfill \\ \hfill & \le \parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }{\left|a\left(x\right)\right|}_{1,1}\hfill \\ \hfill & =c_1{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.\hfill \end{array}$$

(ii).

From Equation (25), we have

$$P_i={\int }_{x_{i-1}}^{x_i}c\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n){\varphi }_{i}dx.$$

If we conduct a summation over all intervals, and derive

$$\begin{array}{cc}\hfill \sum _{i=1}^N\left|P_i\right|& \le \parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }{\parallel a\left(x\right)\parallel }_{\infty }{h}_i\hfill \\ \hfill & =c_2{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.\hfill \end{array}$$

(iii).

From Equation (26), we have

$$Q_i=\frac{{\delta }_i}{h_i}{\int }_{x_{i-1}}^{x_i}(a\left(x\right)c\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n)+a{\left(x\right)}^2{({\tilde{u}}_I^n-{\tilde{u}}^n)}^{\prime })dx.$$

If the inequality ${\overline{C}}_i{h}_i\le {\overline{A}}_i+\epsilon$ holds, then ${\delta }_i=0$ and, as a result, $r_i=0$. On the other hand, if the condition is not satisfied, then we have ${\delta }_i\ne 0$. In this scenario, it can be observed that ${\delta }_i=O\left(h_i\right)$, we obtain

$$\begin{array}{cc}\hfill |Q_i|& \le {C\parallel a\left(x\right)\parallel }_{\infty }\parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }{(\parallel c\left(x\right)\parallel }_{\infty }{h}_i-{2\left|a\left(x\right)\right|}_{1,1})\hfill \\ \hfill & =c_3{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.\hfill \end{array}$$

(iv).

From Equation (27), we have

$$R_i=-\frac{\epsilon {\delta }_i}{h_i}{\int }_{x_{i-1}}^{x_i}a\left(x\right){({\tilde{u}}_I^n-{\tilde{u}}^n)}^{\prime \prime }dx,$$

and, here also we have $r_i=0$ when ${\overline{C}}_i{h}_i\le {\overline{A}}_i+\epsilon$, and if it is not true, we have

$$\begin{array}{cc}\hfill |R_i|& \le {C\left|a\left(x\right)\right|}_{2,1}{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }\hfill \\ \hfill & =c_4{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.\hfill \end{array}$$

(v).

Now, finally from Equation (28), we have

$$S_i=\frac{1}{h_i}{\int }_{x_{i-1}}^{x_i}a\left(x\right)({\tilde{u}}_I^n-{\tilde{u}}^n)dx.$$

Simplifying this, we obtain

$$\begin{array}{cc}\hfill |S_i|& \le {\parallel a\left(x\right)\parallel }_{\infty }{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }\hfill \\ \hfill & =c_5{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.\hfill \end{array}$$

□

Lemma 9.

Consider the problem (7) and let ${\tilde{u}}^n$ be the exact solution to this problem. Furthermore, let ${\tilde{u}}_I^n$ denote the interpolant of ${\tilde{u}}^n$ on a given grid and let ${\tilde{u}}_h^n$ represent the approximate solution at the time level $t_n$. Then, It follows that

$$\parallel {\tilde{u}}^n-{\tilde{u}}_h^n{\parallel }_{\infty }\le C{\parallel {\tilde{u}}_I^n-{\tilde{u}}^n\parallel }_{\infty }.$$

Proof. 

Using Lemmas 6 and 8, we have

$$\begin{array}{cc}\hfill \parallel {\tilde{u}}_I^n-{\tilde{u}}_h^n{\parallel }_{\infty }& \le C\underset{i\in \{1,\dots ,N\}}{max}\left|\sum _{j=i}^N{\left(K({\tilde{u}}_I^n-{\tilde{u}}_h^n)\right)}_j\right|\hfill \\ \hfill & =C\underset{i\in \{1,\dots ,N\}}{max}\left|\sum _{j=i}^N{B}_h({\tilde{u}}_I^n-{\tilde{u}}^n,{\varphi }_j)\right|\hfill \\ \hfill & \le C\left(\left|\sum _{i=1}^N{O}_i\right|+2\left|\sum _{i=1}^N{P}_i\right|+\underset{i\in \{1,\dots ,N\}}{max}(Q_i-Q_N)+\underset{i\in \{1,\dots ,N\}}{max}(R_i-R_N)\right)\hfill \\ & +C\underset{i\in \{1,\dots ,N\}}{max}(S_i-S_N)\hfill \\ \hfill & \le C\left(\sum _{i=1}^N|O_i|+2\sum _{i=1}^N|P_i|+2(\parallel Q_i{\parallel }_{\infty }+\parallel R_i{\parallel }_{\infty }+\parallel S_i{\parallel }_{\infty })\right)\hfill \\ \hfill & \le C\parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }.\hfill \end{array}$$

This implies that the discretization error $e^D$ given in (22) is bounded in terms of the interpolation error $e^I.$ By applying the triangle inequality, we establish the following inequality:

$$\parallel {\tilde{u}}^n-{\tilde{u}}_h^n{\parallel }_{\infty }\le C|{\tilde{u}}_I^n-{\tilde{u}}^n\parallel \infty .$$

Therefore, we have successfully demonstrated the validity of this inequality. □

We now proceed to estimate the upper bound of the interpolation error for the S-mesh.

Theorem 1.

Let ${\tilde{u}}^n$ be the classical solution of the problem (7), ${\tilde{u}}_I^n$ be its interpolant on a Shishkin mesh with the grid points $x_i$ as in (17) and $\epsilon \le \frac{C}{N}.$ Then, the following holds:

$$\begin{array}{c}\hfill \parallel {\tilde{u}}_I^n-{\tilde{u}}^n{\parallel }_{\infty }\le \left\{\begin{array}{c}C{(N^{-1}lnN)}^2,x\in {\Omega }_0,\hfill \\ C{N}^{-2},x\in {\Omega }_s.\hfill \end{array}\right.\end{array}$$

Proof. 

We perform a separate analysis of the error due to interpolation on the domain ${\Omega }_{-}$ as well as ${\Omega }_{+}$.

First, consider x in the domain ${\Omega }_{-}$. Let ${\tilde{u}}_1^n\left(x\right)$ be the solution to the problem:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{L}^{*}{\tilde{u}}_1^n\left(x\right)=g^n\left(x\right),x\in {\Omega }_{-}=(0,z),\hfill \\ {\tilde{u}}^0\left(0\right)=0,{\tilde{u}}^0\left(z\right)=d.\hfill \end{array}\right.\end{array}$$

For now, assume that d is a constant. Then, the solution ${\tilde{u}}_1^n\left(x\right)$ exists and we can decompose it as a sum of two functions $r_1$ and $s_1$, where $r_1$ and $s_1$ satisfy bounds (8) and (9). Then, we can write

$${\tilde{u}}_1^n\left(x\right)-{\tilde{u}}_{1,I}^n\left(x\right)=r_1\left(x\right)-r_{1,I}\left(x\right)+s_1\left(x\right)-s_{1,I}\left(x\right).$$

By means of classical theory, and from (19) and (8), we have the following estimate for the interpolation error $r_1\left(x\right)-r_{1,I}\left(x\right)$ on the regular part:

$$|r_1\left(x\right)-r_{1,I}\left(x\right)|\le C{h}_i^2\underset{x_{i-1}\le x\le x_i}{max}\left|r_1^{\prime \prime }\right|\le C{h}_i^2\le C{N}^{-2}.$$

For the interpolation error $r_1\left(x\right)-r_{1,I}\left(x\right)$ in the layer part, one can derive:

$$\begin{array}{cc}\hfill |r_1\left(x\right)-r_{1,I}\left(x\right)|& \le C{\epsilon }^2{(N^{-1}lnN)}^{2}exp\left(\frac{-2\alpha }{{\tau }_0\epsilon }(x_{i-1}-z)\right)\hfill \\ \hfill & \le C{(N^{-1}lnN)}^2.\hfill \end{array}$$

The above expression is obtained by utilizing the inequality (20), selecting the transition point ${\sigma }_1$, and assuming the condition $\epsilon \le \frac{C}{N}$ holds. Now, for the interpolation error $s_1\left(x\right)-s_{1,I}\left(x\right)$ on the regular part, we have

$$\begin{array}{cc}\hfill |s_1\left(x\right)-s_{1,I}\left(x\right)|\le 2|s_1\left(x\right)|& \le C\underset{x_{i-1}\le x\le x_i}{max}exp\left(-\frac{\alpha }{\epsilon }(z-x)\right)\hfill \\ \hfill & \le Cexp\left(-\frac{\alpha }{\epsilon }{\sigma }_1\right)\hfill \\ \hfill & =C{N}^{-{\tau }_0}.\hfill \end{array}$$

Regarding the interpolation error for the layer component, i.e., $s_1\left(x\right)-s_{1,I}\left(x\right)$, we apply the classical theory to obtain the following:

$$\begin{array}{cc}\hfill |s_1\left(x\right)-s_{1,I}\left(x\right)|& \le C{h}_i^2\underset{x_{i-1}\le x\le x_i}{max}\left|s_1^{\prime \prime }\right|\hfill \\ \hfill & \le C{(N^{-1}lnN)}^2\underset{x_{i-1}\le x\le x_i}{max}exp\left(\frac{-\alpha }{\epsilon }(z-x)\right)\hfill \\ \hfill & \le C{(N^{-1}lnN)}^2.\hfill \end{array}$$

Here, we used the inequality (20), the first inequality of (9) and (21). Similarly, we can prove the case $x\in {\Omega }_{+}$. □

Theorem 2.

Suppose that $u(x,t)$ represents the classical solution of problem (1)–(3), ${\tilde{u}}_h^n\left(x\right)$ represents the numerical solution of the SDFEM scheme (fully discrete) (13) at the $n^{th}$ time step $t_n$. Under the conditions of Lemma 3 and Theorem 1, we have the error estimate as follows:

$$\begin{array}{c}\hfill \parallel {\tilde{u}}_h^n\left(x\right)-u(x,t_n){\parallel }_{\infty }\le \left\{\begin{array}{c}C({(N^{-1}lnN)}^2+\Delta t),x\in {\Omega }_0,\hfill \\ C(N^{-2}+\Delta t),x\in {\Omega }_s.\hfill \end{array}\right.\end{array}$$

Proof. 

We can prove this error estimate by combining the results from Lemma 3 and Theorem 1. □

6. Numerical Validation

In order to demonstrate the accuracy of the theoretical findings, we provide two illustrative examples in this section. We employed the double mesh principle to estimate the errors and their convergence rates for the test problems, as the exact solutions are unknown. The principle involves obtaining the numerical solution $Y^{2N,k/2}(x_i^n,t_i)$ on a grid ${\Omega }_x^{2N}\times {\Omega }_t^{2M}$, where the spatial direction is divided into 2N intervals, while the temporal direction is divided into 2M intervals. The mesh ${\Omega }_x^{2N}\times {\Omega }_t^{2M}$ is obtained by dividing each segment of the previous mesh ${\Omega }_x^N\times {\Omega }_t^M$, where the spatial direction is divided into N intervals, while the temporal direction is divided into M intervals, in two equal parts. Subsequently, we estimate the maximum point-wise error and convergence rate for each $\epsilon$.

$$\begin{array}{c}\hfill E_{\epsilon }^{N,k}=\underset{i,n}{max}\left|Y^{N,k}(x_i^n,t_n)-Y^{2N,k/2}(x_i^n,t_n)\right|,\end{array}$$

$$\begin{array}{c}\hfill D_x^{N,k}=\underset{\epsilon }{max}{E}_{\epsilon }^{N,k},r^{N,k}=lo{g}_2\left(\frac{D_x^{N,k}}{D_x^{2N,k/2}}\right).\end{array}$$

Example 1.

We examined the one-dimensional parabolic partial differential equation with a line source, described by Equations (1)–(3). The given data for this problem were as follows:

$$\begin{array}{cc}\hfill & u(x,0)=0,x\in \Omega =(0,1),\hfill \\ \hfill & u(0,t)=0,u(1,t)=0,t\in [0,1],\hfill \\ \hfill & a\left(x\right)=x+5;b\left(x\right)=1;f(x,t)=x+t;g\left(t\right)=1;z=\frac{1}{2}.\hfill \end{array}$$

In

Table 1. Maximum error and order of convergence for Example 1 with the number M = 50,100.

$\mathit{\epsilon }$	M	N (Number of Grid Points)
↓	↓	${\mathbf{2}}^{\mathbf{6}}$	${\mathbf{2}}^{\mathbf{7}}$	${\mathbf{2}}^{\mathbf{8}}$	${\mathbf{2}}^{\mathbf{9}}$	${\mathbf{2}}^{\mathbf{10}}$	${\mathbf{2}}^{\mathbf{11}}$	${\mathbf{2}}^{\mathbf{12}}$
$2^{-6}$	50	1.0843 $\times {10}^{-1}$	3.1740 $\times {10}^{-2}$	1.2410 $\times {10}^{-2}$	4.0043 $\times {10}^{-3}$	1.2979 $\times {10}^{-3}$	3.6024 $\times {10}^{-4}$	1.0865 $\times {10}^{-4}$
	100	1.0286 $\times {10}^{-1}$	2.6681 $\times {10}^{-2}$	1.0315 $\times {10}^{-2}$	3.7978 $\times {10}^{-3}$	1.5062 $\times {10}^{-3}$	6.2557 $\times {10}^{-4}$	2.6478 $\times {10}^{-4}$
$2^{-7}$	50	1.1120 $\times {10}^{-1}$	2.8842 $\times {10}^{-2}$	1.1237 $\times {10}^{-2}$	3.5721 $\times {10}^{-3}$	1.1214 $\times {10}^{-3}$	3.9296 $\times {10}^{-4}$	1.6141 $\times {10}^{-4}$
	100	1.0748 $\times {10}^{-1}$	2.8066 $\times {10}^{-2}$	1.1374 $\times {10}^{-2}$	4.6676 $\times {10}^{-3}$	1.8575 $\times {10}^{-3}$	8.9909 $\times {10}^{-4}$	4.2970 $\times {10}^{-4}$
$2^{-8}$	50	1.1289 $\times {10}^{-1}$	2.7113 $\times {10}^{-2}$	1.0495 $\times {10}^{-2}$	3.2872 $\times {10}^{-3}$	1.2105 $\times {10}^{-3}$	4.5429 $\times {10}^{-4}$	1.9634 $\times {10}^{-4}$
	100	1.1002 $\times {10}^{-1}$	2.9689 $\times {10}^{-2}$	1.2033 $\times {10}^{-2}$	5.0398 $\times {10}^{-3}$	2.4171 $\times {10}^{-3}$	1.1748 $\times {10}^{-3}$	5.8837 $\times {10}^{-4}$
$2^{-9}$	50	1.1361 $\times {10}^{-1}$	2.6175 $\times {10}^{-2}$	1.0085 $\times {10}^{-2}$	3.2652 $\times {10}^{-3}$	1.2404 $\times {10}^{-3}$	5.0069 $\times {10}^{-4}$	2.1944 $\times {10}^{-4}$
	100	1.1175 $\times {10}^{-1}$	3.0602 $\times {10}^{-2}$	1.2410 $\times {10}^{-2}$	5.3361 $\times {10}^{-3}$	2.5656 $\times {10}^{-3}$	1.3340 $\times {10}^{-3}$	6.5985 $\times {10}^{-4}$
$2^{-10}$	50	1.1395 $\times {10}^{-1}$	2.5686 $\times {10}^{-2}$	9.8700 $\times {10}^{-3}$	3.3070 $\times {10}^{-3}$	1.2620 $\times {10}^{-3}$	5.1495 $\times {10}^{-4}$	2.3576 $\times {10}^{-4}$
	100	1.1263 $\times {10}^{-1}$	3.1082 $\times {10}^{-2}$	1.2609 $\times {10}^{-2}$	5.5213 $\times {10}^{-3}$	2.6932 $\times {10}^{-3}$	1.3756 $\times {10}^{-3}$	7.0938 $\times {10}^{-4}$
$2^{-11}$	50	1.1411 $\times {10}^{-1}$	2.5438 $\times {10}^{-2}$	9.7602 $\times {10}^{-3}$	3.3285 $\times {10}^{-3}$	1.2731 $\times {10}^{-3}$	5.2436 $\times {10}^{-4}$	2.4101 $\times {10}^{-4}$
	100	1.1307 $\times {10}^{-1}$	3.1327 $\times {10}^{-2}$	1.2711 $\times {10}^{-2}$	5.6167 $\times {10}^{-3}$	2.7588 $\times {10}^{-3}$	1.4141 $\times {10}^{-3}$	7.2555 $\times {10}^{-4}$
$2^{-12}$	50	1.1418 $\times {10}^{-1}$	2.5312 $\times {10}^{-2}$	9.7046 $\times {10}^{-3}$	3.3393 $\times {10}^{-3}$	1.2787 $\times {10}^{-3}$	5.2911 $\times {10}^{-4}$	2.4429 $\times {10}^{-4}$
	100	1.1329 $\times {10}^{-1}$	3.1451 $\times {10}^{-2}$	1.2763 $\times {10}^{-2}$	5.6650 $\times {10}^{-3}$	2.7921 $\times {10}^{-3}$	1.4321 $\times {10}^{-3}$	7.3508 $\times {10}^{-4}$
$2^{-13}$	50	1.1422 $\times {10}^{-1}$	2.5249 $\times {10}^{-2}$	9.6767 $\times {10}^{-3}$	3.3448 $\times {10}^{-3}$	1.2815 $\times {10}^{-3}$	5.3150 $\times {10}^{-4}$	2.4594 $\times {10}^{-4}$
	100	1.1340 $\times {10}^{-1}$	3.1514 $\times {10}^{-2}$	1.2789 $\times {10}^{-2}$	5.6894 $\times {10}^{-3}$	2.8088 $\times {10}^{-3}$	1.4416 $\times {10}^{-3}$	7.3966 $\times {10}^{-4}$
$2^{-14}$	50	1.1424 $\times {10}^{-1}$	2.5217 $\times {10}^{-2}$	9.6627 $\times {10}^{-3}$	3.3475 $\times {10}^{-3}$	1.2829 $\times {10}^{-3}$	5.3270 $\times {10}^{-4}$	2.4679 $\times {10}^{-4}$
	100	1.1345 $\times {10}^{-1}$	3.1545 $\times {10}^{-2}$	1.2803 $\times {10}^{-2}$	5.7016 $\times {10}^{-3}$	2.8172 $\times {10}^{-3}$	1.4469 $\times {10}^{-3}$	7.4190 $\times {10}^{-4}$
$2^{-15}$	50	1.1425 $\times {10}^{-1}$	2.5201 $\times {10}^{-2}$	9.6557 $\times {10}^{-3}$	3.3489 $\times {10}^{-3}$	1.2836 $\times {10}^{-3}$	5.3330 $\times {10}^{-4}$	2.4725 $\times {10}^{-4}$
	100	1.1348 $\times {10}^{-1}$	3.1561 $\times {10}^{-2}$	1.2809 $\times {10}^{-2}$	5.7077 $\times {10}^{-3}$	2.8214 $\times {10}^{-3}$	1.4496 $\times {10}^{-3}$	7.4301 $\times {10}^{-4}$
$D_x^{N,k}$	50	1.1425 $\times {10}^{-1}$	3.1740 $\times {10}^{-2}$	1.2410 $\times {10}^{-2}$	4.0043 $\times {10}^{-3}$	1.2979 $\times {10}^{-3}$	5.3330 $\times {10}^{-4}$	2.4725 $\times {10}^{-4}$
	100	1.1348 $\times {10}^{-1}$	3.1561 $\times {10}^{-2}$	1.2809 $\times {10}^{-2}$	5.7077 $\times {10}^{-3}$	2.8214 $\times {10}^{-3}$	1.4496 $\times {10}^{-3}$	7.4301 $\times {10}^{-4}$
$r^{N,k}$	50	1.8478 $\times {10}^0$	1.3548 $\times {10}^0$	1.6319 $\times {10}^0$	1.6253 $\times {10}^0$	1.2832 $\times {10}^0$	1.1089 $\times {10}^0$	-
	100	1.8462 $\times {10}^0$	1.3010 $\times {10}^0$	1.1662 $\times {10}^0$	1.0165 $\times {10}^0$	9.6076 $\times {10}^{-1}$	9.6419 $\times {10}^{-1}$	-

, we list the order of convergence and the maximum point-wise error corresponding to Example 1. When M=50, the CPU run time for the error given in the table was $7.6844\times {10}^2$ s. Additionally, Figure 1 shows the spatial profile of the solution of the problem described in Example 1 at different times t, for a fixed value of $\epsilon$ and is a visual representation of the solution with a strong boundary layer at $x=1$, where the solution changed rapidly over a very small distance. This boundary layer arose due to the strong boundary condition $\widetilde{u_h^n}(1,t)=0$, which forced the solution to approach the value zero very quickly. The figure also shows the location of a strong interior layer at $x=\frac{1}{2}$, which caused a localized increase in the solution, which arose due to the line source term at the line $x=\frac{1}{2}$. Figure 2 shows the corresponding point-wise maximum error, that is, the maximum error decreased as N increased, irrespective of $\epsilon$.

Example 2.

We examined another one-dimensional parabolic partial differential equation with a line source, described by Equations (1)–(3). The given data for this problem were as follows:

$$\begin{array}{cc}\hfill & u(x,0)=0,x\in \Omega =(0,1),\hfill \\ \hfill & u(0,t)=0,u(1,t)=0,t\in [0,1],\hfill \\ \hfill & a\left(x\right)=\frac{5+x}{5+3x^2}+e^{-1/x^2};b\left(x\right)=1;f(x,t)=e^{-\frac{1}{x}}+\sqrt{t};g\left(t\right)=1;z=\frac{1}{2}.\hfill \end{array}$$

Table 2. Maximum error and order of convergence for Example 2 with the number M = 50,100.

$\mathit{\epsilon }$	M	N (Number of Grid Points)
↓	↓	${\mathbf{2}}^{\mathbf{6}}$	${\mathbf{2}}^{\mathbf{7}}$	${\mathbf{2}}^{\mathbf{8}}$	${\mathbf{2}}^{\mathbf{9}}$	${\mathbf{2}}^{\mathbf{10}}$	${\mathbf{2}}^{\mathbf{11}}$	${\mathbf{2}}^{\mathbf{12}}$
$2^{-6}$	50	2.1266 $\times {10}^{-1}$	1.3349 $\times {10}^{-1}$	7.0424 $\times {10}^{-2}$	3.3712 $\times {10}^{-2}$	1.5223 $\times {10}^{-2}$	7.1724 $\times {10}^{-3}$	3.4900 $\times {10}^{-3}$
	100	2.8109 $\times {10}^{-1}$	2.2380 $\times {10}^{-1}$	1.3696 $\times {10}^{-1}$	7.0987 $\times {10}^{-2}$	3.3312 $\times {10}^{-2}$	1.5960 $\times {10}^{-2}$	7.8452 $\times {10}^{-3}$
$2^{-7}$	50	3.5020 $\times {10}^{-1}$	2.2511 $\times {10}^{-1}$	1.2890 $\times {10}^{-1}$	6.7146 $\times {10}^{-2}$	3.3507 $\times {10}^{-2}$	1.6344 $\times {10}^{-2}$	7.8856 $\times {10}^{-3}$
	100	4.5196 $\times {10}^{-1}$	3.8586 $\times {10}^{-1}$	2.6099 $\times {10}^{-1}$	1.4770 $\times {10}^{-1}$	7.6575 $\times {10}^{-2}$	3.8040 $\times {10}^{-2}$	1.8516 $\times {10}^{-2}$
$2^{-8}$	50	4.4469 $\times {10}^{-1}$	3.0061 $\times {10}^{-1}$	1.7211 $\times {10}^{-1}$	9.0480 $\times {10}^{-2}$	4.6009 $\times {10}^{-2}$	2.3068 $\times {10}^{-2}$	1.1432 $\times {10}^{-2}$
	100	5.8200 $\times {10}^{-1}$	5.2658 $\times {10}^{-1}$	3.6716 $\times {10}^{-1}$	2.1077 $\times {10}^{-1}$	1.1126 $\times {10}^{-1}$	5.6420 $\times {10}^{-2}$	2.8120 $\times {10}^{-2}$
$2^{-9}$	50	5.0400 $\times {10}^{-1}$	3.4278 $\times {10}^{-1}$	1.9902 $\times {10}^{-1}$	1.0568 $\times {10}^{-1}$	5.3720 $\times {10}^{-2}$	2.7130 $\times {10}^{-2}$	1.3529 $\times {10}^{-2}$
	100	6.6504 $\times {10}^{-1}$	6.0947 $\times {10}^{-1}$	4.3599 $\times {10}^{-1}$	2.5668 $\times {10}^{-1}$	1.3564 $\times {10}^{-1}$	6.9369 $\times {10}^{-2}$	3.4832 $\times {10}^{-2}$
$2^{-10}$	50	5.3763 $\times {10}^{-1}$	3.6934 $\times {10}^{-1}$	2.1445 $\times {10}^{-1}$	1.1366 $\times {10}^{-1}$	5.8382 $\times {10}^{-2}$	2.9379 $\times {10}^{-2}$	1.4749 $\times {10}^{-2}$
	100	7.1658 $\times {10}^{-1}$	6.5915 $\times {10}^{-1}$	4.7996 $\times {10}^{-1}$	2.8410 $\times {10}^{-1}$	1.5157 $\times {10}^{-1}$	7.7503 $\times {10}^{-2}$	3.8981 $\times {10}^{-2}$
$2^{-11}$	50	5.5558 $\times {10}^{-1}$	3.8344 $\times {10}^{-1}$	2.2246 $\times {10}^{-1}$	1.1819 $\times {10}^{-1}$	6.0637 $\times {10}^{-2}$	3.0661 $\times {10}^{-2}$	1.5383 $\times {10}^{-2}$
	100	7.4477 $\times {10}^{-1}$	6.8558 $\times {10}^{-1}$	5.0407 $\times {10}^{-1}$	3.0050 $\times {10}^{-1}$	1.6044 $\times {10}^{-1}$	8.2210 $\times {10}^{-2}$	4.1429 $\times {10}^{-2}$
$2^{-12}$	50	5.6486 $\times {10}^{-1}$	3.9071 $\times {10}^{-1}$	2.2677 $\times {10}^{-1}$	1.2054 $\times {10}^{-1}$	6.1870 $\times {10}^{-2}$	3.1290 $\times {10}^{-2}$	1.5730 $\times {10}^{-2}$
	100	7.6104 $\times {10}^{-1}$	6.9976 $\times {10}^{-1}$	5.1666 $\times {10}^{-1}$	3.0910 $\times {10}^{-1}$	1.6535 $\times {10}^{-1}$	8.4716 $\times {10}^{-2}$	4.2767 $\times {10}^{-2}$
$2^{-13}$	50	5.6958 $\times {10}^{-1}$	3.9441 $\times {10}^{-1}$	2.2904 $\times {10}^{-1}$	1.2175 $\times {10}^{-1}$	6.2500 $\times {10}^{-2}$	3.1619 $\times {10}^{-2}$	1.5897 $\times {10}^{-2}$
	100	7.6927 $\times {10}^{-1}$	7.0691 $\times {10}^{-1}$	5.2310 $\times {10}^{-1}$	3.1350 $\times {10}^{-1}$	1.6788 $\times {10}^{-1}$	8.6049 $\times {10}^{-2}$	4.3442 $\times {10}^{-2}$
$2^{-14}$	50	5.7196 $\times {10}^{-1}$	3.9627 $\times {10}^{-1}$	2.3017 $\times {10}^{-1}$	1.2236 $\times {10}^{-1}$	6.2816 $\times {10}^{-2}$	3.1785 $\times {10}^{-2}$	1.5984 $\times {10}^{-2}$
	100	7.7341 $\times {10}^{-1}$	7.1050 $\times {10}^{-1}$	5.2636 $\times {10}^{-1}$	3.1573 $\times {10}^{-1}$	1.6915 $\times {10}^{-1}$	8.6729 $\times {10}^{-2}$	4.3793 $\times {10}^{-2}$
$2^{-15}$	50	5.7315 $\times {10}^{-1}$	3.9720 $\times {10}^{-1}$	2.3074 $\times {10}^{-1}$	1.2267 $\times {10}^{-1}$	6.2975 $\times {10}^{-2}$	3.1868 $\times {10}^{-2}$	1.6027 $\times {10}^{-2}$
	100	7.7548 $\times {10}^{-1}$	7.1230 $\times {10}^{-1}$	5.2803 $\times {10}^{-1}$	3.1685 $\times {10}^{-1}$	1.6980 $\times {10}^{-1}$	8.7074 $\times {10}^{-2}$	4.3971 $\times {10}^{-2}$
$D_x^{N,k}$	50	5.7315 $\times {10}^{-1}$	3.9720 $\times {10}^{-1}$	2.3074 $\times {10}^{-1}$	1.2267 $\times {10}^{-1}$	6.2975 $\times {10}^{-2}$	3.1868 $\times {10}^{-2}$	1.6027 $\times {10}^{-2}$
	100	7.7548 $\times {10}^{-1}$	7.1230 $\times {10}^{-1}$	5.2803 $\times {10}^{-1}$	3.1685 $\times {10}^{-1}$	1.6980 $\times {10}^{-1}$	8.7074 $\times {10}^{-2}$	4.3971 $\times {10}^{-2}$
$r^{N,k}$	50	5.2904 $\times {10}^{-1}$	7.8358 $\times {10}^{-1}$	9.1156 $\times {10}^{-1}$	9.6189 $\times {10}^{-1}$	9.8265 $\times {10}^{-1}$	9.9160 $\times {10}^{-1}$	-
	100	1.2260 $\times {10}^{-1}$	4.3187 $\times {10}^{-1}$	7.3683 $\times {10}^{-1}$	8.9999 $\times {10}^{-1}$	9.6349 $\times {10}^{-1}$	9.8568 $\times {10}^{-1}$	-

presents the maximum point-wise error and the associated convergence rate for Example 2, when M = 50. The CPU run time for the error given in the table was $7.8426\times {10}^2$ s. Additionally, as illustrated in Example 1, the graphs of numerical solution and point-wise maximum error of the problem are illustrated in Figure 3 and Figure 4, respectively.

From the above two Examples, 1 and 2, we see that the solutions of the problems exhibited an interior layer at $x=\frac{1}{2}$ and a boundary layer at $x=1.$

Note: For the numerical computation, a system with the following configuration i7 processor, 8.00 GB RAM was used.

7. Conclusions

In this study, we analyzed 1D parabolic equations that were singularly perturbed and contained a discontinuous source term on an interior line. Our objective was to discretize the time derivative using the Euler backward method and apply SDFEM to the locally one-dimensional stationary problems on Shishkin mesh. We aimed to determine the accuracy order of the spatial arrangement and evaluate test problems in terms of their maximum pointwise errors.

Our results showed that the accuracy order of the spatial arrangement was of second-order nature, but, due to the presence of the first-order term $\Delta t$ in the error bound, the overall accuracy order was limited to first-order accuracy. These findings are presented as Theorems 1 and 2. Examples 1 and 2 represent the test problems, and Figure 1 and Figure 3 depict their solutions. The layer occurred at the interior line $x=z$ and boundary line $x=1$, which is evident from these figures. Figure 2 and Figure 4 represent the convergence of the numerical solutions, that is, they illustrate that a higher value of N corresponded to a lower maximum pointwise error.

We evaluated the test problems in terms of their maximum pointwise errors, presented in Table 1 and Table 2. From the tables we see that the computational order of convergence was almost one, but Theorem 2 showed second-order convergence in space and first-order convergence in time. This was due to the first-order term in the final result. Our results demonstrated that the maximum pointwise error stabilized as parameter $ϵ$ decreased and decreased as parameter N increased.

It is worth noting that the results for SDFEM for ordinary differential equations, together with discontinuous source terms, are already available in the literature, as in, for example, [18], and, for the parabolic PDE without source line, in [3]. However, our results in this paper extend SDFEM to a parabolic PDE with a line source term, which is a significant contribution to the field. The error estimates are presented using the above-defined norm.

In summary, our study successfully achieved the objective of extending SDFEM to a parabolic PDE with a line source term. We determined the accuracy order of the spatial arrangement, evaluated the test problems in terms of their maximum pointwise errors, and demonstrated the stability of the maximum pointwise error as the parameter $ϵ$ decreases.