A Local Radial Basis Function Method for Numerical Approximation of Multidimensional Multi-Term Time-Fractional Mixed Wave-Diffusion and Subdiffusion Equation Arising in Fluid Mechanics

Abstract

This article develops a simple hybrid localized mesh-free method (LMM) for the numerical modeling of new mixed subdiffusion and wave-diffusion equation with multi-term time-fractional derivatives. Unlike conventional multi-term fractional wave-diffusion or subdiffusion equations, this equation features a unique time–space coupled derivative while simultaneously incorporating both wave-diffusion and subdiffusion terms. Our proposed method follows three basic steps: (i) The given equation is transformed into a time-independent form using the Laplace transform (LT); (ii) the LMM is then used to solve the transformed equation in the LT domain; (iii) finally, the time domain solution is obtained by inverting the LT. We use the improved Talbot method and the Stehfest method to invert the LT. The LMM is used to circumvent the shape parameter sensitivity and ill-conditioning of interpolation matrices that commonly arise in global mesh-free methods. Traditional time-stepping methods achieve accuracy only with very small time steps, significantly increasing the computational time. To overcome these shortcomings, the LT is used to provide a more powerful alternative by removing the need for fine temporal discretization. Additionally, the Ulam–Hyers stability of the considered model is analyzed. Four numerical examples are presented to illustrate the effectiveness and practical applicability of the method.

1. Introduction

Time-fractional partial differential equations (TFPDEs) have garnered significant attention in the scientific community due to their ability to model nonlocal and memory-dependent phenomena [1]. TFPDEs find numerous applications across various fields, including physics [2,3], biology [4], chemistry [5], hydrology [6], finance [7], and more. One significant class of TFPDEs is that of time-fractional multi-term diffusion equations (TFDEs), which can be further classified into three different categories according to their fractional order: (i) time-fractional subdiffusion equations (TFSDEs), where the order of the Caputo derivative falls within the interval $(0,1)$; (ii) time-fractional wave-diffusion equations (TFWDE), where the order of the Caputo derivative is in the interval $(1,2)$; and (iii) time-fractional diffusion equations (TFDE), where the order of the Caputo derivative lies within the interval $(0,2).$ However, using TFSDEs or TFWDEs independently may not adequately capture features of some physical processes. To overcome these limitations, a mathematical model that combines TFSDEs and TFWDEs, known as a mixed TFSDE and TFWDE, can be used. These mixed equations can provide enhanced accuracy and adaptability, especially for multi-term time-fractional derivatives [8].

Numerous researchers have studied the analytical solutions to these equations. For instance, Daftardar et al. [9] solved a multi-term TFWDE using the separation of variables method. Baranwal et al. [10] utilized the generalized n-dimensional differential transform method to solve heat and wave-like equations of fractional order. Luchko [11] studied the uniqueness and a priori estimates of the solution for a multi-term TFSDE using a plausible maximum concept. Shah et al. [12] derived a semi-analytical solution to fractional-order partial Benney equation. Jiang et al. [13] discussed three potential boundary conditions for a multi-term TFWDE. Shen et al. [14] used the properties of the multivariate Mittag–Leffler (ML) function and the separation of variables approach to solve the multi-term TFWDE and TFDE.

While multi-term TFDEs have some analytical solutions, these solutions involve ML functions, which are exceedingly complicated and challenging to evaluate. As a result, numerous numerical techniques have been developed to solve multi-term TFDEs and TFWDEs. For instance, Hemati et al. [15] developed a reproducing kernel method to investigate the solution of multi-term TFDEs. Rashidinia and Mohmedi [16] obtained approximate TFDE and TFWDE solutions using the spectral technique. Liu et al [17] studied the approximate solution of a multi-term TWDE using fractional predictor–corrector methods. Feng et al. [18] solved mixed TFSDE and TFWDE equations using the finite element/finite difference methods. Bhrawy et al. [19] applied the spectral tau method with Jacobi operational matrix to solve TFWDEs. Li et al. [20] used the local meshless approach based on a multiquadric (MQ) kernel to investigate the solution of multi-term TFWDEs. Bhardwaj and Kumar [21] studied the solution of a mixed TFDE and TFWDE using the LMM for spatial derivatives and the finite difference method for the time-fractional derivative. Shen et al. [22] developed two numerical schemes for a mixed TFDE and TFWDE equation. Zhang et al. [23] studied the two-grid method for 2D nonlinear multi-term mixed TFSDE and TFWDE equations. Other valuable numerical techniques for solving such TFPDEs are discussed in [24,25] and the references therein.

In this article, our objective is to apply the LMM based on the multiquadric radial basis function (MQRBF) and Gaussian radial basis function (GSRBF) coupled with LT to solve multi-term TFSDEs and TFWDEs. The LMM utilizes a set of scattered nodes instead of conventional domain meshing, resulting in reduced computational cost. Additionally, the solution is computed on the local support domains; these local approximations result in small sparse systems that need to be inverted for each node to obtain the local weights for the differential operators. After the local weights are computed, the global matrix can be assembled. Due to its sparse nature, the global matrix is better conditioned and easier to solve than global meshless methods. The same problem was solved in [26] using the Chebyshev spectral collocation method coupled with LT. Spectral methods are highly efficient for problems defined on regular domains due to their fast convergence properties; however, they have limitations on irregular domains, where it becomes hard to generate global basis functions. On the other hand, the LMM requires scattered nodes, eliminating the need for structured meshing. This flexibility makes it more versatile and advantageous, especially for complex geometries where spectral methods may struggle.

The LT is used to avoid the finite difference method (FDM) for time discretization, as the FDM does not always yield a stable solution. With the FDM, stability is achieved only if the errors decay or remain constant throughout the computation. Furthermore, the optimal numerical solution with the FDM often requires very small time steps, which increases the computational time and reduces efficiency, especially for fractional-order problems. However, one of the main challenges with the LT is computation of its inverse. The literature offers various numerical approaches for inverting the LT. In this work, we employ the improved Talbot method [27] and the Stehfest method [28] for LT inversion. We consider the general version of multi-term TFDEs introduced by Feng et al. [18], known as the multi-term mixed TFSDE and TFWDE. The main contributions of this work are as follows:

• A new multi-term fractional-order mixed subdiffusion and diffusion-wave equation is considered. This equation features a unique time–space coupled derivative, and incorporates both wave-diffusion and subdiffusion terms.
• A hybrid method combining the LT method and the LMM is employed to approximate the solution of the considered equation. This method employs the LT and its inverse to replacing finite difference time-stepping, thereby avoiding the costly convolution integral calculation typically associated with time-fractional derivative approximations. The LMM is employed for spatial discretization, effectively circumventing issues related to shape parameter sensitivity and ill-conditioning of interpolation matrices encountered in global mesh-free methods.
• Five examples in one, two, and three dimensions are solved using the proposed scheme. The numerical results demonstrate that the method is robust, stable, and effective. Moreover, the method is easy to use and can be applied to other fractional-order problems, such as the generalized Oldroyd-B fluid model [29], time-fractional telegraph equation [30], and others.

The rest of this article is organized as follows: in Section 2, some basic definitions are provided; Section 3 covers the Ulam–Hyers stability of the considered model; Section 4.1 discusses the LT method; Section 4.2 covers the LMM method for spatial discretization; Section 4.3 addresses numerical inversion of the LT method; Section 4.4 presents the error analysis; the Stability of the proposed numerical scheme is examined in Section 5; Section 6 presents the numerical results; finally, our conclusions is provided in Section 7.

2. Basic Definitions

Below, we present some definitions from fractional calculus relevant to our study, along with important lemmas and theorems.

Definition 1.

The Gamma function is defined by

$$\mathsf{\Gamma }\left(\varsigma \right)={\int }_0^{\infty }{\vartheta }^{\varsigma -1}exp(-\vartheta )d\vartheta ,\varsigma >0$$

and $\mathsf{\Gamma }(1+\varsigma )=\varsigma \mathsf{\Gamma }\left(\varsigma \right).$

Definition 2

([31]). The integral of fractional order $\alpha \ge 0$ in the Riemann–Liouville sense of a function $\mathcal{W}$ is defined as

