1 Introduction↩︎

Singularly perturbed problems commonly emerge across various fields of science and engineering, such as fluid dynamics, ocean circulation, heat transfer, optimal control, and financial mathematics. One of the main characteristics of these problems is the involvement of a small parameter \(\varepsilon > 0\) in the highest derivative, which strongly determines the nature of the solution. Solutions of such equations exhibits multiscale phenomena in general i.e. in certain thin region(s) of the domain the solution and/or its derivatives changes rapidly and outside of it behaves regularly and varies slowly. These regions, characterized by sharp gradients, are referred to as boundary layers or interior layers.

Several numerical methods are available in the literature to solve singularly perturbed differential equations, for details one can refer the books [1]–[3]. Because of the formation of boundary and/or interior layers, numerical methods face major computational difficulties. Finite difference and classical finite element methods (FDMs/FEMs) yield inaccurate solutions and produce non physical oscillations on uniform meshes. Therefore, one requires to develop \(\varepsilon-\)uniform methods to solve such problems. One can find various works on solving singular perturbation problems (SPPs) with boundary layers. Beckett and Mackenzie in [4], derive \(\varepsilon-\)uniform error estimates for two first-order upwind discretization of a SPP in 1D over a non-uniform mesh. Clavero and Gracia in [5] considered a 1D parabolic SPP of reaction-diffusion type. They implemented backward-Euler in time over a uniform mesh and a HODIE FDM in space over a generalized Shishkin mesh and proved first order convergence in time and almost fourth order convergence in space. In [6] Avijit and Natesan studied the convergence behavior of the streamline diffusion FEM (SDFEM) applied to time-dependent singularly perturbed convection-diffusion problems in 1D. In [7], Kaland and Roos implemented Galerkin FEM with \(\theta-\)scheme and discontinuous Galerkin (dG) for time discretization for a parabolic IBVP in 1D. They obtain the convergence of order almost one for both the schemes in time. Zarin and Roos in [8] used the \(h\)-version of non-symmetric dG-FEM with interior penalities (NIPG) over Shishkin-type mesh to solve convection-diffusion problems in 2D. Singh and Natesan [9], considered a parabolic IBVP and discretized the temporal variable using dG on a uniform mesh and spatial variable using NIPG on a Shishkin mesh. They stuidied the superconvergence properties of the proposed scheme. The literature offers limited insight into numerical approaches for solving singularly perturbed turning point problems (SPTPPs) as such problems are harder to solve than nonturning point problems.

Natesan et al. in [10] introduced an \(\varepsilon-\)uniform method for SPTPPs exhibiting boundary layers. Majumdar and Natesan in [11] presented a numerical strategy exhibiting uniform convergence for a class of convection-diffusion SPPs in 2D exhibiting interior layer due to non-smooth convection and source term. Lin\(\ss\) in [12] analyzed the convergence of an upwind method on Shishkin meshes for a class of semilinear SPPs with boundary turning point. In [13], Dunne et al. developed an \(\varepsilon-\)robust numerical method for parabolic singularly perturbed boundary turning point problems (SPBTPPs) by discretizing the temporal variable using the backward-Euler scheme over a uniform grid combined with spatial discretization via the standard upwind approach on a fitted mesh. They proved an almost first order convergence of the proposed method in the maximum norm. Chen et al. in [14] considered a 1D SPBTPP and implemented a variant of SDFEM. Stability and convergence estimates of the proposed scheme were proved on Shishkin-type mesh and Bakhvalov-type mesh. Some other works on solving SPBTPPs can be found in [15], [16] and references therein.

