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

**R. Soundararajan 1,2, V. Subburayan 2,\* and Patricia J. Y. Wong <sup>3</sup>**


**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.

**Keywords:** singularly perturbed problem; parabolic differential equation; convection–diffusion problem; line discontinuous source term; streamline–diffusion finite element method; Shishkin mesh; uniformly convergent

**MSC:** 34K26; 35B25; 65M22; 65M50; 65N22

#### **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–5]. Such problems, but with non-smooth data functions, can be seen in [6–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

**Citation:** Soundararajan, R.; Subburayan, V.; Wong, P. J.Y. Streamline Diffusion Finite Element Method for Singularly Perturbed 1D-Parabolic Convection Diffusion Differential Equations with Line Discontinuous Source. *Mathematics* **2023**, *11*, 2034. https://doi.org/ 10.3390/math11092034

Academic Editor: Gabriel Eduard Vilcu

Received: 22 March 2023 Revised: 21 April 2023 Accepted: 24 April 2023 Published: 25 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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 *ε* 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 initialboundary 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 Sections 3 and 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

$$\mathfrak{M}u := u\_t - \varepsilon u\_{xx} + a(\mathbf{x})u\_x + b(\mathbf{x})u = \delta(\mathbf{x} - z)\lg(t) + f(\mathbf{x}, t), \ (\mathbf{x}, t) \in \Omega^\* \times (0, T], \tag{1}$$

$$
u(\mathbf{x},0) = \mu\_0(\mathbf{x}), \ \mathbf{x} \in \Omega,\tag{2}$$

$$
u(0,t) = 0 = \boldsymbol{u}(1,t), \; t \in [0,T],\tag{3}$$