$$J_{\theta }^{\alpha }\mathcal{W}\left(\theta \right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^{\theta }{(\theta -\tau )}^{\alpha -1}\mathcal{W}\left(\tau \right)d\tau ,\alpha >0,\theta >0.$$

The operator $J_{\theta }^{\alpha }$ possesses the following properties:

$$\left(i\right)J_{\theta }^{\alpha }{J}_{\theta }^{\varpi }\mathcal{W}\left(\theta \right)=J_{\theta }^{\alpha +\varpi }\mathcal{W}\left(\theta \right),$$

(1)

$$\left(ii\right)J_{\theta }^{\alpha }{J}_{\theta }^{\varpi }\mathcal{W}\left(\theta \right)=J_{\theta }^{\varpi }{J}_{\theta }^{\alpha }\mathcal{W}\left(\theta \right),$$

(2)

$$\left(iii\right)J_{\theta }^{\alpha }{\theta }^{\beta }={\displaystyle \frac{\mathsf{\Gamma }(\beta +1)}{\mathsf{\Gamma }(\alpha +\beta +1)}}{\theta }^{\alpha +\beta }.$$

(3)

Definition 3.

The Caputo time-fractional derivative $D_{\tau }^{\mu }u(\overline{\mathit{\varrho }},\tau )$ of fractional orders μ is defined as

$${}_0{}^{c}D_{\tau }^{\mu }u(\overline{\mathit{\varrho }},\tau )={\displaystyle \frac{{\partial }^{\mu }u(\overline{\mathit{\varrho }},\tau )}{\partial {\tau }^{\mu }}}\left\{\begin{array}{cc}{\displaystyle \frac{1}{\mathsf{\Gamma }(m-\mu )}}{\int }_0^{\tau }{\displaystyle \frac{{\partial }^{m}u(\overline{\mathit{\varrho }},\upsilon )}{\partial {\upsilon }^m}}{(\tau -\upsilon )}^{m-\mu -1}d\upsilon ;\hfill & m-1<\mu <m,\hfill \\ {\displaystyle \frac{{\partial }^{m}u(\overline{\mathit{\varrho }},\tau )}{\partial {\tau }^m}},\hfill & \mu =m\in \mathbb{N}.\hfill \end{array}\right.$$

3. Ulam–Hyers Stability

In this section, we establish the Ulam–Hyers stability for the multi-term TFSDE and TFWDE. Let $C(\mathsf{\Lambda },\mathbb{R})$ denote the Banach space of continuous functions from $\mathsf{\Lambda }=\mathsf{\Omega }\times I$ into $\mathbb{R},$ equipped with the maximum norm ${\parallel .\parallel }_{\infty }$ defined as ${\parallel u\parallel }_{\infty }=sup\left\{\right|u(\overline{\mathit{\varrho }},\tau );(\overline{\mathit{\varrho }},\tau )\in \mathsf{\Lambda }\left|\right\}$ [32]. We focus on the following multidimensional multi-term time-fractional mixed subdiffusion and diffusion-wave equation with Dirichlet boundary conditions:

$$\left\{\begin{array}{c}\sum _{i=1}^k{c}_{1,i}{D}_{\tau }^{{\nu }_i}u(\overline{\mathit{\varrho }},\tau )+c_2{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+\sum _{j=1}^l{c}_{3,j}{D}_{\tau }^{{\mu }_j}u(\overline{\mathit{\varrho }},\tau )+c_{4}u(\overline{\mathit{\varrho }},\tau )\hfill \\ =c_5{\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},\tau )+c_6{D}_{\tau }^{\alpha }\left({\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},\tau )\right)+f(\overline{\mathit{\varrho }},\tau ),(\overline{\mathit{\varrho }},\tau )\in \mathsf{\Omega }\times [0,T]\hfill \\ u(\overline{\mathit{\varrho }},0)={\varphi }_1\left(\overline{\mathit{\varrho }}\right)u_{\tau }(\overline{\mathit{\varrho }},0)={\varphi }_2\left(\overline{\mathit{\varrho }}\right)\overline{\mathit{\varrho }}\in \overline{\mathsf{\Omega }},\hfill \\ {\mathcal{L}}_{\mathit{b}}u(\overline{\mathit{\varrho }},\tau )={\psi }_1(\overline{\mathit{\varrho }},\tau )\overline{\mathit{\varrho }}\in \partial \mathsf{\Omega },\tau \in [0,T],\hfill \end{array}\right.$$

(4)

where $c_{1,i}\ge 0$ $c_{3,j}\ge 0,$ $c_s,s=2,4,5,6,$ wherein none of the coefficients are equal to zero simultaneously; moreover, $D_{\tau }^{{\nu }_i},$ $D_{\tau }^{{\mu }_j}$, and $D_{\tau }^{\alpha }$ are Caputo’s time fractional derivatives of orders $\nu$, $\mu$, and $\alpha$ such that $1<{\nu }_1<{\nu }_2<....<{\nu }_k<2,$ $0<{\mu }_1<{\mu }_2<...<{\mu }_l<1,0<\alpha <1$, $u(\overline{\mathit{\varrho }},\tau )\in C^2(\mathsf{\Omega }\times [0,T])$ is an unknown function to be determined,  ${\mathcal{L}}_{\mathit{g}}=\mathsf{\Delta }$ is the Laplace operator, ${\mathcal{L}}_{\mathit{b}}$ is the boundary operator, and $\mathsf{\Omega }$ is the domain with boundary $\partial \mathsf{\Omega }$. The forcing term $f(\overline{\mathit{\varrho }},\tau )$ is assumed to possess sufficient smoothness, and ${\varphi }_1\left(\overline{\mathit{\varrho }}\right)$, ${\varphi }_2\left(\overline{\mathit{\varrho }}\right)$, and ${\psi }_1(\overline{\mathit{\varrho }},\tau )$ are continuous functions.

Lemma 1

([33]). A function $u\in C(\mathsf{\Lambda },\mathbb{R})$ is considered a solution of the fractional integral equation

$$\begin{array}{cc}\hfill u(\overline{\mathit{\varrho }},\tau )=& {\varphi }_1+\tau {\varphi }_2+J_{\tau }^{{\nu }_i}(\left({\delta }_4\nabla u(\overline{\mathit{\varrho }},\tau )+{\delta }_5{D}_{\tau }^{\alpha }(\nabla u(\overline{\mathit{\varrho }},\tau ))+{\delta }_{6}f(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & -\left({\delta }_1{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+{\delta }_2{D}_{\tau }^{{\mu }_j}\left(u(\overline{\mathit{\varrho }},\tau )\right)+{\delta }_{3}u(\overline{\mathit{\varrho }},\tau )\right))\hfill \end{array}$$

(5)

if and only if u is the solution of the problem in (4), where $f(\overline{\mathit{\varrho }},\tau ):C\left(\mathsf{\Lambda }\right)\to C\left(\mathsf{\Lambda }\right)$ and $u(\overline{\mathit{\varrho }},\tau ):C(\mathsf{\Lambda },\mathbb{R})\to C(\mathsf{\Lambda },\mathbb{R})$ are continuous and

$${\delta }_1={\displaystyle \frac{c_2}{{\sum }_{i=1}^k{c}_{1,i}}},{\delta }_2={\displaystyle \frac{{\sum }_{j=1}^l{c}_{3,j}}{{\sum }_{i=1}^k{c}_{1,i}}},{\delta }_3={\displaystyle \frac{c_4}{{\sum }_{i=1}^k{c}_{1,i}}},$$

$${\delta }_4={\displaystyle \frac{c_5}{{\sum }_{i=1}^k{c}_{1,i}}},{\delta }_5={\displaystyle \frac{c_6}{{\sum }_{i=1}^k{c}_{1,i}}},{\delta }_6={\displaystyle \frac{1}{{\sum }_{i=1}^k{c}_{1,i}}},{\delta }_7={\delta }_7.$$

We consider the following inequality, which will be used to introduce Ulam–Hyers stability:

$$|D_{\tau }^{{\nu }_i}\overline{u}(\overline{\mathit{\varrho }},\tau )-F\left(\overline{u}(\overline{\mathit{\varrho }},\tau ),{\displaystyle \frac{\partial \overline{u}(\overline{\mathit{\varrho }},\tau )}{\partial \tau }},D_{\tau }^{{\mu }_j}\overline{u}(\overline{\mathit{\varrho }},\tau ),\mathsf{\Delta }\overline{u}(\overline{\mathit{\varrho }},\tau ),D_{\tau }^{\alpha }\left(\mathsf{\Delta }\overline{u}(\overline{\mathit{\varrho }},\tau )\right),f(\overline{\mathit{\varrho }},\tau )\right)|<ϵ.$$

(6)

Definition 4

([34]). A solution of the problem in (4) admits Ulam–Hyers stability if there exist a positive real number ${\kappa }_1$ such that for every solution $\overline{u}$ of the inequality in (6), there exist an exact solution u such that

$$\parallel u-\overline{u}{\parallel }_{\infty }<{\kappa }_1ϵ.$$

Before proceeding to prove the main results, we introduce the following hypothesis.

For any $(\overline{\mathit{\varrho }},\tau )\in C\left(\mathsf{\Lambda }\right),$ there exist positive real constants $m_1,m_2,m_3,m_4$ such that:

$\left(A_1\right)$

$$|\nabla u(\overline{\mathit{\varrho }},\tau )-\nabla \overline{u}(\overline{\mathit{\varrho }},\tau )|\le m_1|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|,$$

$\left(A_2\right)$

$$|D_{\tau }^{\alpha }(\nabla u(\overline{\mathit{\varrho }},\tau ))-D_{\tau }^{\alpha }(\nabla \overline{u}(\overline{\mathit{\varrho }},\tau ))|\le m_2|\nabla u(\overline{\mathit{\varrho }},\tau )-\nabla \overline{u}(\overline{\mathit{\varrho }},\tau )|\le m_1{m}_2|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|,$$

$\left(A_3\right)$

$$|D_{\tau }^{{\mu }_j}\left(u(\overline{\mathit{\varrho }},\tau )\right)-D_{\tau }^{{\mu }_j}\left(\overline{u}(\overline{\mathit{\varrho }},\tau )\right)|\le m_3|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|,$$

$\left(A_4\right)$

$$|{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}-{\displaystyle \frac{\partial \overline{u}(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}|\le m_4|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|.$$

Theorem 1.

If assumptions ($A_1$–$A_4$) hold, then the problem in (4) is Ulam–Hyers stable.

Proof. 

Let the exact solution to the problem in (4) be expressed as

$$\begin{array}{cc}\hfill u(\overline{\mathit{\varrho }},\tau )=& {\varphi }_1+\tau {\varphi }_2+J_{\tau }^{{\nu }_i}(\left({\delta }_4\nabla u(\overline{\mathit{\varrho }},\tau )+{\delta }_5{D}_{\tau }^{\alpha }(\nabla u(\overline{\mathit{\varrho }},\tau ))+{\delta }_{6}f(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & -\left({\delta }_1{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+{\delta }_2{D}_{\tau }^{{\mu }_j}\left(u(\overline{\mathit{\varrho }},\tau )\right)+{\delta }_{3}u(\overline{\mathit{\varrho }},\tau )\right)),\hfill \end{array}$$

(7)

and let $\overline{u}(\overline{\mathit{\varrho }},\tau )$ represent the approximate solution to the problem in (4), expressed as

$$\begin{array}{cc}\hfill \overline{u}(\overline{\mathit{\varrho }},\tau )=& {\varphi }_1+\tau {\varphi }_2+J_{\tau }^{{\nu }_i}(\left({\delta }_4\nabla \overline{u}(\overline{\mathit{\varrho }},\tau )+{\delta }_5{D}_{\tau }^{\alpha }(\nabla \overline{u}(\overline{\mathit{\varrho }},\tau ))+{\delta }_{6}f(\overline{\mathit{\varrho }},\tau )+{\mathsf{\delta }}_7\mathsf{\Phi }(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & -\left({\delta }_1{\displaystyle \frac{\partial \overline{u}(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+{\delta }_2{D}_{\tau }^{{\mu }_j}\left(\overline{u}(\overline{\mathit{\varrho }},\tau )\right)+{\delta }_3\overline{u}(\overline{\mathit{\varrho }},\tau )\right)).\hfill \end{array}$$

(8)

From Equations (7) and (8), we have

$$\begin{array}{cc}\hfill |u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|=& \left|J_{\tau }^{{\nu }_i}\right(\left({\delta }_4\nabla u(\overline{\mathit{\varrho }},\tau )+{\delta }_5{D}_{\tau }^{\alpha }(\nabla u(\overline{\mathit{\varrho }},\tau ))+{\delta }_{6}f(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & -\left({\delta }_1{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+{\delta }_2{D}_{\tau }^{{\mu }_j}\left(u(\overline{\mathit{\varrho }},\tau )\right)+{\delta }_{3}u(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & -\left({\delta }_4\nabla \overline{u}(\overline{\mathit{\varrho }},\tau )+{\delta }_5{D}_{\tau }^{\alpha }(\nabla \overline{u}(\overline{\mathit{\varrho }},\tau ))+{\delta }_{6}f(\overline{\mathit{\varrho }},\tau )+{\delta }_7\mathsf{\Phi }(\overline{\mathit{\varrho }},\tau )\right)\hfill \\ \hfill & +\left({\delta }_1{\displaystyle \frac{\partial \overline{u}(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}+{\delta }_2{D}_{\tau }^{{\mu }_j}\left(\overline{u}(\overline{\mathit{\varrho }},\tau )\right)+{\delta }_3\overline{u}(\overline{\mathit{\varrho }},\tau )\right)\left)\right|,\hfill \\ \hfill & \le J_{\tau }^{{\nu }_i}\left(\right|{\delta }_4\left|\right|\nabla u(\overline{\mathit{\varrho }},\tau )-\nabla \overline{u}(\overline{\mathit{\varrho }},\tau )|+|{\delta }_5\left|\right|D_{\tau }^{\alpha }(\nabla u(\overline{\mathit{\varrho }},\tau ))-D_{\tau }^{\alpha }(\nabla \overline{u}(\overline{\mathit{\varrho }},\tau ))|\hfill \\ \hfill & +|{\delta }_1\left|\right|{\displaystyle \frac{\partial u(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}-{\displaystyle \frac{\partial \overline{u}(\overline{\mathit{\varrho }},\tau )}{\partial \tau }}|+|{\delta }_2\left|\right|D_{\tau }^{{\mu }_j}\left(u(\overline{\mathit{\varrho }},\tau )\right)-D_{\tau }^{{\mu }_j}\left(\overline{u}(\overline{\mathit{\varrho }},\tau )\right)|\hfill \\ \hfill & +|{\delta }_3\left|\right|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|+|{\delta }_7\left|\right|\mathsf{\Phi }(\overline{\mathit{\varrho }},\tau )\left|\right).\hfill \end{array}$$

(9)

Let $\delta =max\left(|{\delta }_1|,|{\delta }_2|,|{\delta }_3|,|{\delta }_4|,|{\delta }_5|,|{\delta }_6|,|{\delta }_7|\right)$ and $\left|\mathsf{\Phi }\right(\overline{\mathit{\varrho }},\tau \left)\right|<m_5,$ where $m_5>0$; we have

$$\begin{array}{cc}\hfill |u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|\le & \delta J_{\tau }^{{\nu }_i}(m_1|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|+m_1{m}_2|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|+m_3|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|\hfill \\ \hfill & +m_4|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|+|u(\overline{\mathit{\varrho }},\tau )-\overline{u}(\overline{\mathit{\varrho }},\tau )|+m_5),\hfill \end{array}$$

(10)

and taking the norm on both sides, we obtain

$$\begin{array}{cc}\hfill \parallel u-\overline{u}{\parallel }_{\infty }\le & \delta J_t^{{\nu }_i}\left(m_1\parallel u-\overline{u}{\parallel }_{\infty }+m_1{m}_2\parallel u-\overline{u}{\parallel }_{\infty }+m_3\parallel u-\overline{u}{\parallel }_{\infty }+m_4\parallel u-\overline{u}{\parallel }_{\infty }+{\parallel u-\overline{u}\parallel }_{\infty }+m_5\right),\hfill \\ \hfill & \le \delta {\displaystyle \frac{{\tau }^{{\nu }_i}}{\mathsf{\Gamma }(1+{\nu }_i)}}\left(\left(m_1+m_1{m}_2+m_3+m_4+1\right){\parallel u-\overline{u}\parallel }_{\infty }+m_5\right),\hfill \\ \hfill & \le \delta {\displaystyle \frac{T^{{\nu }_i}}{\mathsf{\Gamma }(1+{\nu }_i)}}\left(\left(m_1+m_1{m}_2+m_3+m_4+1\right){\parallel u-\overline{u}\parallel }_{\infty }+m_5\right),\mathrm{for}0\le \tau \le T.\hfill \end{array}$$

(11)

Finally, we obtain

$$\begin{array}{c}\hfill \parallel u-\overline{u}{\parallel }_{\infty }\le \mathcal{Q}{m}_5,\end{array}$$

(12)

where $\mathcal{Q}={\displaystyle \frac{\mathcal{H}}{1-\mathcal{M}}},$  $\mathcal{M}=\mathcal{H}(m_1+m_1{m}_2+m_3+m_4+1)$, and $\mathcal{H}=\delta {\displaystyle \frac{T^{{\nu }_i}}{\mathsf{\Gamma }(1+{\nu }_i)}},$ which proves the required result.    □

4. Proposed Numerical Method

4.1. Laplace Transform

This section focuses on the time discretization of the problem in Equation (4) via the LT. By employing the LT, we can transform the time-dependent problem to an equivalent time-independent problem. The LT of $u(\overline{\mathit{\varrho }},\tau )$ is defined as

$$\widehat{u}(\overline{\mathit{\varrho }},s)=\mathcal{L}\left\{u(\overline{\mathit{\varrho }},\tau )\right\}={\int }_0^{\infty }{e}^{-s\tau }u(\overline{\mathit{\varrho }},\tau )d\tau ,$$

while the LT of ${}_0{}^{c}D_{\tau }^{\mu }u(\overline{\mathit{\varrho }},\tau )$ and ${}_0{}^{c}D_{\tau }^{\nu }u(\overline{\mathit{\varrho }},\tau )$ is

$$\mathcal{L}\left\{{}_0{}^{c}D_{\tau }^{{\mu }_j}u(\overline{\mathit{\varrho }},\tau )\right\}=s^{{\mu }_j}\widehat{u}(\overline{\mathit{\varrho }},s)-s^{{\mu }_j-1}u(\overline{\mathit{\varrho }},0),{\mu }_j\in (0,1).$$

and

$$\mathcal{L}\left\{{}_0{}^{c}D_{\tau }^{{\nu }_i}u(\overline{\mathit{\varrho }},\tau )\right\}=s^{{\nu }_i}\widehat{u}(\overline{\mathit{\varrho }},s)-s^{{\nu }_i-1}u(\overline{\mathit{\varrho }},0)-s^{{\nu }_i-2}{u}_{\tau }(\overline{\mathit{\varrho }},0),{\nu }_i\in (1,2)$$

The LT of Equation (4) implies

$$\begin{array}{cc}\hfill & \mathcal{L}\left[\sum _{i=1}^k{c}_{1,i}{D}_{\tau }^{{\nu }_i}u(\overline{\mathit{\varrho }},\tau )+c_2{D}_{\tau }u(\overline{\mathit{\varrho }},\tau )+\sum _{j=1}^l{c}_{3,j}{D}_{\tau }^{{\mu }_j}u(\overline{\mathit{\varrho }},\tau )+c_{4}u(\overline{\mathit{\varrho }},\tau )\right]\hfill \\ \hfill & =\mathcal{L}\left[c_5{\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},\tau )+c_6{D}_{\tau }^{\alpha }\left({\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},\tau )\right)+f(\overline{\mathit{\varrho }},\tau )\right]\hfill \end{array}$$

and

$$\mathcal{L}\left[{\mathcal{L}}_{\mathit{b}}u(\overline{\mathit{\varrho }},\tau )\right]=\mathcal{L}\left[{\psi }_1(\overline{\mathit{\varrho }},\tau )\right],$$

which implies

$$\begin{array}{cc}\hfill & \sum _{i=1}^k{c}_{1,i}\left(s^{{\nu }_i}\widehat{u}(\overline{\mathit{\varrho }},s)-s^{{\nu }_i-1}u(\overline{\mathit{\varrho }},0)-s^{{\nu }_i-2}{u}_{\tau }(\overline{\mathit{\varrho }},0)\right)+c_2\left(s\widehat{u}(\overline{\mathit{\varrho }},s)-u(\overline{\mathit{\varrho }},0)\right)+\hfill \\ \hfill & \sum _{j=1}^l{c}_{3,j}\left(s^{{\mu }_j}\widehat{u}(\overline{\mathit{\varrho }},s)-s^{{\mu }_j-1}u(\overline{\mathit{\varrho }},0)\right)+c_4\widehat{u}(\overline{\mathit{\varrho }},s)=c_5{\mathcal{L}}_{\mathit{g}}\widehat{u}(\overline{\mathit{\varrho }},s)+\hfill \\ \hfill & c_6\left(s^{\alpha }{\mathcal{L}}_{\mathit{g}}\widehat{u}(\overline{\mathit{\varrho }},s)-s^{\alpha -1}{\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},0)\right)+\widehat{f}(\overline{\mathit{\varrho }},s)\hfill \end{array}$$

and

$${\mathcal{L}}_{\mathit{b}}\widehat{u}(\overline{\mathit{\varrho }},s)={\psi }_1(\overline{\mathit{\varrho }},s),$$

Simplifying the above system, we have

$$\left(\sum _{i=1}^k{c}_{1,i}{s}^{{\nu }_i}I+c_{2}sI+\sum _{j=1}^l{c}_{3,j}{s}^{{\mu }_j}I+c_{4}I-c_5{\mathcal{L}}_{\mathit{g}}-c_6{s}^{\alpha }{\mathcal{L}}_{\mathit{g}}\right)\widehat{u}(\overline{\mathit{\varrho }},s)=\widehat{F}(\overline{\mathit{\varrho }},s),$$

(13)

$${\mathcal{L}}_{\mathit{b}}\widehat{u}(\overline{\mathit{\varrho }},s)={\psi }_1(\overline{\mathit{\varrho }},s),$$

(14)

where

$$\begin{array}{cc}\hfill & \widehat{F}(\overline{\mathit{\varrho }},s)=\sum _{i=1}^k{c}_{1,i}\left(s^{{\nu }_i-1}u(\overline{\mathit{\varrho }},0)+s^{{\nu }_i-2}{u}_{\tau }(\overline{\mathit{\varrho }},0)\right)+c_{2}u(\overline{\mathit{\varrho }},0)+\hfill \\ \hfill & \sum _{j=1}^l{c}_{3,j}\left(s^{{\mu }_j-1}u(\overline{\mathit{\varrho }},0)\right)-c_6{s}^{\alpha -1}{\mathcal{L}}_{\mathit{g}}u(\overline{\mathit{\varrho }},0)+\widehat{f}(\overline{\mathit{\varrho }},s),\hfill \end{array}$$

where I denotes the identity operator and ${\mathcal{L}}_{\mathit{g}}=\mathsf{\Delta }$. In the suggested scheme, the operators ${\mathcal{L}}_{\mathit{g}}$ and ${\mathcal{L}}_{\mathit{b}}$ are first discretized via the LMM, then the system in (13) and (14) is solved in parallel for each point s in Laplace space (for example, see [35,36]). Finally, the inverse Laplace transform is utilized to obtained the approximate solution to the actual problem defined in Equation (4). The next section provides a detailed description of the LMM.

4.2. Localized Meshless Method

In the LMM, the considered boundary value problem is interpolated on $N_{glob}$ points, the distribution of which may be uniform or non-uniform. For each point ${\overline{\mathit{\zeta }}}_i(i=1,2,\dots ,N_{glob})$ in the subdomain, ${\mathsf{\Omega }}_i={\left\{{\overline{\mathit{\varrho }}}_j^i\right\}}_{j=1}^{n_{loc}}(i=1,2,\dots ,N_{glob}),$ $n_{loc}$ denotes the nodes in ${\mathsf{\Omega }}_i.$ The interpolant for $\widehat{u}\left(\overline{\mathit{\varrho }}\right)$ via LMM has the form

$$\widehat{u}\left({\overline{\mathit{\varrho }}}_i\right)=\sum _{j=1}^{n_{loc}}{\chi }_j^i\kappa (\parallel {\overline{\mathit{\varrho }}}_i-{\overline{\mathit{\varrho }}}_j^i\parallel ),$$

(15)

where ${\mathit{\chi }}^i={\left\{{\chi }_j^i\right\}}_{j=1}^{N_{glob}}$ are the unknown coefficients, $\kappa \left(r\right)$ is an RBF, and $r= \parallel {\overline{\mathit{\varrho }}}_i-{\overline{\mathit{\varrho }}}_j^i\parallel$. In the literature, there are several RBFs; here, we use multiquadric radial basis function (MQRBF) and Gaussian radial basis function (GSRBF), respectively defined as follows:

$$\kappa \left(r\right)=\sqrt{1+{\epsilon }^2{r}^2},$$

and

$$\kappa \left(r\right)=e^{-{\epsilon }^2{r}^2},$$

where $\epsilon$ denotes the shape parameter. Every point ${\overline{\mathit{\varrho }}}_j^i$ and its $n_{loc}-1$ neighboring points are called a stencil. We obtain the matrix $n_{loc}\times n_{loc}$ at each stencil as follows:

$${\widehat{\mathit{u}}}^i={\mathit{B}}^i{\mathit{\chi }}^i,i=1,2,3,\dots ,N_{glob}$$

(16)

where ${\widehat{\mathit{u}}}^i={[u\left({\overline{\mathit{\varrho }}}_1^i\right),u\left({\overline{\mathit{\varrho }}}_2^i\right),\dots ,u\left({\overline{\mathit{\varrho }}}_{n_{loc}}^i\right)]}^T,$ ${\mathit{B}}^i=\left[\kappa \right(\parallel {\overline{\mathit{\varrho }}}_k^i-{\overline{\mathit{\varrho }}}_j^i{\parallel \left)\right]}_{1\le k,j\le n_{loc}},$ and ${\mathit{\chi }}^i=[{\chi }_1^i,$ ${\chi }_2^i,\dots ,{\chi }_{n_{loc}}^i{]}^T.$

The coefficients ${\mathit{\chi }}^i$ in Equation (16) are computed as

$${\mathit{\chi }}^i={\left({\mathit{B}}^i\right)}^{-1}{\widehat{\mathit{u}}}^i,$$

(17)

for which the entries are

$$\left(\begin{array}{c}{\chi }_1^i\\ {\chi }_2^i\\ .\\ .\\ .\\ {\chi }_{n_{loc}}^i\end{array}\right)={\left(\begin{array}{cccccc}\kappa (\parallel {\overline{\mathit{\varrho }}}_1^i-{\overline{\mathit{\varrho }}}_1^i\parallel )& \kappa (\parallel {\overline{\mathit{\varrho }}}_1^i-{\overline{\mathit{\varrho }}}_2^i\parallel )& .& .& .& \kappa (\parallel {\overline{\mathit{\varrho }}}_1^i-{\overline{\mathit{\varrho }}}_{n_{loc}}^i\parallel )\\ \kappa (\parallel {\overline{\mathit{\varrho }}}_2^i-{\overline{\mathit{\varrho }}}_1^i\parallel )& \kappa (\parallel {\overline{\mathit{\varrho }}}_2^i-{\overline{\mathit{\varrho }}}_2^i\parallel )& .& .& .& \kappa (\parallel {\overline{\mathit{\varrho }}}_2^i-{\overline{\mathit{\varrho }}}_{n_{loc}}^i\parallel )\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ \kappa (\parallel {\overline{\mathit{\varrho }}}_{n_{loc}}^i-{\overline{\mathit{\varrho }}}_1^i\parallel )& \kappa (\parallel {\overline{\mathit{\varrho }}}_{n_{loc}}^i-{\overline{\mathit{\varrho }}}_2^i\parallel )& .& .& .& \kappa (\parallel {\overline{\mathit{\varrho }}}_{n_{loc}}^i-{\overline{\mathit{\varrho }}}_{n_{loc}}^i\parallel )\end{array}\right)}^{-1}\left(\begin{array}{c}\widehat{u}\left({\overline{\mathit{\varrho }}}_1^i\right)\\ \widehat{u}\left({\overline{\mathit{\varrho }}}_2^i\right)\\ .\\ .\\ .\\ \widehat{u}\left({\overline{\mathit{\varrho }}}_{n_{loc}}^i\right)\end{array}\right).$$

The linear operator ${\mathcal{L}}_{\mathit{g}}$ has an approximation of the form

$${\mathcal{L}}_{\mathit{g}}\widehat{u}=\sum _{j=1}^{n_{loc}}{\chi }_j^i{\mathcal{L}}_{\mathit{g}}\kappa (\parallel {\overline{\mathit{\varrho }}}_i-{\overline{\mathit{\varrho }}}_j^i\parallel ),$$

which can be written as

$${\mathcal{L}}_{\mathit{g}}\widehat{u}={\mathit{\psi }}^i{\left({\mathit{B}}^i\right)}^{-1}{\widehat{\mathit{u}}}^i,$$

where ${\mathit{\psi }}^i= [{\mathcal{L}}_{\mathit{g}}\kappa (\parallel {\overline{\mathit{\varrho }}}_k^i-{\overline{\mathit{\varrho }}}_j^i{\parallel \left)\right]}_{1\le k,j\le n_{loc}}.$ Hence, we have

$${\mathcal{L}}_{\mathit{g}}\widehat{u}={\mathbf{\Psi }}^i{\widehat{\mathit{u}}}^i,$$

(18)

in which ${\mathbf{\Psi }}^i={\mathit{\psi }}^i{\left({\mathit{B}}^i\right)}^{-1}.$ The global form of Equation (18), may be obtained by extending ${\mathbf{\Psi }}^i$ to $\mathbf{\Theta }$ by adding zeros to the corresponding positions in each row. Hence, at each node ${\overline{\mathit{\varrho }}}_i$, the operator ${\mathcal{L}}_{\mathit{g}}$ is approximated as

$${\mathcal{L}}_{\mathit{g}}\widehat{u}=\mathbf{\Theta }\widehat{u}.$$

(19)

The matrix ${\mathbf{\Theta }}_{N_{glob}\times N_{glob}}$ is sparse, with $N_{glob}-n_{loc}$ zero entries and $n_{loc}$ nonzero entries. The boundary operator ${\mathcal{L}}_{\mathit{b}}$ can be estimated using the same process:

$${\mathcal{L}}_{\mathit{b}}\widehat{u}=\mathbf{\Phi }\widehat{u}.$$

(20)

By substituting Equations (19) and (20) into Equations (13) and (14), we obtain the following fully discrete system:

$$\left(\sum _{i=1}^k{c}_{1,i}{s}^{{\nu }_i}I+c_{2}sI+\sum _{j=1}^l{c}_{3,j}{s}^{{\mu }_j}I+c_{4}I-c_5\mathbf{\Theta }-c_6{s}^{\alpha }\mathbf{\Theta }\right)\widehat{u}(\overline{\mathit{\varrho }},s)=\widehat{F}(\overline{\mathit{\varrho }},s),$$

(21)

$$\mathbf{\Phi }\widehat{u}(\overline{\mathit{\varrho }},s)={\psi }_1(\overline{\mathit{\varrho }},s).$$

(22)

By solving the system in (21) and (22) at each node s in LT domain, we obtain the desired solution $\widehat{u}(\overline{\mathit{\varrho }},s).$

Shape Parameter

The optimal value of the parameter $\epsilon$ is derived using the uncertainty principle following [37]. The key steps are outlined in Algorithm 1.

First, ${\left({\mathit{B}}^i\right)}^{-1}$ is evaluated via the svd as ${\left({\mathit{B}}^i\right)}^{-1}={\left({\mathbf{USV}}^{\mathbf{T}}\right)}^{-1}={\mathbf{VS}}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}$ (see  [38]). Thus, we are able to evaluate ${\mathbf{\Psi }}^i$ in (18).

Algorithm 1 Algorithm for selecting the optimal shape parameter

1: Input: ${{\mathit{CN}}_D}_{min},{{\mathit{CN}}_D}_{max},{\epsilon }_{Increment}$ 2: i: ${\mathit{CN}}_D=1$ 3: ii: choose ${10}^{+12}<{\mathit{CN}}_D<{10}^{+16}$ 4: iii: $\mathrm{while}{\mathit{CN}}_D>{{\mathit{CN}}_D}_{max} \mathrm{and} {\mathit{CN}}_D<{{\mathit{CN}}_D}_{min}$ 5: iv: Construct the system matrix ${\mathit{B}}^i$ 6: v: $\mathbf{U},\mathbf{S},\mathbf{V}=svd\left({\mathit{B}}^i\right)$ 7: vi: ${\mathit{CN}}_D={\displaystyle \frac{{\rho }_{max}}{{\rho }_{min}}}$ 8: vii: $\mathrm{if}{\mathit{CN}}_D<{{\mathit{CN}}_D}_{min},\epsilon =\epsilon -{\epsilon }_{Increment}$ 9: viii: $\mathrm{if}{\mathit{CN}}_D>{{\mathit{CN}}_D}_{max},\epsilon =\epsilon +{\epsilon }_{Increment}$ 10: Output: $\epsilon \left(\mathrm{Best}\right)=\epsilon$.

4.3. Inverse Laplace Transform

We obtain the solution to the problem in (4) using the inversion of LT, as follows:

$$u(\overline{\mathit{\varrho }},\tau )={\displaystyle \frac{1}{2\pi i}}{\int }_{\rho -i\infty }^{\rho +i\infty }{e}^{s\tau }\widehat{u}(\overline{\mathit{\varrho }},s)ds={\displaystyle \frac{1}{2\pi i}}{\int }_{\mathcal{H}}{e}^{s\tau }\widehat{u}(\overline{\mathit{\varrho }},s)ds,\rho >{\rho }_0$$

(23)

where the integral (23) is known as the Bromwich integral and $\mathcal{H}$ represents the contour of integration. When dealing with a complex function $\widehat{u}(\overline{\mathit{\varrho }},s)$, approximating the Bromwich integral as defined in Equation (23) can be quite difficult; for example, the closed-form solution to Equation (23) is generally not obtainable if the transform $\widehat{u}(\overline{\mathit{\varrho }},s)$ is highly oscillatory. Many approaches are available in the literature for approximating Equation (23) [39]. In this work, we use the Stehfest method (STM) and the improved Talbot method (ITM).

4.3.1. Improved Talbot Method (ITM)

The ITM employs numerical quadrature to evaluate the integral in (23). Two effective techniques utilized in this context are the midpoint and trapezoidal rule, both of which are applied in conjunction with contour deformation [40]. The main objective of deforming the contour is to effectively handle the factor $e^{s\tau }$. Specifically, the integration path can be adjusted to form a Hankel contour, which starts at $-\infty$ in the third quadrant, encircles all singularities of the transform function, and returns to $-\infty$ in the second quadrant. Because of this arrangement, the exponential factor can decline rapidly, rendering the integral in (23) suitable for approximation via the trapezoidal and midpoint techniques. This deformation is supported by Cauchy’s theorem as long as the contour stays inside the area where the transform function $\widehat{u}(\overline{\mathit{\varrho }},s)$ is analytic. Moreover, some gentle constraints on the decay behavior of $\widehat{u}(\overline{\mathit{\varrho }},s)$ in the left half plane are needed [27]. In the context of the improved Talbot method, we present the Hankel contour in parametric form as follows [27]:

$$\mathcal{H}:s=s\left(\xi \right),-\pi \le \xi \le \pi ,Res(\pm \pi )=-\infty$$

(24)

obtaining

$$s\left(\xi \right)={\displaystyle \frac{M_T}{\tau }}\zeta \left(\xi \right),\zeta \left(\xi \right)=-{\kappa }_1+\sigma \xi cot\left(\mu \xi \right)+\nu i\xi ,$$

(25)

where $\mu ,\nu ,\kappa$ are parameters. Substituting Equation (25) into Equation (23), we obtain

$$u(\overline{\mathit{\varrho }},\tau )={\displaystyle \frac{1}{2\pi i}}{\int }_{-\pi }^{\pi }{e}^{s\left(\xi \right)\tau }\widehat{u}(\overline{\mathit{\varrho }},s\left(\xi \right))s^{\prime }\left(\xi \right)d\xi .$$

(26)

The midpoint rule with a uniform step size of $h={\displaystyle \frac{2\pi }{M_T}}$ is used for approximating Equation (26), as follows:

$$\tilde{u}(\overline{\mathit{\varrho }},\tau )\approx {\displaystyle \frac{1}{M_{T}i}}\sum _{k=1}^{M_T}{e}^{s\left({\xi }_k\right)\tau }\widehat{u}(\overline{\mathit{\varrho }},s\left({\xi }_k\right))s^{\prime }\left({\xi }_k\right),{\xi }_k=-\pi +(k-{\displaystyle \frac{1}{2}})h.$$

(27)

4.3.2. Error Analysis

The subsequent theorem establishes the basis for conducting the error analysis of ITM.

Theorem 2

([27]). Let ${\xi }_k$ be defined as in (27). Consider $f:\mathsf{\Lambda }\to \mathbb{C}$ as an analytic function on the set

$$\mathsf{\Lambda }=\{\xi \in \mathbb{C}:-\pi <Re\xi <\pi ,and-d<Im\xi <c\}.$$

When $c,d>0,$ the following holds:

$${\int }_{-\pi }^{\pi }f\left(\xi \right)d\xi -{\displaystyle \frac{2\pi }{M_T}}\sum _{j=1}^{M_T}f\left({\xi }_k\right)=H_{-}\left(\sigma \right)+H_{+}\left(\varsigma \right),$$

where

$$H_{+}\left(\sigma \right)={\displaystyle \frac{1}{2}}\left({\int }_{-\pi }^{-\pi +i\sigma }+{\int }_{-\pi +i\sigma }^{\pi +i\sigma }+{\int }_{\pi +i\sigma }^{\pi }\right)\left(1+itan\left({\displaystyle \frac{M_T\xi }{2}}\right)\right)f\left(\xi \right)d\xi$$

and

$$H_{-}\left(\varsigma \right)={\displaystyle \frac{1}{2}}\left({\int }_{-\pi }^{-\pi -i\varsigma }+{\int }_{-\pi -i\varsigma }^{\pi -i\varsigma }+{\int }_{\pi -i\varsigma }^{\pi }\right)\left(1-itan\left({\displaystyle \frac{M_T\xi }{2}}\right)\right)f\left(\xi \right)d\xi ,$$

with $\forall 0<\sigma <c, and0<\varsigma <dand{M}_T$ being even; for odd $M_T$, we can substitute $tan\left({\displaystyle \frac{M_T\xi }{2}}\right)$ with $-cot\left({\displaystyle \frac{M_T\xi }{2}}\right)$ if $f\left(\xi \right)$ is real-valued, meaning that $f\left(\overline{\xi }\right)=\overline{f\left(\xi \right)}$. If c and d are taken to be equal, then

$$H\left(\varsigma \right)=H_{+}\left(\varsigma \right)+H_{-}\left(\varsigma \right)=Re{\int }_{-\pi +i\varsigma }^{\pi +i\varsigma }\left(1+itan\left({\displaystyle \frac{M_T\xi }{2}}\right)\right)f\left(\xi \right)d\xi .$$

Thus, we can conclude that the complex tangent function may be bounded in the following ways by examining its behavior:

$$\left|H\left(\varsigma \right)\right|\le {\displaystyle \frac{4\pi \mathcal{C}}{exp\left(c{M}_T\right)-1}}.$$

The above-described procedure is applicable even for $M_T$, $\mathcal{C},c\in {\mathbb{R}}^{+}$, and a similar approach can be applied for odd $M_T.$

To obtain the best results, it is important to identify the ideal contour of integration, which can be achieved by using the optimal values of the parameters in Equation (25). The following values have been proposed by the authors of [27] as the optimal parameters:

$${\kappa }_1=0.61220,\sigma =0.50170,\nu =0.26450,\mathrm{and}\mu =0.64070.$$

With the above optimal parameters, the error estimate is of $O\left(e^{-1.3580M_T}\right).$

4.4. Stehfest Method (STM)

The Gaver–Stehfest approach, which was developed in the late 1960s, has a rich history as one of the most efficient and simple methods for numerical inversion of the LT. Due to its efficiency and simplicity, it has gained widespread use in numerous domains such as finance, computational physics, chemistry, and economics. In 1966, Gaver [41] introduced simple though slowly converging approximations for the inverse LT. Later, in 1970, Stehfest [28] applied convergence acceleration techniques to Gaver’s method, giving rise to the Gaver–Stehfest method. The Gaver–Stehfest algorithm approximates the function $u(\overline{\mathit{\varrho }},\tau )$ using a series of functions, as follows:

$$\tilde{u}(\overline{\mathit{\varrho }},\tau )={\displaystyle \frac{ln2}{\tau }}\sum _{j=1}^{M_S}{\alpha }_j\widehat{u}\left(\overline{\mathit{\varrho }},{\displaystyle \frac{ln2}{\tau }}j\right)$$

(28)

where ${\alpha }_j$ are provided as follows:

$${\alpha }_j={(-1)}^{{\displaystyle \frac{M_S}{2}}+j}\sum _{\gamma =\lfloor {\displaystyle \frac{j+1}{2}}\rfloor }^{min(j,{\displaystyle \frac{M_S}{2}})}{\displaystyle \frac{{\gamma }^{\frac{M_S}{2}}\left(2\gamma \right)!}{(\frac{M_S}{2}-\gamma )!\gamma !(\gamma -1)!(j-\gamma )!(2\gamma -j)!}}.$$

(29)

The system defined in Equations (21) and (22) is solved using the LT parameters $s={\left\{{\displaystyle \frac{ln2}{\tau }}j\right\}}_{j=1}^{M_S}$. The numerical solution $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ of Equation (4) can be derived using Equation (28). The Gaver–Stehfest algorithm has gained significant popularity among researchers due to several qualities: first, it is linear and exact for constant functions, and all the coefficients can be computed manually; most importantly, the algorithm works without requiring complex numbers, and relies solely on the values of the LT along the positive real line. In the literature, this algorithm has been utilized by numerous practitioners, including [42,43], showing that the technique converges rapidly to $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ provided that $u(\overline{\mathit{\varrho }},\tau )$ exhibits non-oscillatory behavior.

Convergence of STM

In [42], the author established the convergence of $\tilde{u}(\overline{\mathit{\varrho }},\tau )$, with the results are founded on the following theorem.

Theorem 3.

Assume that $u:(0,\infty )\to \mathbb{R}$ is a locally integrable function such that the LT $\widehat{u}(\overline{\mathit{\varrho }},s)$ exists for all $s>0,$ where s is the Laplace parameter. Moreover, let $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ represent the approximate solution provided by Equation (28). Then:

1.

$\tilde{u}(\overline{\mathit{\varrho }},\tau )$ converges for a given $u(\overline{\mathit{\varrho }},\tau )$ as τ approaches a specific value.

2.

Let σ and ς be real numbers such that $0<\varsigma <0.25$, and consider the integral

$${\int }_0^{\varsigma }|u(\overline{\mathit{\varrho }},-\tau lo{g}_2(1/2+\xi ))+u(\overline{\mathit{\varrho }},-\tau lo{g}_2(1/2-\xi ))-2\sigma |{\xi }^{-1}d\xi <\infty ;$$

then, as $M_S\to +\infty ,$ it follows that $\tilde{u}(\overline{\mathit{\varrho }},\tau )\to \sigma .$

3.

Assume that $u(\overline{\mathit{\varrho }},\tau )$ is of bounded variation in the vicinity of τ; then, we have

$$\tilde{u}(\overline{\mathit{\varrho }},\tau )\to {\displaystyle \frac{u(\overline{\mathit{\varrho }},\tau +0)+u(\overline{\mathit{\varrho }},\tau -0)}{2}},\mathit{as}{M}_S\to +\infty .$$

Corollary 1.

Based on the assumptions in the above theorem, if 

$$u(\overline{\mathit{\varrho }},\tau +\xi )-u(\overline{\mathit{\varrho }},\tau )=O\left({\left|\xi \right|}^{\vartheta }\right),$$

$\forall \xi$ and for some $\vartheta ,$ then it follows that $\tilde{u}(\overline{\mathit{\varrho }},\tau )\to u(\overline{\mathit{\varrho }},\tau )$, as $M_S\to +\infty .$

Furthermore, in [44] the authors performed various experiments to study the impact of parameters on the numerical scheme’s accuracy. They concluded that “for a specific number of significant digits ${\eta }_1$, one should select $M_S\in {\mathbb{Z}}^{+}$ such that $M_S=\lceil 2.2{\eta }_1\rceil$. Set the precision of the system to ${\eta }_2=\lceil 1.1M_S\rceil$, then compute ${\left\{{\alpha }_j\right\}}_{j=1}^{M_S}$ using (29). Subsequently, calculate $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ using the transformed function and $\tau$ as provided in (28)”. These findings suggest that the error is bounded by

$${10}^{-({\eta }_1+1)}\le {\displaystyle \frac{\tilde{u}(\overline{\mathit{\varrho }},\tau )-u(\overline{\mathit{\varrho }},\tau )}{u(\overline{\mathit{\varrho }},\tau )}}\le {10}^{-{\eta }_1},$$

where $M_S=\lceil 2.2{\eta }_1\rceil$ [45].

The key steps of the LT-based LMM are outlined in the following Algorithm 2.    

Algorithm 2 The numerical scheme

1: Input: The domain, order of the Caputo derivative, initial-boundary conditions, nonhomogeneous term, initial shape parameter, best values of the parameters for the inverse LT method, and final time. 2: Step: Apply the LT to Equation (4) to derive the time-independent elliptic problems defined in Equations (13) and (14). 3: Step: Discretize the operators ${\mathcal{L}}_{\mathit{g}}$ and ${\mathcal{L}}_{\mathit{b}}$ via the LMM, obtain its approximations $\mathbf{\Theta }$ and $\mathsf{\Phi }$, and obtain the fully discrete system in (21) and (22). 4: Step: Solve the system in (21) and (22) for every point “s” in parallel and obtain the LT domain for the numerical solution $\widehat{u}(\overline{\mathit{\varrho }},s)$. 5: Step: Use the numerical inversion technique in (27) or (28) to derive the approximate solution to Equation (4). 6: Step: $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ is the approximate solution.

5. Stability

In this section, we discuss the stability of the fully discrete system. The system in (21) and (22) is written as follows:

$$\mathcal{G}\widehat{\mathit{u}}=\mathit{Q}$$

(30)

where G is the interpolation matrix obtained via the LMM. The stability constant for the system in (30) has the following form:

$$C_S=\underset{\widehat{\mathit{u}}\ne 0}{sup}{\displaystyle \frac{\parallel \widehat{\mathit{u}}\parallel }{\parallel \mathcal{G}\widehat{\mathit{u}}\parallel }}$$

(31)

where $C_S$ is finite for any norm $\parallel .\parallel$ on ${\mathbb{R}}^{N_{glob}}$. From Equation (31), we obtain

$${\parallel \mathcal{G}\parallel }^{-1}\le {\displaystyle \frac{\parallel \widehat{\mathit{u}}\parallel }{\parallel \mathcal{G}\widehat{\mathbf{u}}\parallel }}\le C_S.$$

(32)

We may also write

$$\parallel {\mathcal{G}}^{†}\parallel =\underset{\mathit{\varpi }\ne 0}{sup}{\displaystyle \frac{\parallel {\mathcal{G}}^{†}\mathit{\varpi }\parallel }{\parallel \mathit{\varpi }\parallel }},$$

(33)

where ${\mathcal{G}}^{†}$ is the pseudoinverse of $\mathcal{G}.$ As a result,

$$\parallel {\mathcal{G}}^{†}\parallel \ge \underset{\mathit{\varpi }=\mathcal{G}\widehat{\mathit{u}}\ne 0}{sup}{\displaystyle \frac{\parallel {\mathcal{G}}^{†}\mathcal{G}\widehat{\mathit{u}}\parallel }{\parallel \mathcal{G}\widehat{\mathit{u}}\parallel }}=\underset{\widehat{\mathit{u}}\ne 0}{sup}{\displaystyle \frac{\parallel \widehat{\mathit{u}}\parallel }{\parallel \mathcal{G}\widehat{\mathit{u}}\parallel }}=C_S.$$

(34)

The bounds for the stability constant $C_S$ are provided in Equations (32) and (33). For numerical approximation of the system in Equation (30), the pseudoinverse may be challenging to compute; however, it ensures the system’s numerical stability. The MATLAB command “condest” provides an estimate of the $L_{\infty }$-norm of ${\parallel {\mathcal{G}}^{-1}\parallel }_{\infty };$ thus, we have

$$C_S={\displaystyle \frac{condest\left({\mathcal{G}}^{\prime }\right)}{{\parallel \mathcal{G}\parallel }_{\infty }}}.$$

(35)

This approach works effectively within a short time frame for our sparse interpolation matrix $\mathcal{G}$. The plots of the stability constant $C_S$ are depicted in Figure 1a,b.

6. Examples

This section presents the simulation results, assessing the performance of the method through four examples. We conducted the numerical experiments used MATLAB R2019a on a Windows 10 (64-bit) PC configured with an Intel(R) Core(TM) i5-3317U 1.70 GHz CPU and 12 GB of RAM. The performance of the LMM is quantified using four error norms: the absolute error, root mean square error, maximum absolute error, and relative $L2$ error, respectively defined as follows:

$$L_{Abs}=|u({\overline{\mathit{\varrho }}}_q,\tau )-\tilde{u}({\overline{\mathit{\varrho }}}_q,\tau )|,$$

$$RMS=\sqrt{{\displaystyle \frac{{\sum }_{q=1}^{N_{glob}}{\left(u({\overline{\mathit{\varrho }}}_q,\tau )-\tilde{u}({\overline{\mathit{\varrho }}}_q,\tau )\right)}^2}{N_{glob}}}},$$

$$L_{\infty }=\underset{1\le q\le N_{glob}}{max}|u({\overline{\mathit{\varrho }}}_q,\tau )-\tilde{u}({\overline{\mathit{\varrho }}}_q,\tau )|,$$

$$L_2=\sqrt{{\displaystyle \frac{{\sum }_{q=1}^{N_{glob}}{\left(u({\overline{\mathit{\varrho }}}_q,\tau )-\tilde{u}({\overline{\mathit{\varrho }}}_q,\tau )\right)}^2}{{\sum }_{q=1}^{N_{glob}}{\left(u({\overline{\mathit{\varrho }}}_q,\tau )\right)}^2}}},$$

where $u(\overline{\mathit{\varrho }},\tau )$ and $\tilde{u}(\overline{\mathit{\varrho }},\tau )$ are the exact and approximation solutions, respectively, $N_{glob}$ denotes the number of nodes in the global domain, $n_{loc}$ denotes the number of nodes in the local domain, $M_S$ represents the nodes used in STM, and $M_T$ represents the nodes used in ITM. The problem is solved using the following parameters: $c_{1,1}=c_2=c_{3,1}=c_4=c_5=c_6=1.$ The source term along with the initial and boundary conditions for each example is determined using the exact solution.

6.1. Example 1

In the first example, we examine Equation (4) with analytical solution $u(x,\tau )={\tau }^{(3+\mu +\alpha +\nu )}{e}^{-x^2}.$ The results obtained using the LMM based on MQRBF coupled with ITM and STM for various values of $N_{glob},n_{loc},M_S,$ and $M_T$ are presented in

Table 1. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 1, with $\nu =1.6,\mu =0.7,\alpha =0.8$.

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	250	3	20	2.09$\times {10}^{-2}$	1.71$\times {10}^{-3}$	8.37$\times {10}^{-5}$	1.14	0.337616
			22	2.46$\times {10}^{-3}$	2.00$\times {10}^{-4}$	9.87$\times {10}^{-6}$	1.14	0.338195
			24	2.96$\times {10}^{-4}$	2.32$\times {10}^{-5}$	1.18$\times {10}^{-6}$	1.14	0.375847
	150	5	24	5.06$\times {10}^{-4}$	5.93$\times {10}^{-5}$	3.37$\times {10}^{-6}$	7.38	0.229518
	170			7.44$\times {10}^{-4}$	8.12$\times {10}^{-5}$	4.37$\times {10}^{-6}$	8.37	0.196635
	185			9.55$\times {10}^{-4}$	9.93$\times {10}^{-5}$	5.16$\times {10}^{-6}$	9.11	0.216754
STM			$M_S$
	250	3	16	1.26$\times {10}^{-1}$	1.03$\times {10}^{-2}$	5.07$\times {10}^{-4}$	3.6	0.257375
			18	1.19$\times {10}^{-2}$	9.91$\times {10}^{-4}$	4.76$\times {10}^{-5}$	3.6	0.327381
			20	1.74$\times {10}^{-3}$	2.09$\times {10}^{-4}$	6.97$\times {10}^{-6}$	3.6	0.281991
	150	5	18	8.45$\times {10}^{-3}$	8.80$\times {10}^{-4}$	5.63$\times {10}^{-5}$	4.13	0.182686
	170			9.58$\times {10}^{-3}$	9.65$\times {10}^{-4}$	5.63$\times {10}^{-5}$	4.69	0.184717
	195			9.90$\times {10}^{-3}$	9.07$\times {10}^{-4}$	5.07$\times {10}^{-5}$	5.38	0.201450
[18]				1.40$\times {10}^{-3}$

. Similarly, the results obtained using the LMM-based GSRBF coupled with STM and ITM for various values of $N_{glob},n_{loc},M_S,$ and $M_T$ are presented in

Table 2. Obtained errors using the LMM based on GSRBF coupled with ITM and STM for Example 1, with $\nu =1.6,\mu =0.7,\alpha =0.8$.

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	250	3	20	2.09$\times {10}^{-2}$	1.71$\times {10}^{-3}$	8.36$\times {10}^{-5}$	9.6	0.318618
			22	2.43$\times {10}^{-3}$	1.98$\times {10}^{-4}$	9.74$\times {10}^{-6}$	9.6	0.351536
			24	2.63$\times {10}^{-4}$	2.13$\times {10}^{-5}$	1.05$\times {10}^{-6}$	9.6	0.331514
	150	5	24	6.55$\times {10}^{-4}$	6.97$\times {10}^{-5}$	4.37$\times {10}^{-6}$	8.27	0.174798
	170			8.46$\times {10}^{-4}$	8.52$\times {10}^{-5}$	4.98$\times {10}^{-6}$	9.38	0.212009
	185			9.04$\times {10}^{-4}$	8.76$\times {10}^{-5}$	4.89$\times {10}^{-6}$	10	0.212871
STM			$M_S$
	250	3	16	1.26$\times {10}^{-1}$	1.03$\times {10}^{-2}$	5.06$\times {10}^{-4}$	3	0.258385
			18	1.09$\times {10}^{-2}$	9.05$\times {10}^{-4}$	4.39$\times {10}^{-5}$	3	0.250448
			20	1.85$\times {10}^{-3}$	1.68$\times {10}^{-4}$	7.43$\times {10}^{-6}$	3	0.271159
	150	5	18	9.26$\times {10}^{-3}$	9.66$\times {10}^{-4}$	6.17$\times {10}^{-5}$	4.65	0.182282
	170			9.42$\times {10}^{-3}$	9.31$\times {10}^{-4}$	5.54$\times {10}^{-5}$	5.28	0.180676
	195			9.92$\times {10}^{-3}$	9.58$\times {10}^{-4}$	5.36$\times {10}^{-5}$	5.74	0.200793
[18]				1.40$\times {10}^{-3}$

. The numerical solution is compared to the exact solution in Figure 2a, while the profiles of the numerical solution for $x\in [0,2]$ and $\tau =[0.6,0.7,0.8,0.9,1]$ are presented in Figure 2b, showing that the solution profile exhibits asymptotic behavior.

Furthermore, we compare the $L_2,$ $L_{\infty },$ and $RMS$ errors using the LMM based on MQRBF and GSRBF coupled with ITM and STM for LT inversion with the $M_T$ and $M_S$ values and $N_{glob}=250,n_{loc}=5$ in Figure 3a,b. The surface plot of $L_{Abs}$ using the LMM based on MQRBF coupled with ITM is presented in Figure 4a, while that using the LMM based on GSRBF coupled with STM is presented in Figure 4b.

Figure 5a, Figure 5b, Figure 6a, Figure 6b, Figure 7a and Figure 7b show a comparison of the $L_2,$ $L_{\infty },$ and $RMS$ errors computed using the LMM based on MQRBF coupled with ITM for LT inversion and the LMM based on GSRBF combined with STM for LT inversion for $\tau ,\mu \in [0.1,1]$ and $\nu \in [1.1,2]$, respectively. Overall, the findings indicate that the performance of the LMM based on MQRBF coupled with either ITM or STM significantly outperforms that of the LMM based on GSRBF.

Next, we consider two dimensional problems which are solved on three different domains, as presented in Figure 8a–c.

6.2. Example 2

In the second example, we analyze Equation (4) with exact analytical $u(x,y,\tau )=({\tau }^3+1)sin\left(\pi x\right)sin\left(\pi y\right).$ The numerical solution is investigated on two domains, namely, ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2.$ The $L_2$, $L_{\infty },$ and $RMS$ errors obtained using the LMM coupled with ITM and STM on these domains are detailed in

Table 3. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 2 with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_1.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	1225	55	12	8.16$\times {10}^{-2}$	4.74$\times {10}^{-3}$	6.66$\times {10}^{-5}$	4.88	12.822572
			14	6.32$\times {10}^{-3}$	3.11$\times {10}^{-4}$	5.16$\times {10}^{-6}$	4.88	7.879157
			16	2.28$\times {10}^{-3}$	1.55$\times {10}^{-4}$	1.86$\times {10}^{-6}$	4.88	8.687972
	625	50	18	1.67$\times {10}^{-3}$	1.61$\times {10}^{-4}$	2.68$\times {10}^{-6}$	3.42	2.528842
	676			1.67$\times {10}^{-3}$	1.53$\times {10}^{-4}$	2.47$\times {10}^{-6}$	3.56	2.819362
	1600			6.56$\times {10}^{-3}$	3.62$\times {10}^{-4}$	4.10$\times {10}^{-6}$	5.56	13.686297
[26]	961		20	6.56$\times {10}^{-6}$	4.36$\times {10}^{-7}$	1.70$\times {10}^{-7}$	—-	4.409938
[18]				1.40$\times {10}^{-3}$
STM			$M_S$
	1225	55	8	1.18$\times {10}^{+0}$	6.96$\times {10}^{-2}$	3.38$\times {10}^{-2}$	4.8	4.460539
			10	9.90$\times {10}^{-2}$	5.77$\times {10}^{-3}$	2.82$\times {10}^{-3}$	4.8	4.432449
			12	1.50$\times {10}^{-3}$	1.24$\times {10}^{-4}$	4.31$\times {10}^{-5}$	4.8	4.867652
	576	50	14	1.86$\times {10}^{-3}$	3.76$\times {10}^{-4}$	7.78$\times {10}^{-5}$	3.2	1.606402
	625			2.32$\times {10}^{-3}$	4.83$\times {10}^{-4}$	9.28$\times {10}^{-5}$	3.4	1.789516
	676			2.56$\times {10}^{-3}$	6.64$\times {10}^{-4}$	9.87$\times {10}^{-5}$	3.5	1.935365
[26]	961		10	6.99$\times {10}^{-2}$	5.97$\times {10}^{-3}$	2.33$\times {10}^{-3}$	—-	0.581491
[18]				1.40$\times {10}^{-3}$

,

Table 4. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 2 with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_2$.

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	1611	58	12	1.18$\times {10}^{-1}$	4.94$\times {10}^{-3}$	7.33$\times {10}^{-5}$	10	9.280522
			14	1.15$\times {10}^{-2}$	5.22$\times {10}^{-4}$	7.15$\times {10}^{-6}$	10	9.494589
			16	6.03$\times {10}^{-4}$	6.41$\times {10}^{-5}$	3.74$\times {10}^{-7}$	10	10.208593
	2279	45	26	2.37$\times {10}^{-3}$	1.05$\times {10}^{-4}$	1.04$\times {10}^{-6}$	10	24.286671
		55		1.72$\times {10}^{-3}$	7.90$\times {10}^{-5}$	7.58$\times {10}^{-7}$	10	24.833455
		58		1.58$\times {10}^{-3}$	7.34$\times {10}^{-5}$	6.96$\times {10}^{-7}$	10	27.383034
STM			$M_S$
	1366	54	10	1.27$\times {10}^{-1}$	5.80$\times {10}^{-3}$	3.45$\times {10}^{-3}$	5.8	4.430179
			12	1.15$\times {10}^{-3}$	1.66$\times {10}^{-4}$	3.12$\times {10}^{-5}$	5.8	4.212679
			14	3.57$\times {10}^{-3}$	2.07$\times {10}^{-4}$	9.67$\times {10}^{-5}$	5.8	6.014748
	2299	45	12	9.05$\times {10}^{-3}$	3.83$\times {10}^{-4}$	1.88$\times {10}^{-4}$	6.9	8.843242
		47		7.19$\times {10}^{-3}$	2.97$\times {10}^{-4}$	1.50$\times {10}^{-4}$	7.2	8.970401
		54		1.73$\times {10}^{-3}$	1.00$\times {10}^{-4}$	3.60$\times {10}^{-5}$	7.5	9.941594
[18]				1.40$\times {10}^{-3}$

and

Table 5. Obtained errors using the LMM based on GSRBF coupled with ITM and STM for Example 2 with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_2.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	1611	58	12	1.19$\times {10}^{-1}$	5.25$\times {10}^{-3}$	7.42$\times {10}^{-5}$	10	11.694997
			14	1.34$\times {10}^{-2}$	8.16$\times {10}^{-4}$	8.33$\times {10}^{-6}$	10	12.896010
			16	3.80$\times {10}^{-3}$	8.84$\times {10}^{-4}$	2.36$\times {10}^{-6}$	10	14.197171
	2279	45	26	3.77$\times {10}^{-2}$	6.53$\times {10}^{-3}$	1.65$\times {10}^{-5}$	10	44.586312
		55		1.28$\times {10}^{-1}$	5.80$\times {10}^{-2}$	5.61$\times {10}^{-5}$	10	43.674940
		58		1.36$\times {10}^{-1}$	3.72$\times {10}^{-2}$	5.97$\times {10}^{-5}$	10	46.742999
STM			$M_S$
	1366	54	10	1.42$\times {10}^{-1}$	2.82$\times {10}^{-2}$	3.84$\times {10}^{-3}$	9.1	5.152113
			12	5.67$\times {10}^{-2}$	2.60$\times {10}^{-2}$	1.53$\times {10}^{-3}$	9.1	5.450431
			14	5.66$\times {10}^{-2}$	2.61$\times {10}^{-2}$	1.53$\times {10}^{-3}$	9.1	5.538955
	2299	45	12	5.66$\times {10}^{-2}$	4.62$\times {10}^{-3}$	1.18$\times {10}^{-3}$	10	10.863883
		47		3.72$\times {10}^{-3}$	3.28$\times {10}^{-3}$	7.76$\times {10}^{-4}$	10	11.165862
		54		6.05$\times {10}^{-2}$	9.86$\times {10}^{-3}$	1.26$\times {10}^{-3}$	10	12.157722
[18]				1.40$\times {10}^{-3}$

. Figure 9a and Figure 9b present the computed approximate solutions on the domains ${\mathsf{\Omega }}_1$ and ${\mathsf{\Omega }}_2,$ respectively. The contour plot of $L_{Abs}$ using the LMM based on MQRBF coupled with ITM with $N_{glob}=2500,n_{loc}=60,$ and $M_T=26,$ is depicted in Figure 10a. The corresponding result using the LMM based on GSRBF coupled with STM with $N_{glob}=3600,n_{loc}=35,$ and $M_S=12,$ is shown in Figure 10b.

Figure 11a compares the $L_2$, $L_{\infty },$ and $RMS$ errors using the LMM combined with ITM for different values of $M_T$ with $N_{glob}=2366$ and $n_{loc}=58$ on ${\mathsf{\Omega }}_2.$ Similarly, Figure 11b presents the same error matrices using the LMM combined with STM for varying $M_S$ with $N_{glob}=2500$ and $n_{loc}=20$ on ${\mathsf{\Omega }}_1.$

Figure 12a, Figure 12b, Figure 13a, Figure 13b, Figure 14a and Figure 14b provide a detailed comparison of the $L_2,$ $L_{\infty },$ and $RMS$ errors on ${\mathsf{\Omega }}_2$ when using the LMM based on MQRBF combined with STM and using the LMM based on GSRBF combined with ITM. These results demonstrate the stability of the proposed method across varying values of $\mu ,\alpha ,\nu .$

The surface plot of $L_{Abs}$ using the LMM based on GSRBF coupled with STM on domain ${\mathsf{\Omega }}_1$ is shown in Figure 15a, with $N_{glob}=2500,n_{loc}=50,M_S=14.$ The corresponding surface plot of $L_{Abs}$ on ${\mathsf{\Omega }}_2$ using the LMM based on MQRBF coupled with ITM is depicted in Figure 15b, with $N_{glob}=2279,n_{loc}=58,M_T=26.$ From the obtained results, it is evident that while the LMM based on MQRBF and GSRBF achieves excellent accuracy, the LMM based on MQRBF coupled with ITM or STM shows superior performance in terms of both accuracy and stability.

Additionally, we compared our results with two existing methods. The first comparison was with a finite difference method [18], where our results outperformed the previously reported results. The second comparison involved a Chebyshev method combined with ITM and STM [26]. Our results aligned closely with those obtained using the previous method, although it is worth noting that the results from the reference paper using ITM were slightly better than ours. Nevertheless, few studies have applied spectral collocation methods to irregular geometries, whereas our local RBF method demonstrates strong capability to efficiently solve problems defined on such complex geometries.

6.3. Example 3

In the third example, we analyze Equation (4) with the analytical solution $u(x,y,\tau )=({\tau }^3+1)(1-x^2-y^2).$ The numerical solution is investigated over two domains, namely, ${\mathsf{\Omega }}_2$ and ${\mathsf{\Omega }}_3.$ For this example, the $L_2$, $L_{\infty },$ and $RMS$ errors obtained using the LMM coupled with ITM and STM on both domains are presented in

Table 6. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 3, with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_2.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	1631	75	16	1.08$\times {10}^{-3}$	9.30$\times {10}^{-5}$	2.69$\times {10}^{-5}$	7.2	13.962074
			18	2.61$\times {10}^{-4}$	5.01$\times {10}^{-5}$	6.47$\times {10}^{-6}$	7.2	13.389609
			20	2.44$\times {10}^{-4}$	4.62$\times {10}^{-5}$	6.06$\times {10}^{-6}$	7.2	14.136106
	1631	58	26	4.97$\times {10}^{-4}$	1.03$\times {10}^{-4}$	1.23$\times {10}^{-5}$	6.5	14.534135
		65		6.83$\times {10}^{-4}$	7.93$\times {10}^{-5}$	1.69$\times {10}^{-5}$	6.7	15.105693
		70		6.35$\times {10}^{-4}$	7.23$\times {10}^{-5}$	1.57$\times {10}^{-5}$	7.0	15.680614
STM			$M_S$
	1356	35	12	3.43$\times {10}^{-3}$	1.60$\times {10}^{-4}$	9.34$\times {10}^{-5}$	4.6	2.603911
			16	1.17$\times {10}^{-3}$	1.64$\times {10}^{-4}$	3.18$\times {10}^{-5}$	4.6	2.905730
			18	8.02$\times {10}^{-4}$	2.50$\times {10}^{-4}$	2.18$\times {10}^{-5}$	4.6	3.085991
	1356	43	18	1.60$\times {10}^{-3}$	4.89$\times {10}^{-4}$	4.36$\times {10}^{-5}$	5	3.541962
		45		1.35$\times {10}^{-3}$	3.61$\times {10}^{-4}$	3.67$\times {10}^{-5}$	5	3.549959
		48		8.63$\times {10}^{-4}$	1.61$\times {10}^{-4}$	2.34$\times {10}^{-5}$	5	3.760876
[18]				1.40$\times {10}^{-3}$

,

Table 7. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 3, with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_3.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	5014	55	12	1.77$\times {10}^{-1}$	4.77$\times {10}^{-3}$	2.50$\times {10}^{-3}$	10	187.506487
			14	1.41$\times {10}^{-2}$	5.87$\times {10}^{-4}$	2.00$\times {10}^{-4}$	10	206.602699
			16	5.98$\times {10}^{-3}$	1.86$\times {10}^{-4}$	8.44$\times {10}^{-5}$	10	235.925125
	3804	70	18	8.40$\times {10}^{-3}$	2.62$\times {10}^{-4}$	1.36$\times {10}^{-4}$	9.3	132.534736
	4357			2.09$\times {10}^{-3}$	1.01$\times {10}^{-4}$	3.16$\times {10}^{-5}$	10	185.875164
	4452			6.33$\times {10}^{-3}$	3.01$\times {10}^{-4}$	9.49$\times {10}^{-5}$	10	248.616861
STM			$M_S$
	5064	55	8	2.58$\times {10}^{+0}$	6.68$\times {10}^{-2}$	3.63$\times {10}^{-2}$	7.2	59.997980
			10	2.16$\times {10}^{-1}$	5.77$\times {10}^{-3}$	3.03$\times {10}^{-3}$	7.2	67.779044
			12	4.64$\times {10}^{-3}$	2.45$\times {10}^{-4}$	6.52$\times {10}^{-5}$	7.2	71.659161
	3677	70	14	2.12$\times {10}^{-2}$	5.18$\times {10}^{-3}$	3.49$\times {10}^{-4}$	6.3	42.595716
	3814		14	9.74$\times {10}^{-3}$	2.97$\times {10}^{-4}$	1.57$\times {10}^{-4}$	6.4	45.219558
	4028			8.04$\times {10}^{-3}$	2.30$\times {10}^{-4}$	1.26$\times {10}^{-4}$	6.5	51.109118
[18]				1.40$\times {10}^{-3}$

and

Table 8. Obtained errors using the LMM based on GSRBF coupled with ITM and STM for Example 3, with $\nu =1.6,\mu =0.7,\alpha =0.8$ on ${\mathsf{\Omega }}_3.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	5014	55	12	2.69$\times {10}^{-1}$	3.72$\times {10}^{-2}$	3.80$\times {10}^{-3}$	1.3	175.935646
			14	2.52$\times {10}^{-1}$	3.46$\times {10}^{-2}$	3.56$\times {10}^{-3}$	1.3	207.152231
			16	2.56$\times {10}^{-1}$	3.43$\times {10}^{-2}$	3.62$\times {10}^{-3}$	1.3	232.079818
	3804	70	18	4.66$\times {10}^{-3}$	7.61$\times {10}^{-2}$	7.57$\times {10}^{-3}$	1.34	128.708713
	4357			4.95$\times {10}^{+0}$	2.07$\times {10}^{-1}$	7.51$\times {10}^{-2}$	1.40	182.794493
	4452			7.01$\times {10}^{+1}$	3.21$\times {10}^{+1}$	1.05$\times {10}^{+0}$	1.41	197.282165
STM			$M_S$
	5064	55	8	2.26$\times {10}^{+0}$	7.00$\times {10}^{-2}$	3.18$\times {10}^{-2}$	1.29	51.834455
			10	6.09$\times {10}^{-1}$	2.70$\times {10}^{-2}$	8.56$\times {10}^{-3}$	1.29	57.752426
			12	7.10$\times {10}^{-1}$	2.98$\times {10}^{-2}$	9.98$\times {10}^{-3}$	1.29	65.394584
	3677	70	14	6.90$\times {10}^{+0}$	1.73$\times {10}^{+0}$	1.13$\times {10}^{-1}$	1.29	40.368546
	3814			3.42$\times {10}^{-1}$	3.60$\times {10}^{-2}$	5.55$\times {10}^{-3}$	1.35	41.720485
	4028			7.77$\times {10}^{-1}$	3.66$\times {10}^{-2}$	1.22$\times {10}^{-2}$	1.38	50.557339
[18]				1.40$\times {10}^{-3}$

. Figure 16a and Figure 16b respectively show the computed approximate solutions on domains ${\mathsf{\Omega }}_2$ and ${\mathsf{\Omega }}_3$.

The contour plot of $L_{Abs}$ using the LMM based on MQRBF combined with ITM for LT inversion with $N_{glob}=1600,n_{loc}=20,$ and $M_T=26$ is depicted in Figure 17a. The corresponding contour plot using the LLM based on GSRBF coupled with STM with $N_{glob}=1600,n_{loc}=20,$ and $M_S=14$ is shown in Figure 17b.

In Figure 18a, the $L_2$, $L_{\infty },$ and $RMS$ errors using the LMM coupled with ITM are compared for different values of $M_T,$ with $N_{glob}=2309$ and $n_{loc}=55$ on ${\mathsf{\Omega }}_2$. Similarly, Figure 18b presents a comparison of the same errors using the LMM coupled with STM for different values of $M_S,$ with $N_{glob}=2908$ and $n_{loc}=70$ on ${\mathsf{\Omega }}_3$.

Figure 19a, Figure 19b, Figure 20a, Figure 20b, Figure 21a and Figure 21b provide a detailed comparison of the $L_2,$ $L_{\infty },$ and $RMS$ errors on ${\mathsf{\Omega }}_3$ using the LMM based on GSRBF combined with STM and using the LMM based on MQRBF combined with ITM. These results demonstrate the stability of the method across varying values of $\mu ,\alpha ,\nu .$

The surface plot of $L_{Abs}$ using the LMM based on MQRBF combined with ITM on domain ${\mathsf{\Omega }}_2$ is shown in Figure 22a, with $N_{glob}=1631,n_{loc}=65$, and $M_T=26.$ The corresponding surface plot using the LMM based on GSRBF combined with STM on domain ${\mathsf{\Omega }}_3$ is displayed in Figure 22b, with $N_{glob}=2898,n_{loc}=40$, and $M_S=16$. The results show that both methods produce stable results, with acceptable accuracy even for a large number of nodes; however, the combination of the LMM based on MQRBF with ITM consistently outperforms the LMM based on GSRBF with STM in terms of both stability and accuracy.

6.4. Example 4

In the fourth example, we examine Equation (4) with the exact solution $u(x,y,z,\tau )={\tau }^2{e}^{(x+y+z)}.$ A uniform distribution of nodes over the domain ${[0,1]}^3$ is shown in Figure 23a. The approximate solution for this problem is depicted in Figure 23b. The $L_2$, $L_{\infty },$ and $RMS$ errors using the ITM and STM are reported in

Table 9. Obtained errors using the LMM based on MQRBF coupled with ITM and STM for Example 4, with $\nu =1.6,\mu =0.7,\alpha =0.8.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	2744	70	12	1.01$\times {10}^{-1}$	6.37$\times {10}^{-3}$	1.93$\times {10}^{-3}$	1.5	41.022716
			14	1.78$\times {10}^{-2}$	2.35$\times {10}^{-3}$	3.40$\times {10}^{-4}$	1.5	42.909508
			16	1.51$\times {10}^{-2}$	2.37$\times {10}^{-3}$	2.88$\times {10}^{-4}$	1.5	44.699805
	3375	80	18	4.53$\times {10}^{-3}$	4.57$\times {10}^{-4}$	7.81$\times {10}^{-5}$	1.9	69.752230
	4096			5.42$\times {10}^{-3}$	5.12$\times {10}^{-4}$	8.47$\times {10}^{-5}$	2.0	103.137213
	4913			5.46$\times {10}^{-3}$	5.11$\times {10}^{-4}$	7.79$\times {10}^{-5}$	2.1	142.095793
STM			$M_S$
	1000	60	10	2.39$\times {10}^{-2}$	6.74$\times {10}^{-3}$	7.55$\times {10}^{-4}$	0.9	6.421360
			12	2.32$\times {10}^{-2}$	6.60$\times {10}^{-3}$	7.34$\times {10}^{-4}$	0.9	6.577626
			14	2.22$\times {10}^{-2}$	6.06$\times {10}^{-3}$	7.02$\times {10}^{-4}$	0.9	6.846112
	1331	70	16	4.04$\times {10}^{-3}$	6.28$\times {10}^{-4}$	1.10$\times {10}^{-4}$	1.2	12.352297
	1728			4.79$\times {10}^{-3}$	7.96$\times {10}^{-4}$	1.15$\times {10}^{-4}$	1.3	17.042883
	2197			7.19$\times {10}^{-3}$	1.24$\times {10}^{-3}$	1.53$\times {10}^{-4}$	1.4	25.053652

and

Table 10. Obtained errors using the LMM based on GSRBF coupled with ITM and STM for Example 4, with $\nu =1.6,\mu =0.7,\alpha =0.8.$

ITM	${\mathit{N}}_{\mathit{glob}}$	${\mathit{n}}_{\mathit{loc}}$	${\mathit{M}}_{\mathit{T}}$	${\mathit{L}}_2$	${\mathit{L}}_{\mathit{\infty }}$	$\mathit{RMS}$	$\mathit{\epsilon }$	C.Time(s)
	2744	70	12	1.71$\times {10}^{-1}$	2.90$\times {10}^{-2}$	3.27$\times {10}^{-3}$	2.1	43.519792
			14	2.02$\times {10}^{-1}$	2.58$\times {10}^{-2}$	3.85$\times {10}^{-3}$	2.1	42.735453
			16	2.06$\times {10}^{-1}$	2.55$\times {10}^{-2}$	3.94$\times {10}^{-3}$	2.1	44.671895
	3375	80	18	1.74$\times {10}^{-1}$	1.24$\times {10}^{-2}$	3.01$\times {10}^{-3}$	2.6	69.157279
	4096			2.12$\times {10}^{-1}$	1.58$\times {10}^{-2}$	3.41$\times {10}^{-3}$	2.8	100.112121
	4913			1.78$\times {10}^{-1}$	1.10$\times {10}^{-2}$	2.54$\times {10}^{-3}$	2.9	146.783429
STM			$M_S$
	1000	60	10	1.84$\times {10}^{-1}$	7.05$\times {10}^{-2}$	5.82$\times {10}^{-3}$	1.2	7.577898
			12	1.84$\times {10}^{-1}$	7.07$\times {10}^{-2}$	5.84$\times {10}^{-3}$	1.2	6.926979
			14	1.86$\times {10}^{-1}$	7.13$\times {10}^{-2}$	5.90$\times {10}^{-3}$	1.2	7.190806
	1331	70	16	1.15$\times {10}^{-1}$	1.12$\times {10}^{-2}$	3.17$\times {10}^{-3}$	1.6	12.808148
	1728			1.33$\times {10}^{-1}$	1.17$\times {10}^{-2}$	3.21$\times {10}^{-3}$	1.8	17.709936
	2197			1.37$\times {10}^{-1}$	1.17$\times {10}^{-2}$	2.93$\times {10}^{-3}$	1.9	26.220721

, respectively.

In Figure 24a, the $L_2$, $L_{\infty },$ and $RMS$ errors are compared using the LMM coupled with ITM for different values of $M_T,$ with $N_{glob}=3375$ and $n_{loc}=80.$ Similarly, Figure 24b compares the same errors using the LMM coupled with STM for different of $M_S,$ with $N_{glob}=2197$ and $n_{loc}=70$.

A slice plot of $L_{Abs}$ using the LMM based on GSRBF coupled with ITM is shown in Figure 25a, with $N_{glob}=4913,n_{loc}=80$, and $M_T=28.$ Figure 25b depicts the corresponding slice plot of $L_{Abs}$ using the LMM based on MQRBF coupled with STM, with $N_{glob}=2197,n_{loc}=70$, and $M_S=18$.

The contour slice plot of $L_{Abs}$ using the LMM based on GSRBF coupled with ITM is shown in Figure 26a, with $N_{glob}=5832,n_{loc}=80,$ and $M_T=30.$ The corresponding contour slice plot of $L_{Abs}$ using the LMM based on MQRBF coupled with STM is shown in Figure 26b, with $N_{glob}=4096,n_{loc}=70,$ and $M_S=18.$

The results demonstrate that both methods provide accurate approximations of the exact solution for three-dimensional problems. The combination of MQRBF-based and GSRBF-based LMM with ITM and STM consistently yields stable results across various parameter settings, showcasing the robustness of the proposed method for complex high-dimensional problems.

7. Conclusions

In this article, we have developed an effective hybrid numerical technique coupling the LT and LMM to address fractional-order mixed multi-term subdiffusion and wave-diffusion equations. We employ the LT technique to effectively handle the time-fractional derivative, while the LMM allows for discretization of the spatial derivatives, resulting in a fully discrete system. We solve the discrete system at each point within the complex domain and obtain the solution to the considered problem through LT inversion. Two widely used algorithms, namely, the Stehfest algorithm and improved Talbot algorithm, are utilized for approximating the inverse LT. Our comparative analysis reveals that ITM performs better than STM. Furthermore, we explored the robustness of the method using MQRBF and GSRBF, demonstrating that MQRBF yields superior results. The convergence of each numerical Laplace inversion algorithm and the stability of thhe fully discretized sparse system are thoroughly discussed. Numerical experiments conducted on different domains (including square, star, and circular with uniform node distributions) were carried out to validate the effectiveness of the proposed numerical scheme. Our results indicate that the proposed method effectively solves problems on irregular domains while mitigating issues related to temporal instability The proposed numerical scheme proves effective in solving problems on irregular domains while avoiding issues related to time instability. The precision and stability of the suggested scheme have been rigorously validated through numerical experiments, showing good agreement with the corresponding analytical solutions.