Recently, a non-conforming FEM known as the weak Galerkin FEM (WG-FEM) is becoming popular among the researchers of SPP community. Originally developed in [17], this approach was designed to handle second-order elliptic PDEs. The idea was to approximate the classical gradient operator with weak gradient in the sense of distributions. The adoption of discontinuous basis functions enables the flexible and efficient construction of higher-order elements on arbitrary partitions of the computational domain, thereby accommodating complex mesh configurations with ease. WG-FEM has been widely explored in addressing problems encountered in physics and engineering such as for Stokes equations [18], Maxwell’s equations [19], Brinkman equations [20], etc. Within the literature, WG-FEM have been extensively applied to SPPs, encompassing both elliptic and parabolic types. In [21], Zhu and Xie proposed a WG-FEM for one-dimensional SPPs, utilizing Shishkin mesh to resolve boundary layers. In [22], Zhang and Liu proposed a WG-FEM to address SPPs in 2D over a Shishkin mesh. Tao et al. [23] established a WG-FEM framework for solving 1D reaction-diffusion SPPs on a Shishkin mesh. Toprakseven in [24] analyzed the superconvergence properties of a modified WG-FEM on a Shishkin mesh. Authors in [25] proposed a WG-FEM approximation for parabolic fourth-order SPPs in 2D over a Shishkin mesh. High-dimensional parabolic equations often suffer from the curse of dimensionality, leading to significant computational challenges. To mitigate this issue researchers are employing operator-splitting techniques [26]–[28]. This strategy not only streamlines the solution process but also lowers computational costs and enhances flexibility when dealing with complex geometries. Despite these advantages, the integration of operator-splitting ADI methods within the WG-FEM framework for tackling high-dimensional SPPs has received limited attention in current research [29]–[31].

The contributions of our work are enlisted as follows:

1. We propose the WG-FEM for SPBTPPs, addressing a previously unstudied class of problems.

2. We develop a numerical scheme combining an implicit \(\theta\)-scheme for temporal discretization on a uniform mesh with WG-FEM spatial discretization on a Shishkin mesh.

3. Stability estimates for semi- and fully-discrete formulations are rigorously established.

4. Error estimates in the energy norm are derived, and \(\varepsilon\)-uniform convergence of the proposed scheme is proved.

5. We establish theoretical foundations for future operator splitting based WG-FEM schemes for SPBTPPs in higher dimensions, leveraging the 1D bounds derived in this work.

The organization of the paper is as follows: Section 2 describes the model problem, presents the solution decomposition, and establishes analytical bounds. In Section 3, we develop the Shishkin mesh framework for spatial discretization. Section 4 details the Weak Galerkin formulation, including both semi- and fully discrete schemes, along with their corresponding stability estimates. Section 5 derives the error equations for both schemes and provides a comprehensive error analysis. Numerical experiments validating the robustness of the proposed scheme are presented in Section 6. Finally, Section 7 concludes the article with a summary of key findings and future research directions. In the text, \(C\) denotes a constant that does not depend on \(\varepsilon\), \(N\) or \(\mathcal{M}\), and it may vary in different contexts.

2 Model Problem↩︎