where 0 <sup>&</sup>lt; *<sup>ε</sup>* 1 is a very small positive parameter, <sup>Ω</sup><sup>∗</sup> <sup>=</sup> <sup>Ω</sup><sup>−</sup> <sup>∪</sup> <sup>Ω</sup>+, <sup>Ω</sup><sup>−</sup> = (0, *<sup>z</sup>*), <sup>Ω</sup><sup>+</sup> = (*z*, 1), <sup>Ω</sup> = [0, 1], the functions *<sup>a</sup>*(*x*) <sup>≥</sup> *<sup>α</sup>* <sup>&</sup>gt; 0, *<sup>b</sup>*(*x*) <sup>&</sup>gt; *<sup>β</sup>* <sup>&</sup>gt; 0, *<sup>g</sup>*(*t*), *<sup>f</sup>*(*x*, *<sup>t</sup>*) and *<sup>u</sup>*0(*x*) are sufficiently differentiable and bounded in their respective domains, *δ*(*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 *ε*, or the discretization parameters, such as *<sup>N</sup>* or *<sup>M</sup>*. For practical purposes, the accepted convention is to assume *<sup>ε</sup>* <sup>≤</sup> *CN*−<sup>1</sup> for the convection coefficient problem. Additionally, the authors use the supremum norm <sup>|</sup>*ψ*|*<sup>D</sup>* <sup>=</sup> sup**x**∈*<sup>D</sup>* <sup>|</sup>*ψ*(**x**)<sup>|</sup> 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 Δ*t*, such that Ω*<sup>M</sup> <sup>t</sup>* <sup>=</sup> {*ti* <sup>=</sup> *<sup>i</sup>*Δ*t*}*<sup>M</sup> <sup>i</sup>*=0, <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>T</sup> <sup>M</sup>* , 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 *tn* in the set Ω*<sup>M</sup> <sup>t</sup>* , where *n* ranges from 1 to *M*:

$$(I + \Delta t L)\bar{u}^n = \bar{u}^{n-1} + \Delta t [f(\mathbf{x}, t\_n) + \delta(\mathbf{x} - z)g(t\_n)], \ \mathbf{x} \in \Omega^\* \tag{4}$$

$$
\mathfrak{u}^n(0) = 0 = \mathfrak{u}^n(1),\tag{5}
$$

where the differential operator *<sup>L</sup>* <sup>=</sup> <sup>−</sup>*<sup>ε</sup> <sup>d</sup>*<sup>2</sup> *dx*<sup>2</sup> <sup>+</sup> *<sup>a</sup>*(*x*) *<sup>d</sup> dx* + *<sup>b</sup>*(*x*) and *<sup>u</sup>*˜*<sup>n</sup>* is the solution of (4) and (5). Note that *<sup>u</sup>*˜0(*x*) = <sup>0</sup> <sup>=</sup> *<sup>u</sup>*0(*x*), *<sup>x</sup>* <sup>∈</sup> <sup>Ω</sup>.

To achieve the stability of the scheme given by (4) and (5), one can easily prove that the operator Q = (*I* + Δ*tL*) satisfies the maximum principle:

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

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

$$\begin{cases} L^\* \bar{u}^n := -\varepsilon \frac{d^2 \bar{u}^n(\mathbf{x})}{d^2 \mathbf{x}} + a(\mathbf{x}) \frac{d \bar{u}^n(\mathbf{x})}{d \mathbf{x}} + c(\mathbf{x}) \bar{u}^n(\mathbf{x}) = \mathcal{g}^n(\mathbf{x}), \ \mathbf{x} \in \Omega^\*,\\ \bar{u}^n(0) = 0 = \bar{u}^n(1) \end{cases} \tag{7}$$

where *c*(*x*) = *b*(*x*) + <sup>1</sup> <sup>Δ</sup>*<sup>t</sup>* and *<sup>g</sup>n*(*x*) = <sup>1</sup> <sup>Δ</sup>*<sup>t</sup> <sup>u</sup>*˜*n*−<sup>1</sup> <sup>+</sup> { *<sup>f</sup>*(*x*, *tn*) + *<sup>δ</sup>*(*<sup>x</sup>* <sup>−</sup> *<sup>z</sup>*)*g*(*tn*)}. The above scheme (7) is an ordinary differential equation in space variable *x* for each time step *tn*.

#### *2.3. Maximum Principle and Derivative Estimates*

**Lemma 1** (Maximum Principle)**.** *Suppose that there exists a function ξ, belonging to the set <sup>C</sup>*0(Ω) <sup>∩</sup> *<sup>C</sup>*2(Ω∗)*, that satisfies the following conditions: <sup>ξ</sup>*(*x*) <sup>≥</sup> <sup>0</sup> *for <sup>x</sup>* <sup>=</sup> 0, 1, *<sup>L</sup>*∗*ξ*(*x*) <sup>≥</sup> <sup>0</sup> *for all x* ∈ Ω<sup>∗</sup> *and ξ* (*z*−) − *ξ* (*z*+) <sup>≥</sup> <sup>0</sup> *Then, it can be concluded that <sup>ξ</sup>*(*x*) <sup>≥</sup> <sup>0</sup> *for all x* <sup>∈</sup> <sup>Ω</sup>*.*

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

**Lemma 2.** *Suppose u*˜*<sup>n</sup> is the solution of problem (7) and is decomposed as u*˜*<sup>n</sup>* = *r* + *s. Then, the derivatives of the regular components satisfy*

$$\left\| \frac{d^l r(\mathbf{x})}{d \mathbf{x}^l} \right\|\_{\Omega^\*} \le \mathbb{C} (1 + \varepsilon^{-l+2}), \ l = 0, 1, 2, 3, \tag{8}$$

*and, the derivatives of the singular components satisfy*

$$\left|\frac{d^l s(\mathbf{x})}{d \mathbf{x}^l}\right| \le C \varepsilon^{-l} \begin{cases} \exp\left(\frac{a(\mathbf{x} - \boldsymbol{z})}{\varepsilon}\right), \ \mathbf{x} \in \Omega\_-, \ l = 0, 1, 2, 3, \\\exp\left(\frac{a(\mathbf{x} - \boldsymbol{z})}{\varepsilon}\right), \ \mathbf{x} \in \Omega\_+. \end{cases} \tag{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,* <sup>|</sup> *<sup>∂</sup><sup>i</sup> u <sup>∂</sup>t<sup>i</sup>* | ≤ *<sup>C</sup>*, *<sup>i</sup>* = 0, 1, 2, 3, *then en*<sup>∞</sup> ≤ *<sup>C</sup>*Δ*<sup>t</sup>* <sup>2</sup> *and En*<sup>∞</sup> <sup>≤</sup> *<sup>C</sup>*Δ*t, where en* <sup>=</sup> *<sup>u</sup>*(*x*, *tn*) <sup>−</sup> *<sup>u</sup>*¯*<sup>n</sup> and En* <sup>=</sup> <sup>∑</sup>*<sup>n</sup> <sup>j</sup>*=<sup>1</sup> *ej.*

**Proof.** As *u*˜*<sup>n</sup>* is the solution of (4) and (5), we have

$$t\left\{(I+\Delta tL)\ddot{u}^n-\Delta t\{f(\mathbf{x},t\_n)+\delta(\mathbf{x}-z)\mathbf{g}(t\_n)\}=u^{n-1}.\tag{10}$$

By means of Taylor expansion, we have

$$u(\mathbf{x}, t\_{n-1}) = u(\mathbf{x}, t\_n) - \Delta t \frac{\partial u}{\partial t}(\mathbf{x}, t\_n) + O(\Delta t^2). \tag{11}$$

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

$$\frac{\partial u}{\partial t}(\mathbf{x}, t\_n) = \left\{ f(\mathbf{x}, t\_n) + \delta(\mathbf{x} - \mathbf{z}) \mathbf{g}(t\_n) \right\} - L u^n. \tag{12}$$

From (10)–(12), we have

$$(I + \Delta t L)e\_{\mathfrak{n}} = O(\Delta t^2), \ e\_{\mathfrak{n}}(0) = e\_{\mathfrak{n}}(1) = 0.$$

From this, we have *en* ≤ *<sup>C</sup>*(Δ*t*)<sup>2</sup> by using the stability result in (6).

The stability result and consistency property of (4) and (5) together implies *En* ≤ *C*Δ*t*.

#### **3. Weak Formulation**

The standard weak formulation of problem (7) for a fixed *n* is given as follows. Find *<sup>u</sup>*˜*<sup>n</sup>* <sup>∈</sup> *<sup>V</sup>* <sup>=</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω), such that

$$B(\mathfrak{u}^n, \upsilon) = \mathfrak{g}^n(\upsilon), \,\forall \upsilon \in V\_{\prime}$$

where

$$B(y,v) = \varepsilon(y',v') + (ay'+cy,v),$$

$$\mathbf{g}''(v) = (\mathbf{g}'',v).$$

Here, (. , .) represents the inner product in *L*2(Ω∗).

Consider the mesh in space as <sup>Ω</sup>*<sup>N</sup> <sup>x</sup>* = {*x*0, *x*1, ... , *xN*}, where N is some positive integer. We define a mesh that includes the point *z* as one of its nodes. Let *φ<sup>i</sup>* be the basis for the finite-dimensional subspace *Vh* of piece-wise linear polynomials. The basis is given by

$$\phi\_i(\mathfrak{x}) = \begin{cases} \frac{\mathfrak{x} - \mathfrak{x}\_{i-1}}{h\_i} & \mathfrak{x}\_{i-1} \le \mathfrak{x} \le \mathfrak{x}\_{i-1} \\\\ 0 & \text{otherwise} \end{cases}$$

where *hi* = *xi* − *xi*−<sup>1</sup> and *<sup>i</sup>* ∈ {1, 2, ... , *<sup>N</sup>*}. We use this basis to ensure that *<sup>z</sup>* is one of the mesh points.

The standard Galerkin method is given as follows: find *u*˜*<sup>n</sup> <sup>h</sup>* <sup>∈</sup> *Vh* <sup>⊂</sup> *<sup>V</sup>* <sup>=</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) such that

$$B\_h(\vec{u}\_h^n, v\_h) = g\_h^n(v\_h)\_{\prime} \,\,\forall v\_h \in V\_h.$$

#### **4. Streamline Diffusion Finite Element Formulation**

The streamline diffusion weak formulation for the problem (7) is to find *u*˜*<sup>n</sup> <sup>h</sup>* ∈ *Vh* ⊂ *V* = *H*<sup>1</sup> <sup>0</sup> (Ω), such that

$$B\_{\mathbf{h}}(\bar{\mu}\_{\mathbf{h}}^n, \upsilon\_{\mathbf{h}}) = \mathbb{g}\_{\mathbf{h}}^n(\upsilon\_{\mathbf{h}}), \,\forall \upsilon \in V\_{\mathbf{h}}.\tag{13}$$

where

$$\begin{aligned} \mathcal{B}\_{\hbar}(\boldsymbol{y}, \boldsymbol{w}) &= \varepsilon(\boldsymbol{y}', \boldsymbol{w}') + (\boldsymbol{a}\boldsymbol{y}' + \boldsymbol{c}\boldsymbol{y}, \boldsymbol{w}) + \sum\_{i=1}^{N} \int\_{\boldsymbol{x}\_{i-1}}^{\boldsymbol{x}\_{i}} \delta\_{i}(-\boldsymbol{\varepsilon}\boldsymbol{y}'' + \boldsymbol{a}\boldsymbol{y}' + \boldsymbol{c}\boldsymbol{y}) a \boldsymbol{w}' d\boldsymbol{x}\_{i} \\\\ \mathcal{g}\_{\hbar}^{\boldsymbol{n}}(\boldsymbol{w}) &= (\boldsymbol{\varrho}^{\boldsymbol{n}}, \boldsymbol{w}) + \sum\_{i=1}^{N} \int\_{\boldsymbol{x}\_{i-1}}^{\boldsymbol{x}\_{i}} \delta\_{i}\boldsymbol{\mathfrak{g}}^{\boldsymbol{n}}(\boldsymbol{x}) a(\boldsymbol{x}) \boldsymbol{w}' d\boldsymbol{x}. \end{aligned}$$

The parameter, known as the Streamline–diffusion (SD) parameter, denoted as *δi*, is decided later.

The SDFEM's relevant difference scheme is presented below.

$$L\_{SD}^{N} := -\varepsilon [D^{+}u\_{i} - D^{-}u\_{i}] + \mathfrak{a}\_{i}D^{-}u\_{i} + \mathfrak{f}\_{i}D^{+}u\_{i} + \gamma\_{i}u\_{i} = \mathfrak{g}\_{h}(\mathfrak{g}\_{i}),\tag{14}$$

where the symbols *D*<sup>+</sup> and *D*<sup>−</sup> are given by the following:

$$D^{+}\boldsymbol{u}\_{i} = \frac{\boldsymbol{u}\_{i+1} - \boldsymbol{u}\_{i}}{h\_{i+1}}, \ D^{-}\boldsymbol{u}\_{i} = \frac{\boldsymbol{u}\_{i} - \boldsymbol{u}\_{i-1}}{h\_{i}}.$$

and

$$\begin{split} a\_{i} &= \int\_{x\_{i-1}}^{x\_{i}} a(\mathbf{x}) \boldsymbol{\phi}\_{i} d\mathbf{x} - h\_{i} \int\_{x\_{i-1}}^{x\_{i}} c(\mathbf{x}) \boldsymbol{\phi}\_{i-1} \boldsymbol{\phi}\_{i} d\mathbf{x} + \frac{\delta\_{i}}{h\_{i}} \int\_{x\_{i-1}}^{x\_{i}} a(\mathbf{x})^{2} \, d\mathbf{x} \\ &- \delta\_{i} \int\_{x\_{i-1}}^{x\_{i}} a(\mathbf{x}) c(\mathbf{x}) \boldsymbol{\phi}\_{i-1} d\mathbf{x}, \\ \boldsymbol{\beta}\_{i} &= \int\_{x\_{i}}^{x\_{i+1}} a(\mathbf{x}) \boldsymbol{\phi}\_{i} d\mathbf{x} + h\_{i+1} \int\_{x\_{i}}^{x\_{i+1}} c(\mathbf{x}) \boldsymbol{\phi}\_{i+1} \boldsymbol{\phi}\_{i} d\mathbf{x} - \frac{\delta\_{i+1}}{h\_{i+1}} \int\_{x\_{i}}^{x\_{i+1}} a(\mathbf{x})^{2} \, d\mathbf{x} \\ &- \delta\_{i+1} \int\_{x\_{i}}^{x\_{i+1}} a(\mathbf{x}) c(\mathbf{x}) \boldsymbol{\phi}\_{i+1} d\mathbf{x}, \\ \boldsymbol{\gamma}\_{i} &= \int\_{x\_{i-1}}^{x\_{i+1}} c(\mathbf{x}) \boldsymbol{\phi}\_{i} d\mathbf{x} + \frac{\delta\_{i}}{h\_{i}} \int\_{x\_{i-1}}^{x\_{i}} a(\mathbf{x}) c(\mathbf{x}) d\mathbf{x} - \frac{\delta\_{i+1}}{h\_{i+1}} \int\_{x\_{i}}^{x\_{i+1}} a(\mathbf{x}) c(\mathbf{x}) d\mathbf{x}. \end{split}$$

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 Mmatrix. To satisfy this condition, the streamline–diffusion parameter *δ<sup>i</sup>* 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 \varepsilon\_{\prime\prime}$$

This can be written as *Cihi* <sup>≤</sup> *Ai* <sup>+</sup> *<sup>ε</sup>*, where *Ai* <sup>=</sup> *xi xi*−<sup>1</sup> *<sup>a</sup>*(*x*)*φidx* and *Ci* <sup>=</sup> *xi xi*−<sup>1</sup> *<sup>c</sup>*(*x*)*φi*−1*φidx*. If this condition is not met, *α<sup>i</sup>* is set to zero for *i* ∈ {1, 2, ... , *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{A} \overline{C}\_i}$$

where, *A*<sup>2</sup> *<sup>i</sup>* <sup>=</sup> *xi xi*−<sup>1</sup> (*a*(*x*))2*dx* and *ACi* <sup>=</sup> *xi xi*−<sup>1</sup> *<sup>a</sup>*(*x*)*c*(*x*)*φi*−1*dx*. If we summarize the above conditions, then we obtain

$$\delta\_i = \begin{cases} 0, & \text{if } \ \overline{\mathsf{C}}\_i h\_i \le \overline{A}\_i + \varepsilon \\ \frac{(h\_i \overline{\mathsf{C}}\_i - \overline{A}\_i) h\_i}{\overline{A}\_i^2 - h\_i \overline{\mathsf{A} \mathsf{C}\_i}}, & \text{otherwise.} \end{cases} \tag{15}$$

As a result, it is evident that *δ<sup>i</sup>* = *O*(*hi*). By choosing *vh* as *φ<sup>i</sup>* as well as *u*˜*<sup>h</sup>* = (*u*˜*<sup>n</sup> <sup>h</sup>* (*x*1), *u*˜*n <sup>h</sup>* (*x*2), ... , *<sup>u</sup>*˜*<sup>n</sup> <sup>h</sup>* (*xN*)) <sup>∈</sup> *<sup>R</sup>N*, we obtain the linear system of algebraic equations associated with the SDFEM scheme

$$\mathsf{K}\vec{u}\_{\mathsf{h}} = \mathsf{g}\_{\mathsf{h}\prime}$$

where K represents a tridiagonal matrix *K* = (*kij*) and *gh* = *g<sup>n</sup> <sup>h</sup>* (*φi*). Upon computing the coefficients,

*ki*,*i*−<sup>1</sup> <sup>=</sup> <sup>−</sup> *<sup>ε</sup>h*−<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>h</sup>*−<sup>1</sup> *i xi xi*−<sup>1</sup> *a*(*x*)*φidx* + *xi xi*−<sup>1</sup> *<sup>c</sup>*(*x*)*φi*−1*φidx* <sup>−</sup> *<sup>δ</sup>ih*−<sup>2</sup> *i xi xi*−<sup>1</sup> *a*(*x*)<sup>2</sup> *dx* + *δih*−<sup>1</sup> *i xi xi*−<sup>1</sup> *<sup>a</sup>*(*x*)*c*(*x*)*φi*−1*dx*, *ki*,*<sup>i</sup>* =*εh*−<sup>1</sup> *<sup>i</sup>* <sup>+</sup> *<sup>ε</sup>h*−<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>+</sup> *<sup>h</sup>*−<sup>1</sup> *i xi xi*−<sup>1</sup> *<sup>a</sup>*(*x*)*φidx* <sup>−</sup> *<sup>h</sup>*−<sup>1</sup> *i*+1 *xi*+<sup>1</sup> *xi c*(*x*)*φidx* + *xi xi*−<sup>1</sup> *c*(*x*)*φ*<sup>2</sup> *<sup>i</sup> dx* + *xi*+<sup>1</sup> *xi c*(*x*)*φ*<sup>2</sup> *<sup>i</sup> dx* + *<sup>δ</sup>ih*−<sup>2</sup> *i xi xi*−<sup>1</sup> *a*(*x*)<sup>2</sup> *dx* + *δi*+1*h*−<sup>2</sup> *i*+1 *xi*+<sup>1</sup> *xi a*(*x*)2*dx* + *δih*−<sup>1</sup> *i xi xi*−<sup>1</sup> *<sup>a</sup>*(*x*)*c*(*x*)*φidx* <sup>−</sup> *<sup>δ</sup>i*+1*h*−<sup>1</sup> *i*+1 *xi*+<sup>1</sup> *xi a*(*x*)*c*(*x*)*φ<sup>i</sup> dx ki*,*i*+<sup>1</sup> <sup>=</sup> <sup>−</sup> *<sup>ε</sup>h*−<sup>1</sup> *<sup>i</sup>*+<sup>1</sup> <sup>+</sup> *<sup>h</sup>*−<sup>1</sup> *i*+1 *xi*+<sup>1</sup> *xi a*(*x*)*φidx* + *xi*+<sup>1</sup> *xi c*(*x*)*φi*+1*φ<sup>i</sup> dx* <sup>−</sup> *<sup>δ</sup>i*+1*h*−<sup>2</sup> *i*+1 *xi*+<sup>1</sup> *xi <sup>a</sup>*(*x*)2*dx* <sup>−</sup> *<sup>δ</sup>i*+1*h*−<sup>1</sup> *i*+1 *xi*+<sup>1</sup> *xi a*(*x*)*c*(*x*)*φi*+1*dx*.

**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*∞, *w*−1,∞) stability of the SDFEM, we introduce the discrete Green's function. The i-th discrete Green's function, denoted as *<sup>λ</sup><sup>i</sup>* <sup>∈</sup> *Vh*, is determined by solving the following problem:

$$\begin{cases} B\_h(\phi\_{j\prime}\lambda^i) = \delta\_{ij\prime} \\ \lambda^i(0) = \lambda^i(1) = 0, \end{cases} \tag{16}$$

where *δij* denotes the Kronecker delta. Additional details on these and related concepts can be found in references such as [24–26].

**Lemma 4** ([25])**.** *If λ<sup>i</sup>* = ∑*<sup>N</sup> <sup>j</sup>*=<sup>1</sup> *λ<sup>i</sup> j φj*,*, then the discrete Green's function exhibits the following conditions:*

*(1) λ<sup>i</sup> <sup>j</sup>* ≥ 0, *i*, *j* ∈ {1, 2, . . . , *N*}, *(2)* <sup>0</sup> <sup>≤</sup> *<sup>λ</sup><sup>i</sup>* <sup>1</sup> <sup>&</sup>lt; ... <sup>&</sup>lt; *<sup>λ</sup><sup>i</sup> <sup>i</sup>* <sup>&</sup>gt; *<sup>λ</sup><sup>i</sup> <sup>i</sup>*+<sup>1</sup> ... <sup>&</sup>gt; *<sup>λ</sup><sup>i</sup> N*.

**Lemma 5.** *If we denote c*<sup>0</sup> = min*i*{*ci*}, *i* = 1, 2, ... , *N where ci* = *ki*+1,*<sup>i</sup>* − *ki*,*i*+<sup>1</sup> *and λ<sup>i</sup>* = ∑*<sup>N</sup> <sup>j</sup>*=<sup>1</sup> *λ<sup>i</sup> j φj*, *then λ<sup>i</sup> <sup>i</sup>* <sup>≤</sup> *<sup>c</sup>*−<sup>1</sup> <sup>0</sup> .

**Proof.** From (16), for any *vh* ∈ *Vh*, we have

$$B\_{\text{li}}(\upsilon\_{\text{h}\prime}\lambda^{\bar{i}}) = \upsilon\_{\text{li}}(x\_{\text{i}}).$$

Considering *vh* = ∑*<sup>i</sup> <sup>j</sup>*=<sup>1</sup> *φj*, one can derive

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

Expanding, we obtain

$$\begin{aligned} (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{aligned}$$

K is an M-matrix implying (*ki*,*i*−<sup>1</sup> + *ki*,*<sup>i</sup>* + *ki*,*i*+1) ≥ 0, for *<sup>i</sup>* ∈ {2, ... , *<sup>N</sup>* − <sup>1</sup>}. Now, by means of Lemma 4 and *<sup>k</sup>*1,1 <sup>&</sup>gt; *<sup>k</sup>*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 *λ<sup>i</sup> <sup>i</sup>* <sup>≤</sup> *<sup>c</sup>*−<sup>1</sup> <sup>0</sup> .

**Lemma 6** ([25])**.** *The SDFEM, equipped with the streamline–diffusion parameter specified in Equation (15), is uniformly stable in the* (*l*∞, *w*−1,∞) *norm, as demonstrated by the following result:*

$$||v\_h||\_{\infty} \le \frac{2}{c\_0} \max\_j \left| \sum\_{k=j}^N (Kv\_h)\_k \right|\_{\prime} \quad \forall v\_h \in V\_{h\prime} \ j = 1, \dots, N.$$

#### *4.2. Shishkin Type Mesh*

We used a general type mesh introduced in [27] with adapted layers at the points *<sup>x</sup>* <sup>=</sup> *<sup>z</sup>* and *<sup>x</sup>* <sup>=</sup> 1. Let N be a positive even integer greater than 4 and *<sup>σ</sup>*<sup>1</sup> <sup>=</sup> min{ *<sup>z</sup>* <sup>2</sup> , *<sup>τ</sup>*<sup>0</sup> *<sup>α</sup> ε*ln *N*}, *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> min{ <sup>1</sup>−*<sup>z</sup>* <sup>2</sup> , *<sup>τ</sup>*<sup>0</sup> *<sup>α</sup> <sup>ε</sup>*ln *<sup>N</sup>*}, *<sup>τ</sup>*<sup>0</sup> <sup>≥</sup> 2. Here, we considered *<sup>σ</sup>*<sup>1</sup> <sup>=</sup> *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> *<sup>τ</sup>*<sup>0</sup> *<sup>α</sup> ε*ln *N*. Let Ω*<sup>s</sup>* = (0, *z* − *σ*1) ∪ (*z*, 1 − *σ*2) and Ω<sup>0</sup> = (*z* − *σ*1, *z*) ∪ (1 − *σ*2, 1). On Ω*<sup>s</sup>* and Ω*o*, the mesh is equidistant and graded respectively. The transition points are chosen to be

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

We chose two mesh generating functions *ϕ*<sup>1</sup> and *ϕ*<sup>2</sup> for the particular layers:

$$\begin{aligned} \varphi\_1(\frac{1}{4}) &= \ln N\_\prime & \varphi\_1(\frac{1}{2}) &= 0, \\\\ \varphi\_2(\frac{3}{4}) &= \ln N\_\prime & \varphi\_2(1) &= 0. \end{aligned}$$

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

$$x\_{i} = \begin{cases} \frac{4i}{N}(z - \sigma\_{1}) & i = 0, \dots, \frac{N}{4}, \\ z - \frac{\pi\_{0}}{\pi} \varepsilon \varphi\_{1}(t\_{i}) & i = \frac{N}{4} + 1, \dots, \frac{N}{2}, \\ z + \frac{4(1 - z - \sigma\_{2})(i - \frac{N}{2})}{N} & i = \frac{N}{2} + 1, \dots, \frac{3N}{4}, \\ -\frac{\pi\_{0}}{\pi} \varepsilon \varphi\_{2}(t\_{i}) + 1 & i = \frac{3N}{4} + 1, \dots, N, \end{cases} \tag{17}$$

where *ti* = *<sup>i</sup> <sup>N</sup>* , *ϕ<sup>i</sup>* = − ln *ψ<sup>i</sup>* for *i* = 1, 2, *ψ*1(*t*) = exp (−(2 − 4*t*)ln *N*) and *ψ*2(*t*) = exp (−(4 − 4*t*)ln *N*) on Shishkin mesh (S-mesh).

max*<sup>i</sup>* |*ψ i* |, *i* = 1, 2 plays a part in error analysis as we see in the upcoming section. we have max |*ψ* | ≤ *C* ln *N* in the S-mesh. Assuming that the mesh generating function *ϕi*, for *i* = 1, 2, obeys the following:

$$\max\{|\varphi\_i'|\} \le N\mathbb{C}\_{0\prime} \text{ for } i = 1, 2. \tag{18}$$

We observe that on the coarse part Ω*s*, the following holds for the S-mesh:

$$h\_i N \le \mathbb{C}.\tag{19}$$

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

$$h\_i \le \mathcal{C} \varepsilon \frac{\ln N}{N} \tag{20}$$

and, by the assumption (18), we obtain

$$\frac{h\_i}{\varepsilon} \le CN^{-1} \max |\varrho'| \le C. \tag{21}$$

#### **5. Error Estimations**

The error at each time level *tn*, where 1 ≤ *n* ≤ *M* − 1, can be expressed as the difference between *u*˜*n*(*xi*) and *u*˜*<sup>n</sup> <sup>h</sup>* (*xi*) for all *xi* ∈ Ω∗. Another way to represent this error is by utilizing the linear interpolant *u*˜*<sup>n</sup> <sup>I</sup>* of *<sup>u</sup>*˜*<sup>n</sup>* as follows:

$$\begin{split} \boldsymbol{\sigma}(\boldsymbol{x}\_{i}) = (\boldsymbol{\tilde{u}}^{\boldsymbol{n}} - \boldsymbol{\tilde{u}}\_{h}^{\boldsymbol{n}})(\boldsymbol{x}\_{i}) &= (\boldsymbol{\tilde{u}}^{\boldsymbol{n}} - \boldsymbol{\tilde{u}}\_{I}^{\boldsymbol{n}})(\boldsymbol{x}\_{i}) + (\boldsymbol{\tilde{u}}\_{I}^{\boldsymbol{n}} - \boldsymbol{\tilde{u}}\_{h}^{\boldsymbol{n}})(\boldsymbol{x}\_{i}) \\ &= \boldsymbol{e}^{I} + \boldsymbol{e}^{D}. \end{split}$$

The errors of interpolation and discretization are denoted by *e<sup>I</sup>* and *eD*, 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

$$\mathfrak{e}^{D} = (\widetilde{u}\_I^n - \widetilde{u}\_h^n)(\mathfrak{x}) = \sum\_{i=0}^{N} \mathfrak{e}\_i \mathfrak{d}\_i. \tag{22}$$

By means of the orthogonality property, we have

$$B\_{\hbar}(\mathfrak{u}^n - \mathfrak{u}\_{\hbar'}^n v\_{\hbar}) = 0, \,\,\forall v\_{\hbar} \in V\_{\hbar}.$$

Using the above property, we obtain

$$\begin{aligned} B\_{\mathbb{H}}(e^D, \phi\_i) &= B\_{\mathbb{H}}(\vec{u}\_I^n - \vec{u}\_{h'}^n, \phi\_i) \\ &= B\_{\mathbb{H}}(\vec{u}\_I^n - \vec{u}^n, \phi\_i), \quad i \in \{1, 2, \dots, N - 1\}. \end{aligned}$$

Therefore, the equation for error is

$$\begin{cases} & B\_h(\boldsymbol{\varepsilon}^D, \boldsymbol{\phi}\_i) = B\_h(\boldsymbol{\tilde{n}}\_I^n - \boldsymbol{\tilde{n}}^n, \boldsymbol{\phi}\_i), \quad i \in \{1, 2, \dots, N - 1\}, \\ & e\_0^D = 0, \; e\_N^D = 0. \end{cases} \tag{23}$$

The following lemma gives the explicit expression for *Bh*(*u*˜*<sup>n</sup> <sup>I</sup>* <sup>−</sup> *<sup>u</sup>*˜*n*, *<sup>φ</sup>i*). **Lemma 7.** *The error Equation (23) involves the bilinear form Bh*(*u*˜*<sup>n</sup> <sup>I</sup>* <sup>−</sup> *<sup>u</sup>*˜*n*, *<sup>φ</sup>i*)*, which can be explicitly expressed as follows:*

$$B\_{\hbar}(\vec{u}\_{I}^{n} - \vec{u}^{n}, \phi\_{\bar{i}}) = O\_{\bar{i}} + (P\_{\bar{i}} + P\_{\bar{i}+1}) + (Q\_{\bar{i}} - Q\_{\bar{i}+1}) + (R\_{\bar{i}} - R\_{\bar{i}+1}) + (S\_{\bar{i}} - S\_{\bar{i}+1}),$$

*where*

$$O\_i = \int\_{x\_{i-1}}^{x\_{i+1}} a'(\mathbf{x}) (\overline{u}\_I^n - \overline{u}^n) \phi\_i d\mathbf{x},\tag{24}$$

$$P\_i = \int\_{x\_{i-1}}^{x\_i} c(\mathbf{x}) (\vec{u}\_I^n - \vec{u}^n) \phi\_i d\mathbf{x},\tag{25}$$

$$Q\_i = \frac{\delta\_i}{h\_i} \int\_{x\_{i-1}}^{x\_i} (a(\mathbf{x})c(\mathbf{x})(\tilde{u}\_I^n - \tilde{u}^n) + a(\mathbf{x})^2(\tilde{u}\_I^n - \tilde{u}^n)')d\mathbf{x},\tag{26}$$

$$R\_i = -\varepsilon \delta\_i \mathfrak{h}\_i^{-1} \int\_{\mathfrak{x}\_{i-1}}^{\mathfrak{x}\_i} a(\mathfrak{x}) (\mathfrak{u}\_I^{\mathfrak{n}} - \mathfrak{u}^{\mathfrak{n}})^{\prime\prime} d\mathfrak{x}\_{\prime} \tag{27}$$

$$S\_i = h\_i^{-1} \int\_{x\_{i-1}}^{x\_i} a(\mathbf{x}) (\tilde{u}\_I^n - \tilde{u}^n) d\mathbf{x}.\tag{28}$$

**Proof.** One can prove the required results by a simple calculation.

**Lemma 8.** *Assume that a*(*x*) <sup>∈</sup> *<sup>W</sup>*2,1, *then the following estimates are true:*

$$\begin{aligned} (i) \quad &\sum\_{i=1}^{N} |O\_{i}| \leq c\_{1} \|\vec{u}^{n} - \vec{u}\_{I}^{n}\|\_{\infty} \quad (ii) \sum\_{i=1}^{N} |P\_{i}| \leq c\_{2} \|\vec{u}^{n} - \vec{u}\_{I}^{n}\|\_{\infty} \\ (iii) \; ||Q\_{i}||\_{\infty} \leq c\_{3} \|\vec{u}^{n} - \vec{u}\_{I}^{n}\|\_{\infty} \quad (iv) \; ||R\_{i}||\_{\infty} \leq c\_{4} \|\|\vec{u}^{n} - \vec{u}\_{I}^{n}\|\_{\infty} \\ (v) \; ||S\_{i}||\_{\infty} \leq c\_{5} \|\vec{u}^{n} - \vec{u}\_{I}^{n}\|\_{\infty} .\end{aligned}$$

**Proof.** We prove this Lemma sequentially:

(i). From (24), we know that

$$O\_i = \int\_{\chi\_{i-1}}^{\chi\_{i+1}} a'(\boldsymbol{x}) (\bar{u}\_I^n - \bar{u}^n) \phi\_i d\boldsymbol{x}.$$

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

$$\begin{aligned} \sum\_{i=1}^N |O\_i| &\le ||\vec{u}\_I^n - \vec{u}^n||\_\infty \sum\_{i=1}^N \int\_{\substack{\vec{n}\_{i-1} \\ \alpha \le \|\vec{n}\_I^n - \vec{u}^n\|\_\infty |a(\mathbf{x})|\_{1,1}}}^N a'(\mathbf{x}) d\mathbf{x} \\ &= c\_1 ||\vec{u}\_I^n - \vec{u}^n||\_\infty . \end{aligned}$$

(ii). From Equation (25), we have

$$P\_i = \int\_{\mathfrak{x}\_{i-1}}^{\mathfrak{x}\_i} c(\mathfrak{x}) (\mathfrak{u}\_I^n - \mathfrak{u}^n) \phi\_i d\mathfrak{x}.$$

If we conduct a summation over all intervals, and derive

$$\sum\_{i=1}^{N} |P\_i| \le ||\vec{u}\_I^n - \vec{u}^n||\_\infty ||a(\mathbf{x})||\_\infty \, h\_i$$

$$= c\_2 ||\vec{u}\_I^n - \vec{u}^n||\_\infty.$$

(iii). From Equation (26), we have

$$Q\_i = \frac{\delta\_i}{\hbar\_i} \int\_{\chi\_{i-1}}^{\chi\_i} (a(\mathbf{x})c(\mathbf{x})(\vec{u}\_I^n - \vec{u}^n) + a(\mathbf{x})^2(\vec{u}\_I^n - \vec{u}^n)')d\mathbf{x}.$$

If the inequality *Cihi* ≤ *Ai* + *ε* holds, then *δ<sup>i</sup>* = 0 and, as a result, *ri* = 0. On the other hand, if the condition is not satisfied, then we have *δ<sup>i</sup>* = 0. In this scenario, it can be observed that *δ<sup>i</sup>* = *O*(*hi*), we obtain

$$\begin{aligned} |Q\_i| &\le C \|a(\boldsymbol{x})\|\_{\infty} \|\boldsymbol{\bar{u}}\_I^n - \boldsymbol{\bar{u}}^n\|\_{\infty} (\|\boldsymbol{c}(\boldsymbol{x})\|\_{\infty} h\_i - 2 |a(\boldsymbol{x})|\_{1,1}) \\ &= c\_3 \|\boldsymbol{\bar{u}}\_I^n - \boldsymbol{\bar{u}}^n\|\_{\infty} .\end{aligned}$$

(iv). From Equation (27), we have

$$R\_i = -\frac{\varepsilon \delta\_i}{h\_i} \int\_{x\_{i-1}}^{x\_i} a(\boldsymbol{x}) (\vec{u}\_I^n - \vec{u}^n)^{\prime\prime} d\boldsymbol{x},$$

and, here also we have *ri* = 0 when *Cihi* ≤ *Ai* + *ε*, and if it is not true, we have

$$\begin{aligned} |\mathcal{R}\_i| &\le \mathcal{C} |a(\mathfrak{x})|\_{2,1} ||\mathfrak{u}\_I^n - \mathfrak{u}^n||\_{\infty} \\ &= c\_4 ||\mathfrak{u}\_I^n - \mathfrak{u}^n||\_{\infty} .\end{aligned}$$

(v). Now, finally from Equation (28), we have

$$S\_{\bar{\imath}} = \frac{1}{h\_{\bar{\imath}}} \int\_{x\_{\bar{\imath}-1}}^{x\_{\bar{\imath}}} a(x) (\bar{u}\_I^n - \bar{u}^n) dx.$$

Simplifying this, we obtain

$$\begin{aligned} |S\_i| &\le ||a(\mathfrak{x})||\_\infty ||\mathfrak{u}\_I^n - \mathfrak{u}^n||\_\infty \\ &= c\_5 ||\mathfrak{u}\_I^n - \mathfrak{u}^n||\_\infty. \end{aligned}$$

**Lemma 9.** *Consider the problem (7) and let u*˜*<sup>n</sup> be the exact solution to this problem. Furthermore, let u*˜*<sup>n</sup> <sup>I</sup> denote the interpolant of <sup>u</sup>*˜*<sup>n</sup> on a given grid and let <sup>u</sup>*˜*<sup>n</sup> <sup>h</sup> represent the approximate solution at the time level tn. Then, It follows that*

$$\left\| \left\| \vec{u}^{n} - \vec{u}\_{h}^{n} \right\| \right\|\_{\infty} \leq C \left\| \left\| \vec{u}\_{I}^{n} - \vec{u}^{n} \right\| \right\|\_{\infty}.$$

**Proof.** Using Lemmas 6 and 8, we have

$$\begin{split} \|\boldsymbol{\bar{u}}\_{I}^{n} - \boldsymbol{\bar{u}}\_{h}^{n}\|\_{\infty} &\leq \mathbb{C} \max\_{i \in \{1, \ldots, N\}} \left| \sum\_{j=i}^{N} (\mathcal{K}(\boldsymbol{\bar{u}}\_{I}^{n} - \boldsymbol{\bar{u}}\_{h}^{n}))\_{j} \right| \\ &= \mathbb{C} \max\_{i \in \{1, \ldots, N\}} \left| \sum\_{j=i}^{N} \mathcal{B}\_{h}(\boldsymbol{\bar{u}}\_{I}^{n} - \boldsymbol{\bar{u}}^{n}, \boldsymbol{\phi}\_{j}) \right| \\ &\leq \mathbb{C} \left( \left| \sum\_{i=1}^{N} \mathcal{O}\_{i} \right| + 2 \left| \sum\_{i=1}^{N} \mathcal{P}\_{i} \right| + \max\_{i \in \{1, \ldots, N\}} (\mathcal{Q}\_{i} - \mathcal{Q}\_{N}) + \max\_{i \in \{1, \ldots, N\}} (\mathcal{R}\_{i} - \mathcal{R}\_{N}) \right) \\ &\quad + \mathbb{C} \max\_{i \in \{1, \ldots, N\}} (\mathcal{S}\_{i} - \mathcal{S}\_{N}) \\ &\leq \mathbb{C} \left( \sum\_{i=1}^{N} |\mathcal{O}\_{i}| + 2 \sum\_{i=1}^{N} |\mathcal{P}\_{i}| + 2 (\|\mathcal{Q}\_{i}\|\_{\infty} + \|\mathcal{R}\_{i}\|\_{\infty} + \|\mathcal{S}\_{i}\|\_{\infty}) \right) \\ &\leq \mathbb{C} \|\boldsymbol{\bar{u}}^{n}\_{I} - \boldsymbol{\bar{u}}^{n}\|\_{\infty}. \end{split}$$

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

$$||\mathfrak{u}^n - \mathfrak{u}\_{\hbar}^n||\_{\infty} \le C|\mathfrak{u}\_I^n - \mathfrak{u}^n||\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 u*˜*<sup>n</sup> be the classical solution of the problem (7), u*˜*<sup>n</sup> <sup>I</sup> be its interpolant on a Shishkin mesh with the grid points xi as in (17) and <sup>ε</sup>* <sup>≤</sup> *<sup>C</sup> <sup>N</sup>* . *Then, the following holds:*

$$||\mathfrak{u}\_I^n - \mathfrak{u}^n||\_{\infty} \le \begin{cases} C(N^{-1} \ln N)^2, & x \in \Omega\_{0, \epsilon} \\ C N^{-2}, & x \in \Omega\_{\epsilon} \end{cases}$$

**Proof.** We perform a separate analysis of the error due to interpolation on the domain Ω<sup>−</sup> as well as Ω+.

First, consider *<sup>x</sup>* in the domain <sup>Ω</sup>−. Let *<sup>u</sup>*˜*<sup>n</sup>* <sup>1</sup> (*x*) be the solution to the problem:

$$\begin{cases} L^\* \ddot{u}\_1^n(\varkappa) = \mathcal{g}^n(\varkappa), & \varkappa \in \Omega\_- = (0, z), \\ \ddot{u}^0(0) = 0, \; \dot{u}^0(z) = d. \end{cases}$$

For now, assume that d is a constant. Then, the solution *u*˜*<sup>n</sup>* <sup>1</sup> (*x*) exists and we can decompose it as a sum of two functions *r*<sup>1</sup> and *s*1, where *r*<sup>1</sup> and *s*<sup>1</sup> satisfy bounds (8) and (9). Then, we can write

$$
\tilde{u}\_1^{\mathfrak{n}}(\mathfrak{x}) - \tilde{u}\_{1,I}^{\mathfrak{n}}(\mathfrak{x}) = r\_1(\mathfrak{x}) - r\_{1,I}(\mathfrak{x}) + s\_1(\mathfrak{x}) - s\_{1,I}(\mathfrak{x}).
$$

By means of classical theory, and from (19) and (8), we have the following estimate for the interpolation error *r*1(*x*) − *r*1,*I*(*x*) on the regular part:

$$|r\_1(\mathbf{x}) - r\_{1,I}(\mathbf{x})| \le C h\_i^2 \max\_{x\_{i-1} \le x \le x\_i} |r\_1''| \le C h\_i^2 \le C N^{-2}.$$

For the interpolation error *r*1(*x*) − *r*1,*I*(*x*) in the layer part, one can derive:

$$\begin{aligned} |r\_1(\mathbf{x}) - r\_{1,l}(\mathbf{x})| &\leq \mathcal{C} \varepsilon^2 (N^{-1} \ln N)^2 \exp\left(\frac{-2\kappa}{\tau\_0 \varepsilon} (\mathbf{x}\_{i-1} - z)\right), \\ &\leq \mathcal{C} (N^{-1} \ln N)^2. \end{aligned}$$

The above expression is obtained by utilizing the inequality (20), selecting the transition point *<sup>σ</sup>*1, and assuming the condition *<sup>ε</sup>* <sup>≤</sup> *<sup>C</sup> <sup>N</sup>* holds. Now, for the interpolation error *s*1(*x*) − *s*1,*I*(*x*) on the regular part, we have

$$\begin{split} |s\_1(\boldsymbol{x}) - s\_{1,l}(\boldsymbol{x})| &\leq 2|s\_1(\boldsymbol{x})| \leq C \max\_{\boldsymbol{x}\_{i-1} \leq \boldsymbol{x} \leq \boldsymbol{x}\_i} \exp\left(-\frac{\alpha}{\varepsilon}(\boldsymbol{z} - \boldsymbol{x})\right) \\ &\leq C \exp\left(-\frac{\alpha}{\varepsilon}\sigma\_1\right) \\ &= C N^{-\tau\_0}. \end{split}$$

Regarding the interpolation error for the layer component, i.e., *s*1(*x*) − *s*1,*I*(*x*), we apply the classical theory to obtain the following:

$$\begin{aligned} |s\_1(\boldsymbol{x}) - s\_{1,l}(\boldsymbol{x})| &\leq Ch\_l^2 \max\_{x\_{i-1} \leq \boldsymbol{x} \leq x\_i} |s\_1''|\\ &\leq C(N^{-1} \ln N)^2 \max\_{x\_{i-1} \leq \boldsymbol{x} \leq x\_i} \exp\left(\frac{-\alpha}{\varepsilon}(\boldsymbol{z} - \boldsymbol{x})\right),\\ &\leq C(N^{-1} \ln N)^2. \end{aligned}$$

Here, we used the inequality (20), the first inequality of (9) and (21). Similarly, we can prove the case *x* ∈ Ω+.

**Theorem 2.** *Suppose that u*(*x*, *t*) *represents the classical solution of problem (1)–(3), u*˜*<sup>n</sup> <sup>h</sup>* (*x*) *represents the numerical solution of the SDFEM scheme (fully discrete) (13) at the nth time step tn. Under the conditions of Lemma 3 and Theorem 1, we have the error estimate as follows:*

$$\|\|\tilde{u}\_h^n(\mathbf{x}) - u(\mathbf{x}, t\_n)\|\|\_{\infty} \le \begin{cases} C((N^{-1} \ln N)^2 + \Delta t), & \mathbf{x} \in \Omega\_{0, \mathbf{x}}, \\ C(N^{-2} + \Delta t), & \mathbf{x} \in \Omega\_{\mathbf{x}}. \end{cases}$$

**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*2*N*,*k*/2(*x<sup>n</sup> <sup>i</sup>* , *ti*) on a grid Ω2*<sup>N</sup> <sup>x</sup>* <sup>×</sup> <sup>Ω</sup>2*<sup>M</sup> <sup>t</sup>* , where the spatial direction is divided into 2*N* intervals, while the temporal direction is divided into 2*M* intervals. The mesh Ω2*<sup>N</sup> <sup>x</sup>* <sup>×</sup> <sup>Ω</sup>2*<sup>M</sup> <sup>t</sup>* is obtained by dividing each segment of the previous mesh Ω*<sup>N</sup> <sup>x</sup>* <sup>×</sup> <sup>Ω</sup>*<sup>M</sup> <sup>t</sup>* , 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 *ε*.

$$\begin{aligned} E\_{\varepsilon}^{N,k} &= \max\_{i,\mathfrak{n}} \left| Y^{N,k}(\mathfrak{x}\_{i}^{\mathfrak{n}}, t\_{\mathfrak{n}}) - Y^{2N,k/2}(\mathfrak{x}\_{i}^{\mathfrak{n}}, t\_{\mathfrak{n}}) \right|, \\\\ D\_{\mathfrak{x}}^{N,k} &= \max\_{\mathfrak{c}} E\_{\mathfrak{c}}^{N,k}, \ r^{N,k} = \log\_2 \left( \frac{D\_{\mathfrak{x}}^{N,k}}{D\_{\mathfrak{x}}^{2Nk/2}} \right). \end{aligned}$$

**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{aligned} u(\mathbf{x},0)&=0, \; \mathbf{x} \in \Omega = (0,1),\\ u(0,t)&=0, \; u(1,t) = 0, \; t \in [0,1],\\ u(\mathbf{x})&=\mathbf{x}+\mathbf{5}; \; b(\mathbf{x}) = 1; \; f(\mathbf{x},t) = \mathbf{x}+\mathbf{t}; \; g(\mathbf{t}) = 1; \; z = \frac{1}{2}.\end{aligned}$$

In Table 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 <sup>×</sup> <sup>10</sup><sup>2</sup> 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 *ε* 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 *u*˜*<sup>n</sup> <sup>h</sup>* (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* = <sup>1</sup> <sup>2</sup> , which caused a localized increase in the solution, which arose due to the line source term at the line *x* = <sup>1</sup> 2 . Figure 2 shows the corresponding point-wise maximum error, that is, the maximum error decreased as *N* increased, irrespective of *ε*.


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

**Figure 1.** Numerical solution of Example 1.

**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{aligned} u(\mathbf{x},0)&=0, \; \mathbf{x} \in \Omega = (0,1),\\ u(0,t)&=0, \; u(1,t) = 0, \; t \in [0,1],\\ u(\mathbf{x})&= \frac{5+\mathbf{x}}{5+3\mathbf{x}^2} + e^{-1/\mathbf{x}^2}; \; b(\mathbf{x}) = 1; \; f(\mathbf{x},t) = e^{-\frac{1}{\mathbf{x}}} + \sqrt{t}; \; g(t) = 1; \; z = \frac{1}{2}.\end{aligned}$$

**Figure 2.** Maximum error of Example 1, when *M* = 50.

Table 2 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 <sup>×</sup> 102 s. Additionally, as illustrated in Example 1, the graphs of numerical solution and point-wise maximum error of the problem are illustrated in Figures 3 and 4, respectively.


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

From the above two Examples, 1 and 2, we see that the solutions of the problems exhibited an interior layer at *x* = <sup>1</sup> <sup>2</sup> 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.

**Figure 3.** Numerical solution of Example 2.

**Figure 4.** Maximum error of Example 2, when *M* = 50.

#### **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 onedimensional 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 secondorder nature, but, due to the presence of the first-order term Δ*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 Figures 1 and 3 depict their solutions. The layer occurred at the interior line *x* = *z* and boundary line *x* = 1, which is evident from these figures. Figures 2 and 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 Tables 1 and 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 abovedefined 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.

**Author Contributions:** Methodology, R.S., V.S.; software, R.S., V.S.; formal analysis, R.S., V.S.; investigation, R.S., V.S.; writing—original draft preparation, R.S., V.S., P.J.Y.W.; writing—review and editing, R.S., V.S., P.J.Y.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
