An Efficient Numerical Solution of a Multi-Dimensional Two-Term Fractional Order PDE via a Hybrid Methodology: The Caputo–Lucas–Fibonacci Approach with Strang Splitting

Abstract

The utilization of time-fractional PDEs in diverse fields within science and technology has attracted significant interest from researchers. This paper presents a relatively new numerical approach aimed at solving two-term time-fractional PDE models in two and three dimensions. We combined the Liouville–Caputo fractional derivative scheme with the Strang splitting algorithm for the temporal component and employed a meshless technique for spatial derivatives utilizing Lucas and Fibonacci polynomials. The rising demand for meshless methods stems from their inherent mesh-free nature and suitability for higher dimensions. Moreover, this approach demonstrates the effective approximation of solutions across both regular and irregular domains. Error norms were used to assess the accuracy of the methodology across both regular and irregular domains. A comparative analysis was conducted between the exact solution and alternative numerical methods found in the contemporary literature. The findings demonstrate that our proposed approach exhibited better performance while demanding fewer computational resources.

1. Introduction

Over the last decade, there has been considerable interest in fractional partial differential equations (PDEs), making it a vibrant research field for scientists and engineers. PDEs possess the capability to describe numerous complex phenomena across diverse fields, including biology, fluid mechanics, plasma physics, optics, acoustics, financial mathematics, climate modeling, materials science, and electromagnetics [1,2]. Fractional calculus, which deals with derivatives and integrals of non-integer order, finds broad applications in science and technology. In physics, it models complex systems with non-local behavior, such as anomalous diffusion in fluid dynamics and memory effects in viscoelastic materials [3]. Engineering benefits from fractional calculus are found in signal processing, control theory, and optimization, offering robustness and efficiency in dynamic systems. Telecommunications utilize it for modeling wireless communication channels and improving channel estimation algorithms. Fractional calculus enables advancements in understanding and optimizing real-world processes and systems across various domains [4]. The fractional spatio-temporal PDE models are crucial in describing complex phenomena [5], whereas the two-scale fractal PDE models are essential for simulating intricate structures in materials science and physics, offering insights into multiscale behaviors and enhancing our understanding of fractal geometries in natural and engineered systems [6]. In particular, the Sobolev model equation is integral across various scientific fields, providing a robust framework for investigating phenomena with spatial variations. Its relevance extends across disciplines, encompassing physics, engineering, biology, and finance. In physics, Sobolev PDEs are commonly utilized to model phenomena such as heat conduction, fluid dynamics, and quantum mechanics [7,8]. Engineers employ them to analyze structural mechanics, electromagnetism, and signal-processing problems [7]. Furthermore, Sobolev PDEs find applications in image processing, medical imaging, and geophysical exploration. Their adaptability lies in their capacity to handle irregular domains and boundary conditions while providing solutions that demonstrate smoothness properties.

Numerous studies explored the applications and theoretical foundations of Sobolev PDEs, underscoring their significance in both theoretical and applied contexts. However, numerous researchers faced challenges in deriving and formulating various complex phenomena within nonlinear PDEs with integer orders [9]. In response, fractional calculus is regarded as a viable solution to this issue, as it incorporates a nonlocal property that is absent in nonlinear PDEs with integer orders [10]. Observations indicate that multi-term time-fractional PDEs were proposed to enhance the modeling accuracy in depicting anomalous diffusion processes, capturing various types of viscoelastic damping, accurately representing power-law frequency dependence, and simulating the flow of a fractional Maxwell fluid [11]. This study focused on the two-term time-fractional Sobolev equation in two and three dimensions, defined as follows:

$$\begin{array}{c}\frac{{\partial }^{{\xi }_1}\mathcal{U}}{\partial t^{{\xi }_1}}+\frac{{\partial }^{{\xi }_2}\mathcal{U}}{\partial t^{{\xi }_2}}-\beta \frac{\partial }{\partial t}\left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}\right)-\delta \left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}\right)+\eta \left(\mathcal{U}\left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}\right)\right)\\ +\eta {\left(\frac{\partial \mathcal{U}}{\partial r}+\frac{\partial \mathcal{U}}{\partial s}\right)}^2+\gamma \mathcal{U}=F(r,s,t),r,s\in \mathrm{\Omega }\subset {\mathbb{R}}^2,0<{\xi }_2\le {\xi }_1\le 1,t>0,\end{array}$$

(1)

and

$$\begin{array}{c}\frac{{\partial }^{{\xi }_1}\mathcal{U}}{\partial t^{{\xi }_1}}+\frac{{\partial }^{{\xi }_2}\mathcal{U}}{\partial t^{{\xi }_2}}-\beta \frac{\partial }{\partial t}\left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}+\frac{{\partial }^2\mathcal{U}}{\partial z^2}\right)-\delta \left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}+\frac{{\partial }^2\mathcal{U}}{\partial z^2}\right)\hfill \\ +\eta \left(\mathcal{U}\left(\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}+\frac{{\partial }^2\mathcal{U}}{\partial z^2}\right)\right)+\eta {\left(\frac{\partial \mathcal{U}}{\partial r}+\frac{\partial \mathcal{U}}{\partial s}+\frac{\partial \mathcal{U}}{\partial z}\right)}^2+\gamma \mathcal{U}\hfill \\ =F(r,s,z,t),r,s,z\in \mathrm{\Omega }\subset {\mathbb{R}}^3,0<{\xi }_2\le {\xi }_1\le 1,t>0,\hfill \end{array}$$

(2)

with the following initial and boundary conditions:

$$\mathcal{U}(\overline{\mathit{r}},0)={\mathcal{U}}_0\left(\overline{\mathit{r}}\right),$$

(3)

$$\mathcal{U}(\overline{\mathit{r}},t)=f_1(\overline{\mathit{r}},t),\overline{\mathit{r}}\in \partial \mathrm{\Omega }\subset {\mathbb{R}}^d,d=2,3.$$

(4)

where the constants $\beta$, $\delta$, $\eta$, and $\gamma$ have known values. Moreover, $\frac{{\partial }^{{\xi }_1}}{\partial t^{{\xi }_1}}$ and $\frac{{\partial }^{{\xi }_2}}{\partial t^{{\xi }_2}}$ denote the Caputo derivatives in accordance with Caputo [12] and applied to $\mathcal{U}(\overline{\mathit{r}},t)$, where $0<{\xi }_2\le {\xi }_1\le 1$.

Several works delved into the existence and uniqueness of solutions for Sobolev equations, as documented in [13]. Ewing [14] applied the finite difference method to tackle 2D Sobolev equations. The study [15] used a computationally appealing and precise local meshless technique to provide numerical solutions for three-dimensional two- and three-term time-fractional PDE models. Various test problems were employed to assess the accuracy and reliability of the proposed approach. The paper by Luo et al. [16] introduces a reduced-order extrapolated finite difference iterative scheme for 2D Sobolev equations, employing proper orthogonal decomposition to construct a reduced-order solution space model, thus reducing the computational costs. Numerical experiments validated the scheme’s efficiency and accuracy, indicating its potential for efficient Sobolev equation solving. An accurate and efficient local meshless technique is used in the article [17] to explore numerical solutions for two-term time-fractional Sobolev models. The method approximates the solution on a uniform or scattered set of nodes, yielding sparse and well-conditioned coefficient matrices. Another paper by Luo et al. [18] introduces a novel method that employs reduced-order techniques and extrapolation within a Crank–Nicolson finite volume framework for 2D Sobolev equations. It aims to improve the computational efficiency and accuracy, showcasing promising outcomes via numerical simulations. The paper by Li et al. [19] introduces an extended mixed finite element approach for solving 2D Sobolev equations, integrating a novel formulation with stable discretization for accurate solutions and computational efficiency. Numerical experiments highlighted its effectiveness and superiority over existing methods, especially in maintaining stability and convergence properties. Heydari et al. [20] introduced a novel method that employs orthonormal Bernoulli polynomials to solve distributed-order time-fractional 2D Sobolev equations, which demonstrated accuracy through four test problems. Gao et al. [21] used the local discontinuous Galerkin finite element method to solve a particular class of 2D Sobolev equations, while another work [22] applied the weak Galerkin finite element method, providing an error estimate. Abu et al. [23] presented a method for solving time- and space-fractional Sobolev equations with Caputo fractional derivatives in n-dimensional space, utilizing the reproducing kernel Hilbert space method, which is particularly effective for Caputo class derivatives. A meshless RBFs technique for solving 2D time-fractional Sobolev equations was presented by Hussain et al. [24]. This method uses a finite difference formula for time-fractional derivative approximation and RBFs for spatial operator approximation. Using a Liouville–Caputo derivative technique for the time derivative, Ahmad et al. [25] suggested an efficient meshless method for estimating the numerical solution of 3D time-fractional Sobolev equations. Furthermore, Zhang et al. [26] used a characteristic splitting mixed finite element approach to solve convection-dominated Sobolev equations.

In response to its extensive applications, researchers explored a multitude of numerical and analytical methodologies [27,28] for addressing complex PDEs. These methods encompass various techniques, such as the alternating direction implicit method [29], the homotopy perturbation method [30], the Laplace transform technique [31], the variational approach [32,33], finite element methods (FEMs) [34], finite difference methods (FDMs) [16,35], gradient descent iterative method [36], spectral methods [37], the exp-function method [38], the alternating direction method [39], and meshless methods [40]. Notably, the straightforward nature of the FDM, FEM, and meshless techniques is noteworthy. Recently, hybrid methodologies have emerged to enhance the efficiency and accuracy of numerical solutions, including meshless methods, the method of lines utilizing Fibonacci polynomials, and combinations of the FDM and FEM.