The problem under consideration is the parabolic singularly perturbed boundary turning point problem (SPBTPP) given by: \[\label{spbtp95md95problem} \left\{\begin{array}{ll} u_{t} -\varepsilon u_{xx} -x^{q}\cdot b(x,t) u_{x} +c(x,t)u = {f}(x,t), \,\, q \geq 1, \, (x,t) \in \,\, \Omega \times (0,T],\\[6pt] u(x,0) = u_{0}(x), \,\, x \in [0,1], \\[6pt] u(0,t) = 0, u(1,t) = 0, \,\, t \in [0,T]. \end{array}\right.\tag{1}\] Here \(0 <\varepsilon \ll 1\) is a singular perturbation parameter and \(\Omega = (0,1)\). Considering that \(b(x,t), c(x,t)\) and \(f(x,t)\) are sufficiently smooth in their respective domains, the uniqueness of the solution to 1 is ensured. Also, we denote the convection coefficient as \(x^qb(x,t) = a(x,t)\). The solution of 1 exhibits a turning point behavior when the convection coefficient becomes zero at some location within the domain \(\Omega\). Specifically, if this coefficient vanishes at the domain boundary, the problem is termed a boundary turning point problem. If it vanishes at an interior point, it is referred to as an interior turning point problem. Our problem 1 contains a turning point at \(x=0\), with multiplicity \(q\). If \(q=1\), the problem corresponds to a simple boundary turning point case and multiple boundary turning point case if \(q > 1\). We also have the following assumptions \[\label{coeff95assumptions} \left\{\begin{array}{ll} b(x,t) \geq \mathfrak{b} > 0, \quad (x,t) \in [0,1] \times [0,T], \\[6pt] c(x,t) \geq \mathfrak{c} > 0, \quad (x,t) \in [0,1]\times [0,T]. \end{array}\right.\tag{2}\] The solution of 1 under the assumptions 2 , shows a parabolic boundary layer forming near the vicinity of \(x=0\) with a width of \(\mathcal{O}(\sqrt{\varepsilon})\).

2.1 Solution decomposition↩︎

For error analysis, it is essential to estimate the solution as well as its derivatives. These estimates are detailed in the following lemma.

In order to drive higher order accurate error estimates, we provide sharper bounds on the solution of 1 in the following lemma.

For more details into these lemmas one can refer to [16].

3 Mesh construction for Spatial Discretization↩︎

In this section, we will discuss the construction of a layer adapted mesh, specifically Shishkin mesh for spatial discretization.

3.1 Construction of Shishkin mesh↩︎

We will present here the construction of Shishkin mesh suitable for our SPBTPP 1 posed in [0,1]. We first choose the mesh transition parameter denoted by \(\tau\) as: \[\label{ch1195transition} \tau = \min\Big\{\dfrac{1}{2}, \frac{(k+1)}{\sqrt{\mathfrak{b}}}\sqrt{\varepsilon}\ln N\Big\},\tag{3}\] where \(k\) corresponds to the order of polynomials chosen for the approximation space. Using the transition parameter \(\tau\) we divide our computational domain \(\overline{\Omega}\) into \([0,\tau]\) and \([\tau,1]\) (See Figure 1).

Then each sub-domain is subsequently split into \(N/2\) equally sized sub-intervals. The mesh nodes are determined as \[\label{ch1195xmesh} x_{i} = \left\{\begin{array}{ll} \dfrac{2i\tau}{N}, \quad &i = 0,\ldots,\frac{N}{2},\\[6pt] \tau + \dfrac{2(i-\frac{N}{2})(1-\tau)}{N}, \quad &i = \frac{N}{2}+1,\ldots,N, \end{array}\right.\tag{4}\] and the step size as \[\label{ch1195xmesh95stepsize} h_{i} = \left\{\begin{array}{ll} \dfrac{2\tau}{N}, \quad &i = 0,\ldots,\frac{N}{2},\\[6pt] \dfrac{2(1-\tau)}{N}, \quad &i = \frac{N}{2}+1,\ldots,N. \end{array}\right.\tag{5}\] The step size \(h_i = x_{i+1} - x_{i}\) has the following order \[\begin{align} h_i &=& \mathcal{O}(\sqrt{\varepsilon}N^{-1}\ln N), \quad i \in \{0,\ldots, N/2\}, \\[2pt] h_i &=& \mathcal{O}(N^{-1}), \quad i \in \{N/2+1,\ldots, N\}. \end{align}\]

4 WG-FEM for Spatial Discretization↩︎

Let \(M_{x_i} = [x_{i-1}, x_i]\) represent a mesh element. The domain \(\Omega\) is divided into elements \(M_{x_i}\), forming the partition \(\mathcal{T}_N = \{ M_{x_i} \mid i=1,\ldots,N \}\). Let \(h = \max\{h_i\}\) denote the mesh size, where \(h_i\) is the length of the local element \(M_{x_i}\) for \(i = 1, 2, \ldots, N\).

A weak function defined on the element \(M_{x_i}\) is given by \(u_h = \{ u_0, u_b \}\), where \(u_0 \in L^2(M_{x_i})\) and \(u_b \in L^\infty(\partial M_{x_i})\). The boundary of the element is denoted by \(\partial M_{x_i} = \{ x_{i-1}, x_i \}\). In this formulation, \(u_0\) represents the function’s value inside the element, while \(u_b\) corresponds to its values at the endpoints.

Specifically, for each element \(M_{x_i} \in \mathcal{T}_N\), the weak function \(u_h\) is defined as: \[u_h = \begin{cases} u_0, & \text{in } M_{x_i}, \\ u_b, & \text{on } \partial M_{x_i}. \end{cases}\] We now introduce the discrete weak function space associated with each element \(M_{x_i}\) as follows: \[\mathscr{W}(M_{x_i}, k) = \left\{ u_h = \{ u_0, u_b \} \, \middle| \, u_0|_{M_{x_i}} \in \mathbb{P}_k(M_{x_i}), \, u_b|_{\partial M_{x_i}} \in \mathbb{P}_0(\partial M_{x_i}) \right\},\] where \(\mathbb{P}_k(M_{x_i})\) denotes the space of polynomials of degree at most \(k \in \mathbb{Z}^+\) defined on \(M_{x_i}\), and \(\mathbb{P}_0(\partial M_{x_i})\) represents the space of functions that are piecewise constant on the boundary \(\partial M_{x_i}\).

We define the global weak Galerkin space \(\mathscr{S}_h\) by assembling the local spaces \(\mathscr{W}(M_{x_i}, k)\) and enforcing continuity at the interior mesh nodes. Specifically, we set \[\mathscr{S}_h = \mathscr{S}_h(\Omega) = \left\{ u_h = \{ u_0, u_b \} \, : \, \{ u_0, u_b \}|_{M_{x_i}} \in \mathscr{W}(M_{x_i}, k), \; \forall M_{x_i} \in \mathcal{T}_N \right\}.\]

We further define a subspace \(\mathscr{S}_h^0 \subset \mathscr{S}_h\) consisting of weak functions whose boundary components vanish at the endpoints of the domain. That is, \[\label{ch1195zero} \mathscr{S}_h^0 = \left\{ u_h \in \mathscr{S}_h \, : \, u_b(0) = u_b(1) = 0 \right\}.\tag{6}\]

Equation 6 shows that boundary conditions are strongly imposed in the proposed WG-FEM.

4.1 Weak Derivative↩︎

For a weak function \(u_h \in \mathscr{W}(M_{x_i}, k)\), the weak derivative \(d_{w, M_{x_i}} u_h \in \mathbb{P}_{k-1}(M_{x_i})\) is the unique polynomial that satisfies \[\label{weakderivative} (d_{w, M_{x_i}} u_h, \phi)_{M_{x_i}} = - (u_0, \phi')_{M_{x_i}} + \langle u_b, \phi n \rangle_{\partial M_{x_i}}, \quad \forall \phi \in \mathbb{P}_{k-1}(M_{x_i}),\tag{7}\] where the inner product and boundary pairing are defined by \[(\kappa_1, \kappa_2)_I := \int_I \kappa_1(x) \kappa_2(x) \, dx, \quad \text{and} \quad \langle \kappa_1, \kappa_2 n \rangle_{\partial I} := \kappa_1(x_i)\kappa_2(x_i) - \kappa_1(x_{i-1})\kappa_2(x_{i-1}),\] for \(I = (x_{i-1}, x_i)\).

To approximate the convection term, we define the weak convection derivative as follows: For any weak function \(u_h \in \mathscr{W}(M_{x_i}, k)\), the weak convection derivative \(d_{w, M_{x_i}}^{a} u_h \in \mathbb{P}_k(M_{x_i})\) is the unique polynomial satisfying \[\label{convectionderivative} (d_{w, M_{x_i}}^{a} u_h, \phi)_{M_{x_i}} = - (u_0, (a\phi)')_{M_{x_i}} + \langle u_b, a\phi n \rangle_{\partial M_{x_i}}, \quad \forall \phi \in \mathbb{P}_k(M_{x_i}).\tag{8}\]

We compute the weak derivatives \(d_w\) and \(d_w^{a}\) elementwise within the global WG finite element space \(\mathscr{S}_h\) on each interval \(M_{x_i}\) for \(i = 1, \dots, N\). Specifically, \[(d_w u_h)|_{M_{x_i}} = d_{w, M_{x_i}}(u_h|_{M_{x_i}}), \quad (d_w^{a} u_h)|_{M_{x_i}} = d_{w, M_{x_i}}^{a}(u_h|_{M_{x_i}}), \quad \forall u_h \in \mathscr{S}_h.\]

4.2 Weak Formulation and Well-Posedness↩︎

In this subsection, we present both the semi-discrete and fully-discrete WG-FEM formulations for problem 1 , employing implicit \(\theta-\)scheme for time discretization.

Let \(u_h \in \mathscr{S}_{h}^{0}\) and \(v_h \in \mathscr{S}_{h}^{0}\) be given. The semi-discrete formulation of 1 is: find \(u_h \in \mathscr{S}_{h}^{0}\) such that \[\label{ch1195semi95discrete} ((u_{0})_t, v_0) + A(u_h, v_h) = (f, v_0),\tag{9}\] with the bilinear form \(A(\cdot,\cdot)\) specified by \[\begin{align} \label{A95bilinear} A(u_h, v_h) &=& \sum_{M_{x_i} \in \mathcal{T}_{N}} \varepsilon (d_{w} u_h, d_{w} v_h)_{M_{x_i}} - \sum_{M_{x_i} \in \mathcal{T}_{N}} (d_{w}^{a} u_h, v_0)_{M_{x_i}} + \sum_{M_{x_i} \in \mathcal{T}_{N}} (c u_0, v_0)_{M_{x_i}} \nonumber \\[2pt] && + s_d(u_h, v_h) + s_c(u_h, v_h), \end{align}\tag{10}\] and the stabilizer terms corresponding to diffusion and convection are given by \[\begin{align} s_d(u_h, v_h) &=& \sum_{i=1}^{N} \langle \vartheta^i (u_0 - u_b), v_0 - v_b \rangle_{\partial M_{x_i}}, \\ s_c(u_h, v_h) &=& \sum_{i=1}^{N} \langle a n (u_0 - u_b), v_0 - v_b \rangle_{\partial_+ M_{x_i}}, \end{align}\] with \[\partial_+ M_{x_i} = \{x \in \partial M_{x_i} : a(x) n_{M_{x_i}}(x) \ge 0\},\] and \(\vartheta^i\) defined as a penalty parameter: \[\begin{align} \label{vartheta} \vartheta^i = \begin{cases} N(\ln N)^{-1}, & i \in \{1, 2, \dots, N/2\}, \\ 1, & i \in \{N/2+1, \ldots, N\}. \end{cases} \end{align}\tag{11}\] Now, we define the energy norm as for any \(v_h \in \mathscr{S}_h\) \[\label{energynormx} \||{v}_h\||^2 \,\, := \, \left(\vert {v}_h \vert_{w}\right)^2+ \Vert \sqrt{\gamma} {v}_0\Vert^2 +\left(\vert {v}_h\vert_{C}\right)^2,\tag{12}\] where \[\begin{align} \left(\vert {v}_h \vert_{w}\right)^2 &:=& \varepsilon\Vert d_w {v}_h\Vert^2 + s_d({v}_h, {v}_h), \\[4pt] c(x) &+& \frac{a'(x)}{2} \geq \gamma > 0, \\[4pt] \left(\vert {v}_h \vert_{C}\right)^2 &:=& \sum _{i=1}^{ N}d_i^*|\sqrt{a}({v}_0 - {v}_b)|^2(x_i^-), \end{align}\] with \[\begin{align} d_i^*=\left\{ \begin{array}{ll} \frac{3}{2},&{}\quad i = N,\\[6pt] 1,&{}\quad i=1,\ldots , N-1. \end{array}\right. \end{align}\] Further, we measure the errors in the discrete \(H^{1}-\)energy norm defined as \[\label{energynormxH1} \||{v}_h\||_{\mathcal{E}}^2 \,\, := \, \left(\vert {v}_h\vert_{1}\right)^2+ \Vert \sqrt{\gamma} {v}_0\Vert^2 + \left(\vert {v}_h\vert_{C}\right)^2,\tag{13}\] with \[\left(\vert {v}_h\vert_{1}\right)^2 \,\, := \, \varepsilon\Vert D {v}_0\Vert^2 + s_d({v}_h, {v}_h),\] where \(D{v}\) refers to the classical derivative of \(v\) taken with respect to the space variable.

Using the stability estimate ?? , one can prove that the semi-discrete formulation 9 admits atmost one solution and hence well-posed.

4.2.1 Fully discrete formulation↩︎

We will introduce an evenly spaced discretization for time domain \((0,T]\). Let \(\mathcal{M}\) be the number of partitions over the time interval. Choosing the step size \(\zeta = T/\mathcal{M}\), we can define the mesh as \[\label{ch1195time95mesh} \mathscr{T}^{\mathcal{M}} = \{t_j: t_j = j\zeta, \,\, j = 0,1,\ldots, \mathcal{M} \}.\tag{14}\] Discretizing the semi-discrete formulation 9 in time via the implicit \(\theta\)-scheme yields the fully-discrete problem:

Find \(u_{h}^{n} \in \mathscr{S}_{h}^{0}\) such that \[\label{fully95discrete95form} \Bigg\{\begin{array}{ll} (\overline{\partial} u_h^{n}, v_0) + A(\theta u_{h}^{n} + (1-\theta)u_{h}^{n-1}, v_{h}) = (\theta f(t_n) + (1-\theta)f(t_{n-1}), v_{0}), \quad \forall v_h \in \mathscr{S}_{h}^{0},\\[6pt] u_{h}^{0} = \mathscr{I}u_{0}, \end{array}.\tag{15}\] where \(\theta \in [\frac{1}{2}, 1]\) and \(\overline{\partial} u_h^{n} := \dfrac{u_{h}^{n} - u_{h}^{n-1}}{\zeta}\).

5 Error Estimates for Semi- and Fully-Discrete WG-FEM↩︎

In this section, we shall present the error estimates of the proposed semi- and fully-discrete schemes. We first introduce interpolation operator as follows: Let the collection of \(k+1\) nodal functions on \({M}_{x_i} \in \mathcal{T}_{N}\) be given by \[N_{0}(s) = s(x_{i-1}),\quad N_{k}(s) = s(x_{i}),\] \[N_{l}(s) = h_{i}^{-l}\int_{{M}_{x_i}}^{} (x-x_{i-1})^{k-1}\,s(x)\,dx,\,\, l = 1,\ldots,k-1.\] Using these nodal functions, we proceed to construct a local interpolation operator \(\mathscr{I} : H^{1}(M_{x_i}) \to \mathbb{P}_{k}(M_{x_i})\) by \[\label{ch1195inter95op} N_{l}(\mathscr{I}s - s) = 0,\,\, l = 0,\ldots,k.\tag{16}\] Since \(\mathscr{I}s\vert_{{M}_{x_i}}\) is continuous, \(\{\mathscr{I}s\vert_{M_{x_i}}, \mathscr{I}s\vert_{\partial M_{x_i}}\}\) is in our weak finite element space \(\mathscr{S}_{h}\).

Using Lemma 2 in conjunction with Lemma 6, we obtain the following bounds for the interpolation error.

In the following, we establish error equations to analyze the consistency error arising from WG-FEM.

Now, we give an estimate for the semi-discrete scheme 9 using Theorem 1.

We now present the following theorem, which provides an error estimate for the fully discrete scheme given in 15 .

6 Numerical Experiments↩︎

To demonstrate the effectiveness and accuracy of the proposed method, we present numerical results that both support and validate the theoretical predictions. We discretize in time via the backward Euler (BE, \(\theta = 1\)) and Crank-Nicolson (CN, \(\theta = 1/2\)) methods. Let \(e_h^n = u^n - u_h^n\) be the error. To compute the convergence order, we apply the formula \[Order = \dfrac{\ln(e_{h}^{n}/e_{2h}^{n})}{\ln(2\ln(N)/\ln(2N))}.\]

To obtain the exact solution \[u(x,t) = (1 - e^{-t})\Big(1 - x + xe^{-1/\sqrt{\varepsilon}} - e^{-x/\sqrt{\varepsilon}}\Big),\] the source term \(f(x,t)\) is appropriately chosen.

We obtain the energy norm error and the corresponding order of convergence of our proposed scheme for various values of \(\varepsilon, N\) and \(k\), which are displayed in Tables 1 and 2. For spatial error to dominate we fixed \(\zeta = 1/5000\). These tables illustrate that the energy norm convergence rate is of order \(\mathcal{O}((N^{-1}\ln N)^{k})\), which closely aligns with the Theorem 4. Figures 2 and 4 each present two subfigures: a surface plot of the numerical solution and a comparison of the exact and numerical solution at the final time \(T=1\), corresponding to \(\theta = 1\) and \(\theta = 1/2\), respectively. To validate the theoretical convergence rates, we plot the energy norm errors with respect to mesh size \(N\) on a logarithmic scale, for varying perturbation parameters \(\varepsilon\) and polynomial degrees \(k = 1, 2\). Figure 3 depicts these results for \(\theta = 1\), while Figure 5 illustrates the corresponding behavior for \(\theta = 1/2\). Both figures highlight the alignment of numerical errors with the predicted \(\mathcal{O}((N^{-1} \ln N)^k)\) rates, demonstrating the schemes’ robustness across parameter regimes.

To highlight the differences in convergence behavior, we fixed a relatively coarse temporal step size (\(\zeta = 1/500\)) and evaluated the errors for \(\theta = 1\) and \(\theta = 1/2\). As presented in Table 3, the scheme with \(\theta = 1\) shows deteriorating convergence rates with increasing \(N\), due to the dominant temporal discretization error. In contrast, the scheme with \(\theta = 1/2\) retains consistent convergence, owing to its higher temporal accuracy.

To demonstrate how the parameter \(q\) affects the solution of the degenerate PDE, we evaluate the errors in the energy norm along with the associated convergence rates for different values of \(q\), keeping \(\varepsilon = 10^{-8}\). The results, presented in Table 4, indicate that the error bound is independent of \(q\). In Table 5, the maximum norm errors are compared with the values reported in [33].

7 Conclusion and future direction↩︎

In conclusion, we have developed and rigorously analyzed a WG-FEM for SPBTPPs, employing an implicit \(\theta\)-scheme for temporal discretization. The proposed scheme exhibits stability as well as uniform convergence with respect to the energy norm, as validated by comprehensive numerical experiments. As part of ongoing work, we are focusing on formulating and analyzing ADI-type operator splitting schemes within the WG-FEM framework for multidimensional SPBTPPs.

Declarations↩︎

Funding↩︎

The first author gratefully acknowledges the financial support provided by the Indian Institute of Technology Guwahati.

Data Availability↩︎

Data sharing is not applicable to this paper because no datasets were created or examined during the current study.

Conflict of interest↩︎

The authors declare that they have no conflict of interest.

References↩︎