This investigation focused on computing numerical solutions for the suggested model using a hybrid methodology based on the Caputo derivative. This approach integrates Fibonacci polynomials with the established Caputo derivative concept, exploiting the relationship between Fibonacci and Lucas polynomials. This integration provides a significant advantage in the straightforward implementation of higher-order derivatives. Moreover, the proposed approach reduces computational costs by enhancing the accuracy, even with a limited number of collocation sites. Remarkably, these polynomials find diverse real-world applications in the realm of differential equations.

Addressing boundary value problems accurately involves exploring the interplay between Chebyshev and Lucas polynomials, as demonstrated in previous studies [41,42]. For instance, the Lucas sequence was utilized to approximate integro-differential equations [43], while Lucas polynomials were applied to solve higher-order differential equations [44]. Additionally, the efficacy of a Fibonacci polynomial methodology in resolving Volterra–Fredholm integral differential equations was demonstrated [45]. Furthermore, a hybrid Taylor–Lucas polynomial approach was introduced for addressing delay difference equations [46]. Notably, novel methodologies for solving time-dependent PDEs were proposed by combining hybrid Fibonacci and Lucas polynomial schemes [47,48]. Moreover, researchers employed finite differences and Lucas polynomials to achieve effective numerical solutions for various PDE models [49,50].

Motivation

The primary objective of this research was to introduce a relatively new numerical approach designed to solve two-term time-fractional PDE models in both two and three dimensions. This method combines the Liouville–Caputo fractional derivative scheme with the Strang splitting algorithm for the temporal component and utilizes a meshless technique for spatial derivatives, incorporating Lucas and Fibonacci polynomials. Below are some of the highlighted key features of the proposed study:

• The complex characteristics of fractional nonlinear PDEs make calculating analytical solutions challenging, driving ongoing research efforts to develop accurate and efficient numerical methodologies, with the two-term fractional order Sobolev model equation in both two and three dimensions holding significant importance across multiple scientific domains.
• This study aimed to introduce an efficient numerical framework specifically designed for solving PDEs with temporal fractions.
• The proposed methodology adopts a hybrid approach, integrating Fibonacci and Lucas polynomials with finite difference techniques, while also addressing the temporal direction through the utilization of the Liouville–Caputo fractional derivative in conjunction with a splitting mechanism.
• Lucas and Fibonacci polynomials, unlike orthogonal counterparts, like Chebyshev polynomials, are non-orthogonal, eliminating the need for interval transformations. Additionally, they facilitate the straightforward approximation of higher-order derivatives for unknown functions.
• Furthermore, the approach is characterized by its simplicity and ability to enhance the accuracy, even in scenarios involving fewer nodal points, with the aim to provide a robust and effective numerical solution to the intricate challenges posed by nonlinear PDEs.

This paper is organized as follows: Section 2 provides an overview of the fundamental terms and concepts. The proposed technique for the underlying model equations is discussed in Section 3, while theoretical results regarding stability and error analysis are presented in Section 4. Section 5 utilizes numerical experiments to validate the method’s efficacy, and finally, Section 6 summarizes the outcomes and presents concluding remarks to finalize the work.

2. Basic Concepts in Fractional Calculus and Polynomial Theory

In fractional calculus, fractional derivatives are crucial. Fractional calculus finds applications in various fields, including physics, engineering, biology, and finance, providing a powerful framework for modeling complex phenomena characterized by memory, non-locality, and fractal behavior. The fundamental definition of fractional calculus involves the generalization of differentiation and integration to non-integer orders. Here are some basic definitions of frequently used fractional derivatives.

Definition 1. 

The Riemann–Liouville derivative [51,52] is a mathematical concept used to generalize the notion of a derivative to functions that are not necessarily differentiable in the traditional sense. It is commonly applied in the field of fractional calculus to describe the behavior of functions with fractional orders of differentiation.

$$\frac{{\partial }^{\xi }\mathcal{U}(\overline{\mathbf{r}},t)}{\partial t^{\xi }}=\frac{1}{\mathrm{\Gamma }\left(1-\xi \right)}\frac{d}{dt}\underset{t}{\overset{T}{\int }}\frac{(\mathcal{U}(\overline{\mathbf{r}},\vartheta )-\mathcal{U}(\overline{\mathbf{r}},T))}{{(\vartheta -t)}^{\xi }}d\vartheta ,0<\xi <1.$$

(5)

Definition 2. 

Caputo’s fractional derivative [12] is a widely used mathematical tool in fractional calculus, providing a powerful approach for describing and analyzing non-local and memory-dependent phenomena in various scientific and engineering fields.

$$\frac{{\partial }^{\xi }\mathcal{U}(\overline{\mathbf{r}},t)}{\partial t^{\xi }}=\frac{1}{\mathrm{\Gamma }\left(1-\xi \right)}\underset{0}{\overset{t}{\int }}\frac{\partial \mathcal{U}(\overline{\mathbf{r}},\zeta )}{\partial \zeta }{\left(t-\zeta \right)}^{-\xi }d\zeta ,0<\xi <1.$$

(6)

Definition 3. 

The Atangana and Baleanu fractional derivative [53], proposed by Atangana and Baleanu, is a non-local operator that extends the concept of fractional calculus. It provides a novel approach for modeling complex systems that exhibit memory and non-local behavior.

$$\underset{a}{A}BC\frac{{\partial }^{\xi }\mathcal{U}(\overline{\mathbf{r}},t)}{\partial t^{\xi }}=\frac{B\left(\xi \right)}{1-\xi }\underset{a}{\overset{t}{\int }}{\mathcal{U}}^{\prime }\left(\overline{\mathbf{r}}\right)E_{\xi }\left(-\frac{\xi {(t-\overline{\mathbf{r}})}^{\xi }}{1-\xi }\right)d\overline{\mathbf{r}},0<\xi <1.$$

(7)

Definition 4. 

He’s fractional derivative [54]:

$$\frac{{\partial }^{\xi }\mathcal{U}(\overline{\mathbf{r}},t)}{\partial t^{\xi }}=\frac{1}{\mathrm{\Gamma }\left(1-\xi \right)}\frac{d}{d\overline{\mathbf{r}}}\underset{t_0}{\overset{t}{\int }}{(t-\zeta )}^{-\xi }[{\mathcal{U}}_0\left(\zeta \right)-\mathcal{U}\left(\zeta \right)],0<\xi <1.$$

(8)

2.1. Fibonacci and Lucas Polynomial Theory

Below are explanations outlining the definitions, practical implications, and utilization of Lucas and Fibonacci polynomials, encompassing tasks such as approximating unfamiliar functions and their derivatives.

• Fibonacci polynomials [49]:

The Fibonacci polynomial, which is an extension of Fibonacci numbers, is defined through a three-term recurrence relation:

$$F_k\left(r\right)=k{F}_{k-1}\left(r\right)+F_{k-2}\left(r\right),k\ge 2,$$

(9)

where the sequence begins with $F_0\left(r\right)=0$ and $F_1\left(r\right)=1$. When $r=1$, Equation (9) generates the familiar sequence of Fibonacci numbers.

• Lucas polynomials [49]:

One can define the Lucas polynomials using the three-term recurrence relation:

$$L_k\left(r\right)=k{L}_{k-1}\left(r\right)+L_{k-2}\left(r\right),k\ge 2,$$

(10)

with $L_0\left(r\right)=2$ and $L_1\left(r\right)=r$ as initial values. This allows Equation (10) to produce a series of Lucas numbers with $r=1$.

Lemma 1 

([49]). The kth Fibonacci polynomial $F_k\left(r\right)$ can be used to define the mth-order derivative of the kth Lucas polynomial $L_k\left(r\right)$:

$$L_k^{\left(m\right)}\left(r\right)=k{F}_k\left(r\right)D^{m-1},D^{m-1}=\underset{(m-1)\mathrm{time}}{\underbrace{D\times D\times D\cdots D}}$$

(11)

where an $(M+1)\times (M+1)$ matrix is represented by D, which can be described as

$$D=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ 0\\ ⋮& & d\\ 0\end{array}\right],$$

where d is computed according to [49]:

$$d_{ij}=\left\{\begin{array}{cc}isin\frac{(j-i)\pi }{2},\hfill & \mathit{if}j>i,\hfill \\ 0,\hfill & \mathit{otherwise}.\hfill \end{array}\right.$$

2.2. Function Approximation

Suppose $\mathcal{U}\in L^2\left(\mathbb{R}\right)$, and let $\mathcal{U}\left(\overline{\mathit{r}}\right)$ be continuous. Then, through the linear combination of the kth Lucas polynomials, it is possible to represent $\mathcal{U}$ as

$$\mathcal{U}\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_k{L}_k\left(\overline{\mathit{r}}\right),$$

(12)

where the Lucas polynomials are denoted as $L_k\left(\overline{\mathit{r}}\right)$, with the corresponding unknown coefficients represented by ${\mathrm{\Lambda }}_k$.

Likewise, $\mathcal{U}\left(\overline{\mathit{r}}\right)$ can also be expanded under the same conditions by using a linear combination of the kth Fibonacci polynomials, as demonstrated below:

$$\mathcal{U}\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_k{F}_k\left(\overline{\mathit{r}}\right).$$

In this context, $F_k\left(\overline{\mathit{r}}\right)$ denotes the Fibonacci polynomials, while ${\mathrm{\Lambda }}_k$ stands for the coefficients yet to be determined.

The Lucas polynomial series can be used to expand the function $\mathcal{U}\left(\overline{\mathit{r}}\right)$ and find the first-order derivative:

$${\mathcal{U}}^{\prime }\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_k{L}_k^{\prime }\left(\overline{\mathit{r}}\right),$$

(13)

and the function $\mathcal{U}\left(\overline{\mathit{r}}\right)$’s associated $mth$-order derivative is given as

$${\mathcal{U}}^m\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_k{L}_k^{\left(m\right)}\left(\overline{\mathit{r}}\right),$$

(14)

where

$${\mathcal{U}}^m\left(\overline{\mathit{r}}\right)=\frac{{\mathrm{d}}^m\mathcal{U}\left(\overline{\mathit{r}}\right)}{\mathrm{d}{\overline{\mathit{r}}}^m},L_k^{\left(m\right)}\left(\overline{\mathit{r}}\right)=\frac{{\mathrm{d}}^m{L}_k\left(\overline{\mathit{r}}\right)}{\mathrm{d}{\overline{\mathit{r}}}^m}.$$

Equations (13) and (14) can be written by incorporating the relation (11):

$${\mathcal{U}}^{\prime }\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_{k}k{F}_k\left(\overline{\mathit{r}}\right),$$

(15)

and similarly, the mth derivative can be computed as

$${\mathcal{U}}^{\left(m\right)}\left(\overline{\mathit{r}}\right)=\sum _{k=0}^{\infty }{\mathrm{\Lambda }}_{k}k{F}_k\left(\overline{\mathit{r}}\right)D^{m-1},$$

(16)

in which D and $D^{m-1}$ are defined above.

Remark 1. 

In numerical computations, the expansion of $\mathcal{U}\left(\overline{\mathbf{r}}\right)$ and its mth derivative often involves utilizing truncated Fibonacci and Lucas polynomial series. In other words, we take

$$\mathcal{U}\left(\overline{\mathbf{r}}\right)\simeq \sum _{k=0}^M{\mathrm{\Lambda }}_k{L}_k\left(\overline{\mathbf{r}}\right),M\in \mathbb{N},$$

and

$${\mathcal{U}}^{\left(m\right)}\left(\overline{\mathbf{r}}\right)\simeq \sum _{k=0}^M{\mathrm{\Lambda }}_k{L}_k^{\left(m\right)}\left(\overline{\mathbf{r}}\right)=\sum _{k=0}^M{\mathrm{\Lambda }}_{k}k{F}_k\left(\overline{\mathbf{r}}\right)D^{m-1},M\in \mathbb{N},$$

(17)

3. Suggested Methodology

This section is devoted to formulating the suggested meshless method for approximating the 2D time-fractional Sobolev model equations given in Equations (20) and (21). To simplify the notation throughout the discussion in this section, we introduce the following:

$${\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)=\mathcal{U}(r,s,t^{n+1}),{\mathcal{U}}_{ij}^{n+1}=\mathcal{U}(r_i,s_j,t^{n+1}),$$

where $r_i$ and $s_j$ are the collocation points specified as follows, and $t^n=n\times \delta t$, with $\delta t$ representing the time step size:

$$r_i=a+i{h}_r,s_j=c+j{h}_s,(i,j = 1, 2, \cdots , M,M\in \mathbb{N}),$$

where $h_r=(b-a)/M$ and $h_s=(d-c)/M$ are the mesh step sizes in the r and s spatial directions throughout the spatial domain $\mathrm{\Xi }=[a,b]\times [c,d]\subset {\mathbb{R}}^2$.

3.1. Time Discretization

The widely recognized $L1$-formula, with $𝒪\left(\delta t^{2-{\xi }_1}\right)$ as the approximation order error, where $0<{\xi }_1\le 1$, is utilized to derive the discrete representation of the $(n+1)$th time level, as outlined in [55,56].

$$\begin{array}{cc}\hfill \frac{{\partial }^{{\xi }_1}\mathcal{U}(\overline{\mathit{r}},t^{n+1})}{\partial t^{{\xi }_1}}& =\frac{1}{\mathrm{\Gamma }(1-{\xi }_1)}{\int }_0^{t^{n+1}}\frac{\partial \mathcal{U}(\overline{\mathit{r}},\zeta )}{\partial \zeta }{(t^{n+1}-\zeta )}^{-{\xi }_1}\mathrm{d}\zeta ,\hfill \\ \\ \hfill & =\frac{1}{\mathrm{\Gamma }(1-{\xi }_1)}\sum _{j=0}^n{\int }_{j\times \delta t}^{(j+1)\times \delta t}\frac{\partial \mathcal{U}(\overline{\mathit{r}},\zeta )}{\partial \zeta }{(t^{n+1}-\zeta )}^{-{\xi }_1}\mathrm{d}\zeta ,\hfill \\ \\ \hfill & =\frac{1}{\mathrm{\Gamma }(1-{\xi }_1)}\sum _{j=0}^n\left[\frac{{\mathcal{U}}^{j+1}\left(\overline{\mathit{r}}\right)-{\mathcal{U}}^j\left(\overline{\mathit{r}}\right)}{\delta t}+O\left(\delta t\right)\right]{\int }_{j\times \delta t}^{(j+1)\times \delta t}{((j+1)\delta t-\zeta )}^{-{\xi }_1}\mathrm{d}\zeta .\hfill \end{array}$$

Following integration, it gives

$$\begin{array}{c}\hfill \frac{{\partial }^{{\xi }_1}\mathcal{U}(s,r,t^{n+1})}{\partial t^{{\xi }_1}}=\left\{\begin{array}{c}{A}_{{\xi }_1}\sum _{j=0}^n{K}_{{\xi }_1}\left(j\right)\left[{\mathcal{U}}^{n-j+1}\left(\overline{\mathit{r}}\right)-{\mathcal{U}}^{n-j}\left(\overline{\mathit{r}}\right)\right]+O\left(\delta t^{2-{\xi }_1}\right),0<{\xi }_1<1,\hfill \\ \\ \frac{{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)}{\delta t}+O\left(\delta t\right),{\xi }_1=1,\hfill \end{array}\right.\end{array}$$

(18)

where $A_{{\xi }_1}=\frac{\delta t^{-{\xi }_1}}{\mathrm{\Gamma }(2-{\xi }_1)}$ and $K_{{\xi }_1}\left(j\right)={(j+1)}^{1-{\xi }_1}-{\left(j\right)}^{1-{\xi }_1}$. Therefore, we may write the following for $0<{\xi }_1<1$ after ignoring the error term:

$$\begin{array}{c}\hfill \frac{{\partial }^{{\xi }_1}\mathcal{U}(\overline{\mathit{r}},t^{n+1})}{\partial t^{{\xi }_1}}=A_{{\xi }_1}\left[{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\right]+A_{{\xi }_1}\sum _{j=1}^n{K}_{{\xi }_1}\left(j\right)\left[{\mathcal{U}}^{n-j+1}\left(\overline{\mathit{r}}\right)-{\mathcal{U}}^{n-j}\left(\overline{\mathit{r}}\right)\right],\end{array}$$

(19)

with $K_{{\xi }_1}\left(j\right)=1$ and $j=0$. The fractional derivative of order ${\xi }_2$ can be found above.

Before applying the $\theta$-weighted rule to the suggested Equations (1)–(4), let us rewrite them in compact form as follows:

$$\begin{array}{c}\hfill \frac{{\partial }^{{\xi }_1}\mathcal{U}(\overline{\mathit{r}},t)}{\partial t^{{\xi }_1}}+\frac{{\partial }^{{\xi }_2}\mathcal{U}(\overline{\mathit{r}},t)}{\partial t^{{\xi }_2}}-\beta \frac{\partial {\nabla }^2\mathcal{U}(\overline{\mathit{r}},t)}{\partial t}-\delta {\nabla }^2\mathcal{U}(\overline{\mathit{r}},t)-\eta \nabla \left(\mathcal{U}(\overline{\mathit{r}},t)\nabla \mathcal{U}(\overline{\mathit{r}},t)\right)+\gamma \mathcal{U}(\overline{\mathit{r}},t)=f(\overline{\mathit{r}},t),\\ \hfill \overline{\mathit{r}}\in \mathrm{\Omega }\subset {\mathbb{R}}^d,d=2,3,0<{\xi }_2\le {\xi }_1\le 1,t>0,\end{array}$$

(20)

with the following conditions:

$$\begin{array}{c}\hfill \mathcal{U}(\overline{\mathit{r}},0)={\mathcal{U}}_0\left(\overline{\mathit{r}}\right),\mathcal{U}(\overline{\mathit{r}},t)=f_1(\overline{\mathit{r}},t),\overline{\mathit{r}}\in \partial \mathrm{\Omega },\end{array}$$

(21)

where the Laplacian and gradient operators are represented by the symbols ${\nabla }^2$ and ∇, respectively.

The $\theta$-weighted rule is now applied, and Equation (20) becomes

$$\begin{array}{cc}\hfill & A_{{\xi }_1}{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)+A_{{\xi }_2}{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-\beta \left(\frac{{\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-{\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)}{\delta t}\right)-\theta \delta {\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)+\theta \gamma {\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\hfill \\ \hfill & -\theta \nabla \left({\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\nabla {\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\right)=A_{{\xi }_1}{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+A_{{\xi }_2}{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+(1-\theta )\delta {\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+(1-\theta )\eta \nabla \left({\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\nabla {\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\right)\hfill \\ \hfill & -(1-\theta )\gamma {\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+f_{\theta }^{n+1}\left(\overline{\mathit{r}}\right)-Q_{\xi }^n\left(\overline{\mathit{r}}\right),\hfill \end{array}$$

(22)

where

$$f_{\theta }^{n+1}\left(\overline{\mathit{r}}\right)=\theta f^{n+1}\left(\overline{\mathit{r}}\right)+(1-\theta )f^n\left(\overline{\mathit{r}}\right),$$

and

$$Q_{\xi }^n\left(\overline{\mathit{r}}\right)=A_{{\xi }_1}\sum _{j=1}^n{K}_{{\xi }_1}\left(j\right)\left[{\mathcal{U}}^{n-j+1}(r,s)-{\mathcal{U}}^{n-j}(r,s)\right]+A_{{\xi }_2}\sum _{j=1}^n{K}_{{\xi }_2}\left(j\right)\left[{\mathcal{U}}^{n-j+1}(r,s)-{\mathcal{U}}^{n-j}(r,s)\right].$$

Since

$$\nabla \left({\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\nabla {\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\right)={\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)+{\left({\mathcal{U}}_x^{n+1}\left(\overline{\mathit{r}}\right)\right)}^2+{\left({\mathcal{U}}_y^{n+1}\left(\overline{\mathit{r}}\right)\right)}^2,$$

the nonlinear terms are linearized as follows [24]:

$${\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\equiv {\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+{\mathcal{U}}^n\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right),$$

(23)

$${\left({\mathcal{U}}_r^{n+1}\left(\overline{\mathit{r}}\right)\right)}^2\equiv 2{\mathcal{U}}_r^n\left(\overline{\mathit{r}}\right){\mathcal{U}}_r^{n+1}\left(\overline{\mathit{r}}\right)-{\left({\mathcal{U}}_r^n\left(\overline{\mathit{r}}\right)\right)}^2,$$

(24)

$${\left({\mathcal{U}}_s^{n+1}\left(\overline{\mathit{r}}\right)\right)}^2\equiv 2{\mathcal{U}}_s^n\left(\overline{\mathit{r}}\right){\mathcal{U}}_s^{n+1}\left(\overline{\mathit{r}}\right)-{\left({\mathcal{U}}_s^n\left(\overline{\mathit{r}}\right)\right)}^2.$$

(25)

The insertion of Equations (23)–(25) and simple rearrangement reduces Equation (22) to

$$\begin{array}{cc}\hfill & \left(A_{{\xi }_1}+A_{{\xi }_2}+\theta \gamma -\theta {\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\right){\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-\left(\frac{\beta }{\delta t}+\theta \delta -\theta \eta {\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\right){\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill =& \left(A_{{\xi }_1}+A_{{\xi }_2}-(1-\theta )\gamma \right){\mathcal{U}}^n\left(\overline{\mathit{r}}\right)-\left(\frac{\beta }{\delta t}-(1-\theta )\delta \right){\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+H_{\theta }^{n+1}\left(\overline{\mathit{r}}\right)-G_{\theta }^n\left(\overline{\mathit{r}}\right),\hfill \end{array}$$

(27)

where

$$\begin{array}{c}\hfill H_{\theta }^{n+1}\left(\overline{\mathit{r}}\right)=\left\{\begin{array}{c}{f}_{\theta }^{n+1}-Q_{\xi }^n{\mathcal{U}}^n\left(\overline{\mathit{r}}\right),n\ge 1,\hfill \\ f_{\theta }^1,n=1,\hfill \end{array}\right.\end{array}$$

(28)

and

$$\begin{array}{c}\hfill G_{\theta }^n\left(\overline{\mathit{r}}\right)=\eta \theta \left(2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right){\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+{\left({\mathcal{U}}_x^n\left(\overline{\mathit{r}}\right)\right)}^2+{\left({\mathcal{U}}_y^n\left(\overline{\mathit{r}}\right)\right)}^2\right)+(1-\theta )\eta \nabla \left({\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\nabla {\mathcal{U}}^n\left(\overline{\mathit{r}}\right)\right).\end{array}$$

Let us label

$$\begin{array}{cc}\hfill & p_{\theta }^n\left(\overline{\mathit{r}}\right)=A_{{\xi }_1}+A_{{\xi }_2}+\theta \gamma -\theta \eta {\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right),q_{\theta }^n=\frac{\beta }{\delta t}+\delta \theta -\theta \eta {\mathcal{U}}^n\left(\overline{\mathit{r}}\right),\hfill \\ \hfill & r_{\theta }=A_{{\xi }_1}+A_{{\xi }_2}-(1-\theta )\gamma ,s_{\theta }=\frac{\beta }{\delta t}-(1-\theta )\delta ,\hfill \end{array}$$

which, when substituted into Equation (26), give us

$$p_{\theta }^n{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)-q_{\theta }^n{\nabla }^2{\mathcal{U}}^{n+1}\left(\overline{\mathit{r}}\right)=r_{\theta }{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+s_{\theta }{\nabla }^2{\mathcal{U}}^n\left(\overline{\mathit{r}}\right)+H_{\theta }^{n+1}\left(\overline{\mathit{r}}\right)-G_{\theta }^n\left(\overline{\mathit{r}}\right).$$

(29)

We have approximated in the time direction at this point using Equation (29). Furthermore, the time direction is coupled with the splitting technique known as the Strang splitting algorithm, which is given as follows [4,57,58,59]:

$$\frac{\partial {\mathcal{U}}_1}{\partial t}={\mathcal{L}}_1{\mathcal{U}}_1\left(t\right)\mathrm{with}t\in [t^n,t^{n+1/2}]\mathrm{and}{\mathcal{U}}_1\left(t^n\right)={\mathcal{U}}_{sp}^n$$

(30)

$$\frac{\partial {\mathcal{U}}_2}{\partial t}={\mathcal{L}}_2{\mathcal{U}}_2\left(t\right)\mathrm{with}t\in [t^n,t^{n+1}]\mathrm{and}{\mathcal{U}}_2\left(t^n\right)={\mathcal{U}}_1\left(t^{n+1/2}\right)$$

(31)

$$\frac{\partial {\mathcal{U}}_3}{\partial t}={\mathcal{L}}_1{\mathcal{U}}_3\left(t\right)\mathrm{with}t\in [t^{n+1/2},t^{n+1}]\mathrm{and}{\mathcal{U}}_3\left(t^{n+1/2}\right)={\mathcal{U}}_2\left(t^{n+1}\right)$$

(32)

where ${\mathcal{L}}_1\mathcal{U}\left(t\right)=\delta {\nabla }^2\mathcal{U}\left(t\right)-\eta \left(\mathcal{U}\ast {\nabla }^2\mathcal{U}\left(t\right)\right)$, ${\mathcal{L}}_2\mathcal{U}\left(t\right)=-\eta {\left(\nabla \mathcal{U}\left(t\right)\right)}^2-\gamma \mathcal{U}\left(t\right)$, and ${\mathcal{U}}_{sp}^n={\mathcal{U}}_0$ is the initial solution defined for Equation (20). The Strang splitting method achieves second-order accuracy. Further information can be found in reference [4].

3.2. Spatial Discretization

We can expand the function $\mathcal{U}\left(\overline{\mathit{r}}\right)$, which is a two-dimensional function denoted as $\mathcal{U}(r,s)$, into the following truncated series of Lucas polynomials, expressed as a product:

$$\mathcal{U}(r,s)\simeq \sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}{L}_k\left(r\right)L_m\left(s\right)=\sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}{\mathcal{L}}_{km}(r,s),$$

(33)

where ${\mathcal{L}}_{km}(r,s)=L_k\left(r\right)L_m\left(s\right)$. Equation (33) can be expressed in compact form as

$$\mathcal{U}(r,s)\simeq {\mathbf{L}}^{\mathtt{T}}(r,s)\mathbf{\Lambda },$$

(34)

with $\mathbf{\Lambda }={[{\mathrm{\Lambda }}_{00},{\mathrm{\Lambda }}_{01},{\mathrm{\Lambda }}_{10},\cdots ,{\mathrm{\Lambda }}_{MM}]}^{\mathtt{T}}$ and ${\mathbf{L}}^{\mathtt{T}}(r,s)=[{\mathcal{L}}_{00}(r,s),{\mathcal{L}}_{01}(r,s),{\mathcal{L}}_{10}(r,s),\cdots ,{\mathcal{L}}_{MM}(r,s)]$.

One can compute and briefly represent the consolidated partial derivatives of $\mathcal{U}(r,s)$ as follows:

$$\frac{\partial \mathcal{U}(r,s)}{\partial r}\simeq \sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}\frac{\mathrm{d}{L}_k\left(r\right)}{\mathrm{d}r}{L}_m\left(s\right)=\sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}k\left(F_k\left(r\right)\right)L_m\left(s\right)={\mathbf{L}}_r^{\mathtt{T}}(r,s)\mathbf{\Lambda },$$

and

$$\frac{{\partial }^2\mathcal{U}(r,s)}{\partial r^2}\simeq \sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}\frac{{\mathrm{d}}^2{L}_k\left(r\right)}{\mathrm{d}{r}^2}{L}_m\left(s\right)=\sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}k\left(F_k\left(r\right)D\right)L_m\left(s\right)={\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)\mathbf{\Lambda },$$

where

$${\mathbf{L}}_r^{\mathtt{T}}(r,s)={\left\{\frac{\partial }{\partial r}{\mathcal{L}}_{km}(r,s)\right\}}_{k,m=0}^M={\left\{\frac{\mathrm{d}{L}_k\left(r\right)}{\mathrm{d}r}{L}_m\left(s\right)\right\}}_{k,m=0}^M={\left\{k\left(F_k\left(r\right)\right)L_m\left(s\right)\right\}}_{k,m=0}^M,$$

$${\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)={\left\{\frac{{\partial }^2}{\partial r^2}{\mathcal{L}}_{km}(r,s)\right\}}_{k,m=0}^M={\left\{\frac{{\mathrm{d}}^2{L}_k\left(r\right)}{\mathrm{d}{r}^2}{L}_m\left(s\right)\right\}}_{k,m=0}^M={\left\{k\left(F_k\left(r\right)D\right)L_m\left(s\right)\right\}}_{k,m=0}^M.$$

Similarly

$$\frac{\partial \mathcal{U}(r,s)}{\partial s}\simeq \sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}{L}_k\left(r\right)\frac{\mathrm{d}{L}_m\left(s\right)}{\mathrm{d}r}=\sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}m{L}_k\left(r\right)\left(F_m\left(s\right)\right)={\mathbf{L}}_s^{\mathtt{T}}(r,s)\mathbf{\Lambda },$$

and

$$\frac{{\partial }^2\mathcal{U}(r,s)}{\partial s^2}\simeq \sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}{L}_k\left(r\right)\frac{{\mathrm{d}}^2{L}_m\left(s\right)}{\mathrm{d}{s}^2}=\sum _{k=0}^M\sum _{m=0}^M{\mathrm{\Lambda }}_{km}m{L}_k\left(r\right)\left(F_m\left(s\right)D\right)={\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\mathbf{\Lambda },$$

with

$${\mathbf{L}}_s^{\mathtt{T}}(r,s)={\{\frac{\partial }{\partial s}{\mathcal{L}}_{km}(r,s)\}}_{k,m=0}^M={\{L_k\left(r\right)\frac{{\mathrm{d}}^2{L}_m\left(s\right)}{\mathrm{d}{s}^2}\}}_{k,m=0}^M=\{{m{L}_k\left(r\right)\left(F_m\left(s\right)\right)\}}_{k,m=0}^M.$$

$${\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)=\{{\frac{{\partial }^2}{\partial s^2}{\mathcal{L}}_{km}(r,s)\}}_{k,m=0}^M={\{L_k\left(r\right)\frac{{\mathrm{d}}^2{L}_m\left(s\right)}{\mathrm{d}{s}^2}\}}_{k,m=0}^M=\{{m{L}_k\left(r\right)\left(F_m\left(s\right)D\right)\}}_{k,m=0}^M.$$

These results can be further refined as

$$\nabla \mathcal{U}(r,s)=\frac{\partial \mathcal{U}}{\partial r}+\frac{\partial \mathcal{U}}{\partial s}\simeq \left({\mathbf{L}}_r^{\mathtt{T}}(r,s)+{\mathbf{L}}_s^{\mathtt{T}}(r,s)\right)\mathbf{\Lambda }.$$

(35)

$${\nabla }^2\mathcal{U}(r,s)=\frac{{\partial }^2\mathcal{U}}{\partial r^2}+\frac{{\partial }^2\mathcal{U}}{\partial s^2}\simeq \left({\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\right)\mathbf{\Lambda }.$$

(36)

It is consequently possible to express the time-dependent function $\mathcal{U}(r,s,t)$ and its compact form partial derivatives as follows:

$$\begin{array}{c}\hfill \mathcal{U}(r,s,t)\approx {\mathbf{L}}^{\mathtt{T}}(r,s)\mathbf{\Lambda }\left(t\right),\frac{{\partial }^2\mathcal{U}(r,s,t)}{\partial r^2}\approx {\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)\mathbf{\Lambda }\left(t\right),\frac{{\partial }^2\mathcal{U}(r,s,t)}{\partial s^2}\approx {\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\mathbf{\Lambda }\left(t\right),\end{array}$$

(37)

with

$$\begin{array}{c}\hfill \nabla \mathcal{U}(r,s,t)\approx \left({\mathbf{L}}_r^{\mathtt{T}}(r,s)+{\mathbf{L}}_s^{\mathtt{T}}(r,s)\right)\mathbf{\Lambda }\left(t\right),{\nabla }^2\mathcal{U}(r,s,t)\approx \left({\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\right)\mathbf{\Lambda }\left(t\right),\end{array}$$

(38)

where $\mathbf{\Lambda }\left(t\right)$ represents the vector of the time-varying unknown coefficients. Once the information obtained from Equations (37) and (38) is combined with Equation (29), we obtain

$$\begin{array}{cc}\hfill & p_{\theta }^n\mathbf{L}(r,s){\mathbf{\Lambda }}^{n+1}-q_{\theta }^n\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\right){\mathbf{\Lambda }}^{n+1}\hfill \\ \hfill & =r_{\theta }\mathbf{L}(r,s){\mathbf{\Lambda }}^n+r_{\theta }\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r,s)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r,s)\right){\mathbf{\Lambda }}^n+{\mathbf{H}}_{\theta }^{n+1}(r,s)-{\mathbf{G}}_{\theta }^n(r,s),\forall (r,s)\in \mathrm{\Omega }.\hfill \end{array}$$

(39)

The boundary conditions are stated as

$$\mathbb{B}\mathbf{L}(r,s){\mathbf{\Lambda }}^{n+1}\approx f_1(r,s,t^{n+1})=f_1^{n+1}(r,s),\forall (r,s)\in \partial \mathrm{\Omega },$$

(40)

where ${\mathbf{\Lambda }}^n=\mathbf{\Lambda }\left(t^n\right)$.

3.3. Full Discretization

We collocate Equations (39) and (40) at discrete mesh points $r_i$ and $s_j$ at their respective time levels to derive the discrete form of the time-fractional differential equation as expressed in Equation (20), resulting in the subsequent set of equations:

$$\begin{array}{cc}\hfill & p_{\theta }^n\mathbf{L}(r_i,s_j){\mathbf{\Lambda }}^{n+1}-q_{\theta }^n\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r_i,s_j)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r_i,s_j)\right){\mathbf{\Lambda }}^{n+1}\hfill \\ \hfill & =s_{\theta }\mathbf{L}(r_i,s_j){\mathbf{\Lambda }}^n+r_{\theta }\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r_i,s_j)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r_i,s_j)\right){\mathbf{\Lambda }}^n+{\mathbf{H}}_{\theta }^{n+1}(r_i,s_j)-{\mathbf{G}}_{\theta }^n(r_i,s_j),\forall (r_i,s_j)\in \mathrm{\Omega },\hfill \end{array}$$

and

$$\mathbb{B}\mathbf{L}(r_i,s_j){\mathbf{\Lambda }}^{n+1}\approx f_1(r_i,s_j,t^{n+1})=f_1^{n+1}(r_i,s_j),\forall (r_i,s_j)\in \partial \mathrm{\Omega },$$

where $i,j=0,1,2,\cdots ,M$. This implies the following matrix–vector form:

$$\mathcal{G}{\mathbf{\Lambda }}^{n+1}=\mathcal{H}{\mathbf{\Lambda }}^n+{\mathcal{Q}}^{n+1}.$$

(41)

In this case, the matrix entries are provided by

$$\begin{array}{c}\hfill {\left(\mathcal{G}\right)}_{ij}=\left\{\begin{array}{cc}{p}_{\theta }^n\mathbf{L}(r_i,s_j)-q_{\theta }^n\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r_i,s_j)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r_i,s_j)\right),\hfill & (r_i,s_j)\in \mathrm{\Omega },\hfill \\ \mathbb{B}\mathbf{L}(r_i,s_j),\hfill & (r_i,s_j)\in \partial \mathrm{\Omega },\hfill \end{array}\right.\end{array}$$

(42)

$$\begin{array}{c}\hfill {\left(\mathcal{H}\right)}_{ij}=\left\{\begin{array}{cc}{s}_{\theta }\mathbf{L}(r_i,s_j)+r_{\theta }\left({\mathbf{L}}_{rr}^{\mathtt{T}}(r_i,s_j)+{\mathbf{L}}_{ss}^{\mathtt{T}}(r_i,s_j)\right),\hfill & (r_i,s_j)\in \mathrm{\Omega },\hfill \\ 0,\hfill & (r_i,s_j)\in \partial \mathrm{\Omega },\hfill \end{array}\right.\end{array}$$

(43)

and

$$\begin{array}{c}\hfill {\left(\mathcal{Q}\right)}_{ij}=\left\{\begin{array}{cc}{\mathbf{H}}_{\theta }^{n+1}(r_i,s_j)-{\mathbf{G}}_{\theta }^n(r_i,s_j),\hfill & (r_i,s_j)\in \mathrm{\Omega },\hfill \\ f_1(r_i,s_j,t^{n+1}),\hfill & (r_i,s_j)\in \partial \mathrm{\Omega }.\hfill \end{array}\right.\end{array}$$

(44)

The construction of the scheme for the two-dimensional model equation is finalized. Simplifying it to a one-dimensional form involves eliminating the variable s so that $\overline{\mathit{r}}=r$. The linear system represented by Equation (41) benefits from the assurances provided by [47], ensuring the existence of a solution. Solving Equation (41) allows for the determination of the vector containing unknown coefficients at every time step. The initial condition allows for initialization in the following manner:

$${\mathbf{L}}^{\mathtt{T}}(r_i,s_j){\mathbf{\Lambda }}^0=\mathcal{U}(r_i,s_j,0)={\mathcal{U}}_0⟹\mathbf{A}{\mathbf{\Lambda }}^0={\mathcal{U}}_{\mathbf{0}}.$$

When solved, it yields ${\mathbf{\Lambda }}^0$. Utilizing Equation (41), the iteration is subsequently repeated up to the relevant time level. Upon obtaining the coefficients vector, the solution at the corresponding time level can be determined by

$$\mathcal{U}(r_i,s_j,t^n)={\mathbf{L}}^{\mathtt{T}}(r_i,s_j){\mathbf{\Lambda }}^n,\forall (r_i,s_j)\in \mathrm{\Omega },(n\ge 0).$$

A similar procedure can be used in the case of a three-dimensional problem.

4. Analyzing Errors

To conduct the error analysis, we referred to the following theorem.

Theorem 1. 

The error can be defined as follows when denoting the exact solution as $\widehat{\mathbb{U}}$ and the computed solution obtained by the suggested meshless method as $\mathcal{U}$ for the underlying problems.

The error term $\mathcal{E}$ is constrained within the following expression:

$$\left|\mathcal{E}\right|\le \frac{4exp\left(2\kappa \right){cosh}^2\left(2\mathcal{P}\right){\kappa }^{2(N+1)}}{{((N+1)!)}^2}.$$

Proof. 

Consider the absolute error between the exact solution and the approximate numerical solution using the proposed hybrid meshless method:

$$\left|\mathcal{E}\right|=|\widehat{\mathbb{U}}-\mathcal{U}|,$$

where $\widehat{\mathbb{U}}={\sum }_{k=0}^{\infty }{\sum }_{m=0}^{\infty }{\mathrm{\Lambda }}_k{\mathrm{\Lambda }}_m{L}_k\left(r\right)L_m\left(s\right)$ and $\mathcal{U}={\sum }_{k=0}^M{\sum }_{m=0}^M{\mathrm{\Lambda }}_k{\mathrm{\Lambda }}_m{L}_k\left(r\right)L_m\left(s\right)$. Next, the truncated term is as follows:

$$\left|\mathcal{E}\right|=\sum _{k=N+1}^{\infty }\sum _{m=N+1}^{\infty }{\mathrm{\Lambda }}_k{\mathrm{\Lambda }}_m{L}_k\left(r\right)L_m\left(s\right)$$

(45)

In [60], it is demonstrated that

$$L_m\left(r\right)\le 2{\vartheta }^m,$$

$$|{\mathrm{\Lambda }}_m|\le \frac{{\mathcal{P}}^{m}cosh\left(2\mathcal{P}\right)}{m!},$$

where the famous golden ratio is denoted by $\vartheta$. Consequently, Equation (45) suggests that

$$\left|\mathcal{E}\right|=4{cosh}^2\left(2\mathcal{P}\right)\sum _{k=N+1}^{\infty }\sum _{m=N+1}^{\infty }\frac{{\kappa }^{k+m}}{k!m!},$$

where $\kappa =\mathcal{P}\vartheta$. This inequality can alternatively be expressed as

$$\left|\mathcal{E}\right|=4exp\left(2\kappa \right){cosh}^2\left(2\mathcal{P}\right){\left[1-\frac{\mathrm{\Gamma }(N+1,\kappa )}{\mathrm{\Gamma }(N+1)}\right]}^2.$$

(46)

The incomplete gamma function in this case is $\mathrm{\Gamma }(N+1,\kappa )$, and the complete gamma function is $\mathrm{\Gamma }(N+1)$ [61]. Equation (46) can be expressed in integral form as follows:

$$\left|\mathcal{E}\right|\le \frac{4exp\left(2\kappa \right){cosh}^2\left(2\mathcal{P}\right)}{{(N!)}^2}{\left[{\int }_0^{\kappa }{t}^{N}exp(-t)dt\right]}^2.$$

Since $exp(-t)<1$, $\forall t>0$, we obtain:

$$\left|\mathcal{E}\right|\le \frac{4exp\left(2\kappa \right){cosh}^2\left(2\mathcal{P}\right){\kappa }^{2(N+1)}}{{((N+1)!)}^2}.$$

□

Stability Analysis

The stability of the meshless method is crucial for ensuring the accuracy and convergence of numerical solutions over time. By mitigating numerical instabilities, the meshless method maintains the solution quality and reliability across diverse computational domains. The stability analysis was carried out in this work by utilizing a matrix technique. The methodology employed was consistent with that described in previous studies [24,62].

Theorem 2. 

Consider w as the computed approximate solution to the underlying problem. In this context, the amplification matrix ϝ is defined as $\mathsf{F}=L{\mathcal{G}}^{-1}\mathcal{H}{L}^{-1}$. It is essential to confirm that the absolute maximum eigenvalue of ϝ, denoted as $\rho \left(\mathsf{F}\right)$, stays less than or equal to 1 to ensure stability.

Proof. 

We aim to show that $\mathsf{F}=L{\mathcal{G}}^{-1}\mathcal{H}{L}^{-1}$. From Equation (17), we write

$${\mathbf{\Lambda }}^n=L^{-1}{\mathcal{U}}^n.$$

(47)

The relationship can be described utilizing Equation (41):

$${\mathcal{U}}^{n+1}=L{\mathcal{G}}^{-1}\mathcal{H}{L}^{-1}{\mathcal{U}}^n+L{\mathcal{G}}^{-1}{\mathcal{Q}}^{n+1}.$$

(48)

Moreover, the error vector

$$\mathcal{E}=\widehat{\mathbb{U}}-\mathcal{U},$$

(49)

satisfies the following condition:

$${\mathcal{E}}^{n+1}=\mathsf{F}{\mathcal{E}}^n.$$

Thus, the amplification matrix is defined as $\mathsf{F}=L{\mathcal{G}}^{-1}\mathcal{H}{L}^{-1}$. Adherence to the Lax–Richtmyer stability criterion [63], which requires that

$$\left|\right|\mathsf{F}\left|\right|\le 1,$$

ensures the stability of the method. An illustration demonstrating that the spectral radius $\rho \left(\mathsf{F}\right)$ must not exceed 1 to fulfill this condition is provided in the subsequent section. □

5. Discussion and Numerical Results

The proposed hybrid meshless method (HMM) was assessed for its effectiveness and precision in solving model Equations (1) and (2). The evaluation includes various scenarios with both uniform and scattered nodes in non-rectangular and rectangular domains. The CPU time was measured in seconds for all cases. All computations were performed on an HP PC laptop equipped with an Intel(R) Core(TM) i5-7200U CPU @ 2.50 GHz 2.71 GHz and 8 GB of RAM, using MATLAB (R2012a) software.

The accuracy was measured using error norms, including the following:

$$\begin{array}{c}\hfill e_{\infty }=\max \left(\right|\widehat{\mathbb{U}}-\mathcal{U}\left|\right),e_{rms}=\sqrt{\frac{{\sum }_{i=1}^N{\left({\widehat{\mathbb{U}}}_i-{\mathcal{U}}_i\right)}^2}{N}},e_{L2}=\sqrt{\delta h\sum _{h=1}^N{\left({\widehat{\mathbb{U}}}_h-{\mathcal{U}}_h\right)}^2},\end{array}$$

(50)

where the exact solution is $\widehat{\mathbb{U}}$, while the approximate solution is $\mathcal{U}$.

Problem 1. 

The exact solution of model (1) with $\delta =1$ and $\eta =\gamma =0$ is

$$\mathcal{U}(\overline{\mathbf{r}},t)=e^{-t}sin\left(\pi r\right)sin\left(\pi s\right),\overline{\mathbf{r}}=(r,s)\in \mathrm{\Omega }.$$

(51)

The numerical results of the proposed hybrid meshless method were used to compute the numerical solution of Problem 1, as shown in

Table 1. Problem 1 numerical results using HMM and [17].

${\mathit{e}}_{\infty }$
Method	$\mathit{\delta }\mathit{t}$	$\mathit{t}=\mathbf{1}$		$\mathit{t}=\mathbf{5}$
		$\mathit{MQ}$[17]	HMM	$\mathit{MQ}$[17]	HMM
Explicit	0.1	1.8985 $\times {10}^{-2}$	8.5002 $\times {10}^{-3}$	1.5971 $\times {10}^{-3}$	8.3765 $\times {10}^{-4}$
	0.05	9.2780 $\times {10}^{-3}$	5.3215 $\times {10}^{-3}$	8.2371 $\times {10}^{-4}$	5.1759 $\times {10}^{-4}$
	0.025	4.5880 $\times {10}^{-3}$	1.0717 $\times {10}^{-3}$	4.1804 $\times {10}^{-4}$	1.0352 $\times {10}^{-4}$
	0.0125	2.2815 $\times {10}^{-3}$	9.0701 $\times {10}^{-3}$	2.1057 $\times {10}^{-4}$	9.0047 $\times {10}^{-5}$
	0.00625	1.1377 $\times {10}^{-3}$	8.0520 $\times {10}^{-4}$	1.0567 $\times {10}^{-4}$	8.6471 $\times {10}^{-5}$
Implicit	0.1	1.7500 $\times {10}^{-2}$	8.3842 $\times {10}^{-3}$	1.7904 $\times {10}^{-3}$	8.5478 $\times {10}^{-4}$
	0.05	8.9133 $\times {10}^{-3}$	5.2547 $\times {10}^{-3}$	8.7217 $\times {10}^{-4}$	5.5662 $\times {10}^{-4}$
	0.025	4.4982 $\times {10}^{-3}$	1.1406 $\times {10}^{-3}$	4.3024 $\times {10}^{-4}$	1.0639 $\times {10}^{-4}$
	0.0125	2.2594 $\times {10}^{-3}$	9.9320 $\times {10}^{-4}$	2.1364 $\times {10}^{-4}$	9.9331 $\times {10}^{-5}$
	0.00625	1.1323 $\times {10}^{-3}$	8.1481 $\times {10}^{-4}$	1.0645 $\times {10}^{-4}$	8.1101 $\times {10}^{-5}$
Crank–Nicolson	0.1	3.0659 $\times {10}^{-4}$	9.0152 $\times {10}^{-5}$	3.0957 $\times {10}^{-5}$	8.3810 $\times {10}^{-6}$
	0.05	8.4249 $\times {10}^{-5}$	1.2486 $\times {10}^{-5}$	8.4718 $\times {10}^{-6}$	1.1754 $\times {10}^{-6}$
	0.025	2.4213 $\times {10}^{-5}$	6.2692 $\times {10}^{-6}$	2.3901 $\times {10}^{-6}$	3.1455 $\times {10}^{-7}$
	0.0125	7.1414 $\times {10}^{-6}$	8.6427 $\times {10}^{-7}$	6.8938 $\times {10}^{-7}$	8.0157 $\times {10}^{-8}$
	0.00625	2.1392 $\times {10}^{-6}$	1.1503 $\times {10}^{-7}$	2.0203 $\times {10}^{-7}$	5.8023 $\times {10}^{-8}$

. The $e_{\infty }$ norm was obtained for the final times $t=1$ and $t=5$. The obtained numerical results were compared with the method reported in [17]. The results using explicit, implicit, and Crank–Nicolson schemes for the proposed method were computed for different time step sizes, and one could conclude that as the time step size $\delta t$ decreased, the accuracy of the proposed method increased. The comparison shows that the HMM was more accurate than the method given in [17]. In

Table 2. Problem 1 numerical results using HMM.

$\mathit{N}$	${\mathit{\xi }}_1={\mathit{\xi }}_2=0.1$		${\mathit{\xi }}_1={\mathit{\xi }}_2=0.5$		${\mathit{\xi }}_1={\mathit{\xi }}_2=0.9$		CPU Time
	${\mathit{e}}_{\infty }$	${\mathit{e}}_{\mathit{ems}}$	${\mathit{e}}_{\infty }$	${\mathit{e}}_{\mathit{ems}}$	${\mathit{e}}_{\infty }$	${\mathit{e}}_{\mathit{ems}}$
10	3.1204 $\times {10}^{-9}$	2.7813 $\times {10}^{-9}$	4.2101 $\times {10}^{-9}$	3.7823 $\times {10}^{-9}$	5.9105 $\times {10}^{-9}$	3.7930 $\times {10}^{-9}$	15.64
20	2.6534 $\times {10}^{-9}$	2.6591 $\times {10}^{-9}$	3.6587 $\times {10}^{-9}$	3.6901 $\times {10}^{-9}$	3.8257 $\times {10}^{-9}$	3.7133 $\times {10}^{-9}$	18.89
25	2.0673 $\times {10}^{-9}$	2.2873 $\times {10}^{-9}$	3.2510 $\times {10}^{-9}$	3.1246 $\times {10}^{-9}$	3.2780 $\times {10}^{-9}$	3.3854 $\times {10}^{-9}$	19.05

, the numerical results are presented, where they were computed by the proposed method using different fractional orders ξ, various nodal points N, $\delta t=0.001$, and $t=1$. The results were computed as $e_{\infty }$ and $e_{rms}$ norms. It can be observed that increasing the number of nodes improved the accuracy to some extent, whereas changing the fractional order did not significantly affect the accuracy. In this case, the HMM yielded good and accurate results. Additionally, the CPU time is provided in this table, demonstrating the efficiency of the suggested method.

Comparative results of the HMM and [17] are presented in Figure 1 for various nodes in the form of different error norms. It can be observed that the results produced by the HMM were superior to those of [17]. Furthermore, Figure 2 depicts the comparison of the numerically computed solution with the exact solution, along with the absolute error.

Problem 2. 

Consider the model Equation (1) with $\beta =1$ and $\gamma =\delta =0$ having the exact solution

$$\mathcal{U}(\overline{\mathbf{r}},t)=e^{r-s-t}sin\left(\pi r\right)sin\left(\pi s\right),\overline{\mathbf{r}}=(r,s)\in \mathrm{\Omega }.$$

(52)

Meshless methods offer the advantage of computing solutions without being constrained by the computational geometry. These methods naturally adjust to the geometry of the domain and offer greater flexibility in handling boundary conditions. Additionally, meshless methods often exhibit excellent scalability and parallelizability, making them suitable for large-scale simulations. To assess the true effectiveness of the proposed meshless method, various irregular geometries were examined. To evaluate the performance of the proposed HMM in non-rectangular domains, as shown in Figure 3 for Problem 2,

Table 3. Problem 2 numerical results of HMM using domains given in Figure 3.

Time	Error Norm	Domain 1	Domain 2	Domain 2	Domain 2	Domain 2
		HMM	HMM	IQ [17]	MQ [17]	IMQ [17]
t = 1	$e_{\infty }$	7.0321 $\times {10}^{-5}$	8.1482 $\times {10}^{-4}$	2.8503 $\times {10}^{-3}$	1.0762 $\times {10}^{-4}$	1.6486 $\times {10}^{-4}$
	$e_{rms}$	6.2456 $\times {10}^{-6}$	8.2453 $\times {10}^{-5}$	5.3969 $\times {10}^{-4}$	1.2498 $\times {10}^{-5}$	1.6465 $\times {10}^{-5}$
	$e_{L2}$	6.2710 $\times {10}^{-5}$	9.3608 $\times {10}^{-4}$	5.8179 $\times {10}^{-3}$	1.3473 $\times {10}^{-4}$	1.7749 $\times {10}^{-4}$
t = 2	$e_{\infty }$	3.1438 $\times {10}^{-5}$	4.1930 $\times {10}^{-5}$	6.9177 $\times {10}^{-3}$	7.4827 $\times {10}^{-5}$	1.2021 $\times {10}^{-4}$
	$e_{rms}$	4.4711 $\times {10}^{-6}$	5.591 $\times {10}^{-6}$	1.1019 $\times {10}^{-3}$	8.9038 $\times {10}^{-6}$	1.2001 $\times {10}^{-5}$
	$e_{L2}$	5.5379 $\times {10}^{-5}$	5.6283 $\times {10}^{-5}$	1.1879 $\times {10}^{-2}$	9.5983 $\times {10}^{-5}$	1.2937 $\times {10}^{-4}$
t = 3	$e_{\infty }$	9.3510 $\times {10}^{-6}$	1.8646 $\times {10}^{-5}$	1.9053 $\times {10}^{-2}$	3.9556 $\times {10}^{-5}$	6.5942 $\times {10}^{-5}$
	$e_{rms}$	1.0392 $\times {10}^{-6}$	2.8366 $\times {10}^{-6}$	2.6854 $\times {10}^{-3}$	4.8116 $\times {10}^{-6}$	6.5818 $\times {10}^{-6}$
	$e_{L2}$	2.0406 $\times {10}^{-5}$	2.9815 $\times {10}^{-5}$	2.8949 $\times {10}^{-2}$	5.1869 $\times {10}^{-5}$	7.0953 $\times {10}^{-5}$

presents the numerical results obtained using ${\xi }_1={\xi }_2=0.5$ and $\delta t=0.0001$, along with various time values, in comparison with the method outlined in [17]. The table illustrates that the HMM achieved greater accuracy compared with [17]. Additionally, in

Table 4. Problem 2 results of HMM and [17] for ${\xi }_1={\xi }_2=0.6$.

${\mathit{e}}_{\mathit{L}2}$
		t = 1	t = 1	t = 2	t = 2
Method	$\mathbf{\delta }\mathbf{t}$	HMM	[17]	HMM	[17]
Explicit	0.1	2.0148 $\times {10}^{-2}$	5.1872 $\times {10}^{-2}$	1.3867 $\times {10}^{-2}$	3.7523 $\times {10}^{-2}$
	0.05	9.5461 $\times {10}^{-3}$	2.5444 $\times {10}^{-2}$	8.4984 $\times {10}^{-3}$	1.8637 $\times {10}^{-2}$
	0.025	8.3610 $\times {10}^{-3}$	1.2612 $\times {10}^{-2}$	5.5718 $\times {10}^{-3}$	9.2941 $\times {10}^{-3}$
	0.125	3.4513 $\times {10}^{-3}$	6.2822 $\times {10}^{-3}$	1.9845 $\times {10}^{-3}$	4.6433 $\times {10}^{-3}$
	0.00625	9.3719 $\times {10}^{-4}$	3.1364 $\times {10}^{-3}$	8.3805 $\times {10}^{-4}$	2.3216 $\times {10}^{-3}$
Implicit	0.1	1.3584 $\times {10}^{-2}$	4.8894 $\times {10}^{-2}$	9.7926 $\times {10}^{-3}$	3.7120 $\times {10}^{-2}$
	0.05	9.4241 $\times {10}^{-3}$	2.4862 $\times {10}^{-2}$	7.3200 $\times {10}^{-3}$	1.8657 $\times {10}^{-2}$
	0.025	7.3721 $\times {10}^{-3}$	1.2528 $\times {10}^{-2}$	5.7313 $\times {10}^{-3}$	9.3447 $\times {10}^{-3}$
	0.0125	2.6219 $\times {10}^{-3}$	6.2845 $\times {10}^{-3}$	9.9378 $\times {10}^{-4}$	4.6732 $\times {10}^{-3}$
	0.00125	8.9866 $\times {10}^{-4}$	3.1458 $\times {10}^{-3}$	7.5566 $\times {10}^{-4}$	2.3357 $\times {10}^{-3}$
Crank–Nicolson	0.1	4.3021 $\times {10}^{-4}$	3.5095 $\times {10}^{-4}$	1.5380 $\times {10}^{-4}$	2.5348 $\times {10}^{-4}$
	0.05	7.6320 $\times {10}^{-5}$	9.2813 $\times {10}^{-5}$	7.3817 $\times {10}^{-5}$	7.1868 $\times {10}^{-5}$
	0.025	9.8431 $\times {10}^{-6}$	4.9749 $\times {10}^{-5}$	8.0116 $\times {10}^{-6}$	3.8280 $\times {10}^{-5}$
	0.0125	4.9669 $\times {10}^{-6}$	2.4637 $\times {10}^{-5}$	2.0167 $\times {10}^{-6}$	1.8688 $\times {10}^{-5}$
	0.00625	7.3499 $\times {10}^{-7}$	1.0941 $\times {10}^{-5}$	7.9173 $\times {10}^{-7}$	8.2362 $\times {10}^{-6}$

, the numerical outcomes derived from the proposed HMM were utilized to compute the $e_{L_2}$ norm for $\delta t$, which was then compared with the results obtained using the method reported in [17]. Employing explicit, implicit, and Crank–Nicolson schemes, the proposed method’s results were calculated for different time step sizes, demonstrating that decreasing $\delta t$ enhanced the accuracy of the proposed approach. This comparison highlights the superior performance of the HMM over the method described in [17].

Figure 4 presents a comparison of the HMM results with those in [17] using various error norms for different fractional orders ξ. Once again, the HMM demonstrated superior performance. Figure 5 illustrates the close agreement between the exact and numerical solutions of Problem 2 with the parameters $\xi =0.5$, $N={10}^2$, $\delta t=0.001$, and time $t=0.5$ in terms of the absolute error.

Problem 3. 

The exact solution of model (2) with $\delta =1$ and $\eta =\gamma =0$ is as follows:

$$\mathcal{U}(\overline{\mathbf{r}},t)=e^{-t}sin\left(\pi r\right)sin\left(\pi s\right)sin\left(\pi z\right),\overline{\mathbf{r}}=(r,s,z)\in \mathrm{\Omega }.$$

(53)

The results for Problem 3 were compared with the exact solution and the approach provided in [15] for different values of $\delta t$, $N=10$, and using $t=1$ and ${\xi }_1={\xi }_2=0.3$. This served to illustrate the accurateness and effectiveness of the suggested strategy. These findings are presented in

Table 5. Problem 3 numerical results of HMM and [15].

$\mathit{\delta }\mathit{t}$	${\mathit{e}}_{\mathit{\infty }}$	${\mathit{e}}_{\mathit{rms}}$	${\mathit{e}}_{\mathit{\infty }}$ [15]	${\mathit{e}}_{\mathit{rms}}$ [15]	CPU Time of HMM
0.1	7.7488 $\times {10}^{-5}$	3.7354 $\times {10}^{-5}$	2.9285 $\times {10}^{-4}$	9.2558 $\times {10}^{-5}$	1.39
0.05	2.6377 $\times {10}^{-5}$	8.5621 $\times {10}^{-6}$	7.3125 $\times {10}^{-5}$	2.3112 $\times {10}^{-5}$	1.54
0.01	6.7313 $\times {10}^{-7}$	3.5478 $\times {10}^{-7}$	2.9195 $\times {10}^{-6}$	9.2271 $\times {10}^{-7}$	2.60
0.005	2.4307 $\times {10}^{-7}$	8.5271 $\times {10}^{-8}$	7.2896 $\times {10}^{-7}$	2.3039 $\times {10}^{-7}$	3.96
0.001	8.1268 $\times {10}^{-9}$	3.2773 $\times {10}^{-9}$	2.9007 $\times {10}^{-8}$	9.1677 $\times {10}^{-9}$	17.73
0.0005	2.4501 $\times {10}^{-10}$	9.9388 $\times {10}^{-10}$	7.2236 $\times {10}^{-9}$	2.2830 $\times {10}^{-9}$	42.14

, along with the CPU time. The table shows that after a few iterations, the suggested meshless technique yielded better results. The accuracy also improved with an increase in the number of iterations, with both error norms reaching as low as ${10}^{-10}$.

The findings of Problem 3 were computed using the proposed HMM in

Table 6. Problem 3 numerical results of HMM and [15].

$\mathit{\xi }$	$\mathit{t}=1$				$\mathit{t}=2$
	${\mathit{e}}_{\infty }$	${\mathit{e}}_{\mathit{rms}}$	${\mathit{e}}_{\infty }$  [15]	${\mathit{e}}_{\mathit{rms}}$  [15]	${\mathit{e}}_{\infty }$	${\mathit{e}}_{\mathit{rms}}$	${\mathit{e}}_{\infty }$  [15]	${\mathit{e}}_{\mathit{rms}}$  [15]
0.2	8.69430 $\times {10}^{-8}$	5.3992 $\times {10}^{-8}$	7.3170 $\times {10}^{-7}$	2.3126 $\times {10}^{-7}$	7.9456 $\times {10}^{-8}$	3.9267 $\times {10}^{-8}$	5.3832 $\times {10}^{-7}$	1.7014 $\times {10}^{-7}$
0.4	9.65010 $\times {10}^{-8}$	6.5001 $\times {10}^{-8}$	7.3053 $\times {10}^{-7}$	2.3089 $\times {10}^{-7}$	8.1655 $\times {10}^{-8}$	4.9126 $\times {10}^{-8}$	5.3746 $\times {10}^{-7}$	1.6987 $\times {10}^{-7}$
0.6	9.27650 $\times {10}^{-8}$	8.4952 $\times {10}^{-8}$	7.2584 $\times {10}^{-7}$	2.2940 $\times {10}^{-7}$	7.8037 $\times {10}^{-8}$	5.7601 $\times {10}^{-8}$	5.3402 $\times {10}^{-7}$	1.6878 $\times {10}^{-7}$
0.8	8.72000 $\times {10}^{-8}$	7.3528 $\times {10}^{-8}$	7.0765 $\times {10}^{-7}$	2.2364 $\times {10}^{-7}$	7.6584 $\times {10}^{-8}$	3.9773 $\times {10}^{-8}$	5.2065 $\times {10}^{-7}$	1.6455 $\times {10}^{-7}$

, with varying values of fractional order ξ for $N=10$, $t=1$, and $t=2$. We evaluated the outcomes in terms of $e_{rms}$ and $e_{\infty }$. Accurate results were found for several time-fractional orders in this problem, and these results were also compared with the findings published in [15].

For Problem 3, we considered some non-rectangular domains shown in Figure 6 and in Figure 7 (left) to evaluate the performance of the suggested HMM. In domain 3, scattered nodes were taken into account, as shown in Figure 6 (left). The results obtained in this case are shown in Figure 8, whereas the results for domain 5 are shown in Figure 7 (right) for final times up to $t=4$. These results provide evidence of the accurate performance of the HMM in non-rectangular domains.

Problem 4. 

The exact solution for Equation (2) with $\beta =1$, $\gamma =1$, and $\delta ={\pi }^2$ is

$$\mathcal{U}(\overline{\mathbf{r}},t)=e^{t}sin\left(\pi r\right)sin\left(\pi s\right)sin\left(\pi z\right),\overline{\mathbf{r}}=(r,s,z)\in \mathrm{\Omega }.$$

(54)

For Problem 4, the behavior of the numerical and exact solutions using $N=21$ and $t=0.01$ are visualized in Figure 9 at $z=0.5$ and ${\xi }_1={\xi }_2=0.5$, along with the absolute error. One can see from these figures that the approximate solution was very compatible with the exact solution. Also, Figure 10 shows the contour graph of the problem at ${\xi }_1={\xi }_2=0.3$ and ${\xi }_1={\xi }_2=0.9$.

6. Conclusions

This paper introduces a relatively new numerical methodology for addressing two-term time-fractional PDE models in two and three dimensions. We combined the Louiville–Caputo fractional derivative scheme with the Strang splitting algorithm for the temporal component and employed a meshless technique for spatial derivatives utilizing Lucas and Fibonacci polynomials. Error norms were utilized to evaluate the accuracy of the methodology across regular and irregular domains. A comparative study was conducted between the exact solution and alternative numerical methods outlined in the contemporary literature. This inquiry illustrated that the proposed approach yielded superior performance while requiring reduced computational resources. Moreover, the study illustrated how minimizing the errors could be accomplished by refining the spatial and temporal stages. With minor modifications, the proposed approach could be applied to various intricate fractional partial differential equations that arise in fields such as fluid dynamics, quantum mechanics, and financial mathematics, where they model phenomena with significant spatial and temporal dependencies.