A hybrid collocation method for the approximation of 2D time fractional diffusion-wave equation

1.   Introduction

In recent years, spectral methods have attracted much attention from researchers. Spectral methods are very popular techniques for approximation of fractional partial differential equations (FPDEs) arising in many fields of engineering and other sciences such as bioengineering ^[1,2], anamolous diffusion ^[3,4], climate prediction models ^[5], and quantum mechanics ^[6]. Spectral methods produces very accurate solutions for the FPDEs defined on regular as well as irregular geometries ^[7,8]. These methods are highly accurate (converge exponentially) and are easily implemented for the discretization of FPDEs. A type of spectral method that is widely applicable and utilized are spectral collocation methods. In spectral collocation methods, the solution of the FPDEs is approximated in such a way that the residual vanishes at a particular set of nodes, which are known as collocation nodes. Various physical models, such as the transmission of heat, propagation of waves, diffusion. etc. consist of only one temporal fractional derivative, which is only restricted to explaining the most fundamental characteristics regarding a particular physical phenomenon. However, a single-time fractional derivative is insufficient to elaborate the complex interacting physical processes. The FPDE models comprising the multi-term time-fractional derivatives play a significant role in explaining various physical phenomena, such as the wave-diffusion in viscous and elastic solids ^[9,10], the anomalous diffusion phenomena ^[11], the oxygen delivery to tissues through a capillary ^[12], the contamination in groundwater^[13], and heat conduction ^[14], etc.

In this article, we study the Laplace transformed Chebyshev collocation method (LTCCM) for the approximation of two-dimensional MT-TFDWE in a bounded domain $\mathrm{D}$ of the form

with boundary conditions

and initial conditions

where $c_{1} > 0, c_{2}\geq 0, \mathrm{D} \subset \mathbb{R}^{2}$ presents the bounded domain with the smooth boundary $\partial \mathrm{D}$, $\mathrm{S}({\boldsymbol{\bar{\vartheta}}}, \tau)$, $\psi_{1}({\boldsymbol{\bar{\vartheta}}})$, $\psi_{2}({\boldsymbol{\bar{\vartheta}}})$, $\mathrm{b}({\boldsymbol{\bar{\vartheta}}}, \tau)$ are the given continuous functions, $\nabla^{2} = \frac{\partial^{2}}{\partial \vartheta^{2}}+\frac{\partial^{2}}{\partial \upsilon^{2}}$, $\mathfrak{B}$ is the linear boundary operator and the multi-term derivative $\mathfrak{M}_{\alpha, \alpha_{1}, \alpha_{2}, ..., \alpha_{p}}(\partial_{\tau}) \mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ is defined as

where $\alpha, \alpha_{1}, \alpha_{2} \ldots \alpha_{p} \in(1, 2),$ and $\partial_{\tau}^{\alpha} \mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ is the Caputo's time fractional derivative ^[15], defined as follows

Exact solutions are in compact form, which provides basis for comparison of the solution, and also to check out the efficiency and accuracy of the numerical scheme. Many authors obtained the exact solution of MT-TFDWEs, such as Jiang et al. ^[16] who utilized the method of separation of variables. The exact solution of MT-TFDWE by the Wavelets expansion based on LT is obtained by Liu et al. ^[17]. A semi-analytical collocation Trefftz scheme for solving MT-TFDWE is utilized in ^[18]. The authors of ^[19] obtained the analytical solution of MT-TFDWE using the finite difference method based on the L1 scheme. Daftardar et al. ^[20] solved MT-TFDWE by the Adomian decomposition method. However, the exact solution of the MT-DWEs is mostly expressed in terms of complex or multi-variate functions, i.e., Fox-H and Mittag-Leffler functions, which are extremely complicated and difficult to evaluate. Consequently, it is preferable to solve these equations numerically, especially in the case when the exact solutions are not easily available. Due to these complexities of finding the analytical solutions of FPDEs, the development of an efficient and easy-to-use numerical algorithm is the need of the day.

Numerous numerical techniques have been presented by the researchers in the last few decades, but still a lot of work needs to be done. Salehi ^[21] used a meshless point collocation method for the numerical simulation of 2D MT-TFDWE. Soltani et al. ^[22] investigated a wavelet approach for the approximation of MT-TFDWE. A new modified homotopy perturbation method (MHPM) for solving MT-TFDWE was elaborated by Jafari et al.^[23]. Wei et al. ^[24] utilized the finite element method (FEM) for the solution of 2-D MT-TFDWE. The authors of ^[25] developed a spatial sixth-order numerical scheme for solving the time-fractional diffusion equation. Jafari et al. ^[26] examined the homotopy analysis method for solving MT-TFDWE. Dehghan et al. ^[27] utilized a high order difference scheme coupled with Galerkin spectral technique for the simulation of MT-TFDWE. A finite difference scheme coupled with Legendre spectral method is utilized by Chen et al. ^[28] for numerical approximation of MT-TFDWE. Numerical methods based on the predictor–corrector method are proposed by Liu et al.^[29] for the approximation of MT-TFDWE. Li et al. ^[30] formulated a mixed finite-element method for multi-term time-fractional diffusion and diffusion–wave equations. A compact difference approach with L1 approximation for temporal and spatial discretization of MT-TFDWE was presented by Ren et al. ^[31]. Saffarian and Mohebbi ^[32] used the Galerkin spectral element method for the numerical solution of 2-D MT-TFDWE. The Galerkin finite-element method for a multi-term time-fractional diffusion equation was presented by the authors of ^[33]. The spectral method for the two-dimensional diffusion-wave equation is discussed in ^[34]. The numerical solution of the multi-term time-fractional diffusion and diffusion-wave models is utilized by Rashidinia et al. ^[35]. A spectral tau methods for distributed-order fractional diffusion equations is utilized by Zaky and Machado in ^[36].

In almost all the aforementioned numerical methods, they utilize finite element/difference schemes for the temporal discretization of FPDEs. These conventional time-stepping approaches ^[37,38] have a far higher computing cost and have critical stability constraints because the precision of these methods can only be attained for very small time steps. Furthermore, these methods are reliable and accurate if the errors remain constant or decreases as the computations move forward. The aim of this study is to implement the LTCCM for the numerical approximation of a two-dimensional multi-term time-fractional diffusion-wave equation with Dirichlet's boundary conditions. The LT is employed for temporal discretization, and the CCM is used for discretization of space variables.

In literature, it has been shown that the LT technique is an efficient alternative to the classical time-marching methods, and many authors have coupled the LT with other spatial discretization methods. The resulting hybrid techniques are computationally low, efficient, and highly accurate as compared to traditional numerical schemes. Sudicky et al. ^[39] coupled the LT with the Galerkin technique for mass transport in discretely fractured porous formations. Kamran et al. ^[40] utilized the numerical inverse LT techniques with radial basis functions for the approximation of advection diffusion models. The boundary particle method with the LT for the solution of the time fractional diffusion equations was proposed by Fu et al. ^[41]. Kuhlman ^[42] reviewed the inverse LT algorithms coupled with the boundary element method. A hybrid LT analytic method for solution of the transport models was proposed in ^[43]. Moridis et al. ^[44] coupled the LT and the finite difference method for approximation of flow models in porous medium.

On the other hand, the CCM uses Chebyshev collocation nodes, which are non-zero over the entire domain, whereas the classical time marching methods utilize the basis functions, which are non-zero over small subintervals ^[45]. Moreover, the CCM utilizes the global interpolation basis such as the Lagrange polynomial for the spatial approximation of FPDEs; thus, the CCM is a highly accurate and efficient approach. After employing the LTCCM, we obtain the solution of the FPDEs in the LT domain; thus, an NILTM is needed to obtain the solution back into the time domain, which is not easy work to tackle. Many researchers, have proposed numerical algorithms for the NILTMs, some frequently applicable and highly accurate techniques are Crump ^[46] proposed a NILTM by utilizing a Fourier series approximation. The authors of ^[47,48,49] proposed some highly accurate and efficient NILTM. Weideman and Trefethen ^[50] used parabolic and hyperbolic contours for calculation of the inverse LT. Weeks ^[51] developed NILTM by using Laguerre functions. McLean et al. ^[52] developed NILTMs for evolution equations by utilizing contour integration method. In this study, we utilize the NILTM proposed in ^[52].

2.   Existence and uniqueness results

In this section, we study the existence and uniqueness of solutions to the problem in Eqs (1)–(3). A detailed discussion on the uniqueness, existence of solution, and Ulam–Hyers stability of MT-TFDWE can be found in ^[53]. The following theorems will help in proving our main results.

Theorem 2.1. (Theorem 2.1, ^[53]) Consider a complete metric space $(Z, d) \; \mathit{\text{such that}} \; Z \neq \phi$. Consider a mapping $\mathrm{V}: Z \rightarrow Z$ such that for $z_{1}, z_{2} \in Z$, then the following inequality holds:

Then $\mathrm{V}$ has a fixed point $z^{*} \in Z$, and the fixed point is unique.

Theorem 2.2. (Theorem 2.2, ^[53]) Consider the operator $\mathcal{Y}: \mathcal{V} \rightarrow \mathcal{V}$ is continuous (completely). If

is bounded set, then, $\mathcal{Y}$ has fixed points in $\mathcal{Y}$.

Theorem 2.3. (Theorem 2.3, ^[53]) Let $C(\mathcal{V}, \mathbb{R})$ be the metric with supremum norm and $\mathcal{V}$ be a compact metric space. Then

is compact set if and only if $\mathcal{Y}$ is closed, bounded, and equi-continuous.

Lemma 2.4. ^[53] Consider the continuous functions $\mathrm{S}({\boldsymbol{\bar{\vartheta}}}, \tau): C(\mathrm{D}) \rightarrow C(\mathrm{D}), \mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau): C(\mathrm{D}, \mathbb{R}) \rightarrow C(\mathrm{D}, \mathbb{R})$, and $\mathbf{U} \in C(\mathrm{D}, \mathbb{R})$ is a solution of the integral equation

if and only if it is the solution of the problem defined in Eqs (1)–(3), where

Define the operator $\mathrm{J}: C(\mathrm{D}, \mathbb{R}) \rightarrow C(\mathrm{D}, \mathbb{R})$, that transforms, the problem to a fixed point problem, as follows:

the fixed point of $\mathrm{J}$, is the solution of the problem defined in Eqs (1)–(3). We consider the following hypotheses before proving our main results, for any $({\boldsymbol{\bar{\vartheta}}}, \tau) \in C(\mathrm{D}), \; \exists \; \; \epsilon_{j} > 0, \; j = 1, 2, 3, ..., 17$ such that

Theorem 2.5. If the assumptions $(\mathbf{H}_{\mathbf{3}}, \mathbf{H}_{\mathbf{4}})$ hold, the considered problem defined in Eqs (1)–(3) has a unique solution.

Proof. There are several steps in proving the theorem. In the first step, we will show that $\mathrm{J}$ is continuous.

Suppose a sequence $\mathbf{U}_{n} \rightarrow \mathbf{U}$, where $\mathbf{U} \in C(\mathrm{D}, \mathbb{R})$, and for $({\boldsymbol{\bar{\vartheta}}}, \tau) \in \mathrm{D}$, we can write

since $\mathbf{U}$ is continuous, therefore, we have

Hence, $\mathrm{J}$ is continuous. In the second step, we will show that $\mathrm{J}$ maps bounded sets to bounded sets, it is sufficient to prove that, for $\gamma > 0, \; \exists \; \zeta > 0$, such that

we have

for $\tau \in[0, T]$, we have

which gives

i.e.,

Which completes the proof.

In the third step, we prove that $\mathrm{J}$ is equi-continuous. For $R_{\gamma}$, a bounded set as defined in the second step, let $R_{\gamma} \subset C(\mathrm{D}, \mathbb{R})$, for $\mathbf{U} \in R_{\gamma}$, and ${\boldsymbol{\bar{\vartheta}}}^{\prime}, \; {\boldsymbol{\bar{\vartheta}}}^{\prime \prime}, \; \tau^{\prime}, \; \tau^{\prime \prime} \in \mathrm{D}$, and ${\boldsymbol{\bar{\vartheta}}}^{\prime} < {\boldsymbol{\bar{\vartheta}}}^{\prime \prime}, \; \tau^{\prime} < \tau^{\prime \prime}$, we have

finally, we obtain

The above equation is free of $\mathbf{U}$. Hence

Hence, by Theorem 2.3, $\mathrm{J}$ is completely continuous.

In the fourth step, we find a priori bound. Define $\chi = \{\mathbf{U} \in C(\mathrm{D}, \mathbb{R}): \mathbf{U} = \varepsilon \mathrm{J}\mathbf{U}, \varepsilon \in(0, 1)\}$. We prove that $\chi$ is bounded. If $\mathbf{U} \in \chi$, then $\mathbf{U} = \varepsilon \mathrm{J} \mathbf{U}$, and $0 < \varepsilon < 1$. Then for $\tau \in[0, T]$, consider

By utilizing Eq (6), we get

which further gives

Which gives that $\chi$ is bounded. Therefore, by Theorem 2.2, the operator $\mathrm{J}$ has a unique fixed point. Thus, the problem in Eqs (1)–(3) has unique a solution. □

Next, we prove that MT-TFDWE defined in Eqs (1)–(3) is Ulam–Hyres stable. The following inequality will help in developing Ulam–Hyers stability

Definition 2.6. ^[54] The solution of the problem in Eqs (1)–(3) is Ulam–Hyers stable, if for any $K \in \mathbb{R}^{+}$, such that for every approximate solution $\tilde{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \tau)$ there exists an exact solution $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ such that

Theorem 2.7. If the assumptions $(\mathbf{H}_{\mathbf{3}}, \mathbf{H}_{\mathbf{4}})$ holds, then the problem in Eqs (1)–(3) is Ulam–Hyers stable.

Proof. Let the exact solution of the problem in Eqs (1)–(3) is

and let $\tilde{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \tau)$ be the approximate solution of the problem in Eqs (1)–(3)

subtracting Eq (8) from Eq (7) and taking norm

which implies

let $\Psi = \max \bigg(|\Psi_{1}|, |\Psi_{2}|, |\Psi_{3}|, |\Psi_{4}|, |\Psi_{5}|\bigg)$, and $|h({\boldsymbol{\bar{\vartheta}}}, \tau)| < \epsilon_{14}$, where $\epsilon_{14} > 0$, then we have

by taking the infinity norm of Eq (9), we have

which further implies that

now for $t \in[0, T]$, we have

which implies

where

$\epsilon_{16} = \epsilon_{15}(\epsilon_{3}+\epsilon_{4}+1)$ and $\epsilon_{15} = \Psi \frac{T^{\alpha_{q}}}{\Gamma(1+\alpha_{q})}$, the above result can be simplified as

where $\epsilon_{17} = \frac{\epsilon_{14} \epsilon_{15}}{1-\epsilon_{16}}$, which completes the proof. □

3.   Methodology

In our numerical scheme, the MT-TFDWE model in Eqs (1)–(3), the time variable is transformed into the Laplace domain via the LT. The spatial discretization of the transformed problem is carried out by employing the CCM. To get the original solution to the problem defined in Eqs (1)–(3), we utilize the well-known numerical inversion approach, the CIM.

3.1.   Discretization of time variable via Laplace transform

For time discretization, we employ LT to Eqs. (1)–(3), such that we get a time-independent 2D MT-TFDWE. Let for $\tau > 0,$ a piecewise continuous function $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$, and let $\widehat{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$, be its LT function, is given as

The LT of Caputo's time fractional derivative ^[15] is defined as

Employing the LT on Eqs (1)–(3), we obtain

and

where

the time-independent system is obtained as

and

where $\widehat{\mathrm{T}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s}) = \mathrm{s}^{\alpha-1} \psi_{1}({\boldsymbol{\bar{\vartheta}}})+\mathrm{s}^{\alpha-2} \psi_{2}({\boldsymbol{\bar{\vartheta}}})+\mathrm{s}^{\alpha_{1}-1} \psi_{1}({\boldsymbol{\bar{\vartheta}}})+\mathrm{s}^{\alpha_{1}-2} \psi_{2}({\boldsymbol{\bar{\vartheta}}})+\ldots+\mathrm{s}^{\alpha_{p}-1} \psi_{1}({\boldsymbol{\bar{\vartheta}}})+\mathrm{s}^{\alpha_{p}-2} \psi_{2}({\boldsymbol{\bar{\vartheta}}})+\widehat{\mathrm{S}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s}), \mathrm{L} = c_{1} \nabla^{2}-c_{2} I$ and $I$ are the linear differential and identity operators, respectively. In the LTCCM, after employing the LT, we discretize $\mathrm{L}$ by utilizing CCM, and the approximate solution of the discrete system defined in Eqs (13)–(14) is obtained, and finally the solution of the original problem defined in Eqs (1)–(3) is obtained via NILTM. The CCM is discussed in next the section.

3.2.   Collocation method

Collocation methods are widely used for solving FPDEs and offer very accurate results. In collocation methods, a global polynomial interpolant is utilized on appropriate nodes, such as Chebychev collocation nodes, to extend the unknown solution of an FPDE. Numerical differentiation is a common method in scientific computing that substitutes the derivative of a function by a discrete approximation ^[55]. The related spatial derivatives are then approximated by using a discrete derivative operator, the differentiation matrices (DM) ^[56]. To provide a comprehensive overview, the major features of DM are reviewed.

In CCM, the solution is considered over [-1, 1] and interpolated $\{(\vartheta_{l}, \widehat{\mathbf{U}}(\vartheta_{l}))\}$ by the Lagrange interpolation polynomial (LP) ${Ł}_{l}(\vartheta)$ of degree less or equal to $\mathrm{n}$ ^[7,57].

where $\widehat{\mathbf{U}}_{l} = \widehat{\mathbf{U}}(\vartheta_{l})$, and ${Ł}_{l}(\vartheta)$ is the LP at collocation points $\{\vartheta_{l}\}_{l = 0}^{\mathrm{n}}$ as follows:

For the discretization of the space variables in $\mathrm{D} = [-1, 1]$, the Chebyshev points are utilized

The first derivative $\widehat{\mathbf{U}}_{\vartheta}$ is approximated by CCM as

where $\mathbf{d}_{\mathrm{n}}$ have entries of the form

and the non-diagonal entries of $\mathbf{d}_{\mathrm{n}}$ are given as

where $\varpi_{l}^{-1} = \prod_{k \neq l}^{n}(\vartheta_{k}-\vartheta_{l})$, and diagonal entries of $\mathbf{d}_{\mathrm{n}}$ are calculated by

The elements of the matrix $\mathbf{d}_{\mathrm{n}}$ of order $\beta$ are analytically calculated as

More efficient elaboration of the differentiation matrices can be found in ^[58]. Welfert ^[56] obtained an easy-to-use recursion relation for the calculation of the differentiation matrix as follows:

Now, we consider $\mathrm{D} = [-1, 1] \times[-1, 1],$ the square domain. Commonly, the Chebyshev collocation nodes are ${\boldsymbol{\bar{\vartheta}}}_{k l}$, and are presented as

The LP associated to ${\boldsymbol{\bar{\vartheta}}}_{k l}$ is written as

where $\mathrm{L}_{k l}\left({\boldsymbol{\bar{\vartheta}}}_{k l}\right) = \bigg\{\varrho_{k l}\bigg\}_{k, l = 0}^{\mathrm{n}}$ (Kronecker product). The derivatives of $2^{\text{nd}}$ order of the LP defined in Eq (17) are calculated as

where $\mathbf{d}_{\mathrm{n}}^{2}$ is DM of $2^{\text {nd }}$ order based on collocation nodes. Employing operator $\mathrm{L}$ on Eq (17) with collocation nodes ${\boldsymbol{\bar{\vartheta}}}_{r s}$, we have

finally the approximation of $\mathrm{L}$ by the CCM is given as

By using (19) in (13), we obtain

In order to incorporate the boundary conditions in Eq (14), we consider $\mathrm{L}_{\mathrm{M}}$ and all collocation nodes ${\boldsymbol{\bar{\vartheta}}}$. Moreover, the rows of $\mathrm{L}_{\mathrm{M}}$ are replaced by the corresponding nodes at the boundary with the vectors whose magnitude is one, which have unity in position corresponding to the diagonal of $\mathrm{L}_{\mathrm{M}}$. So, the boundary conditions $\mathfrak{B} \widehat{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s}) = \widehat{\mathrm{b}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$, in Eq (14) will be executed straightaway ^[7]. After re-organizing all the corresponding columns and rows of $\mathrm{L}_{\mathrm{M}}$, we obtain the matrix

where $Z$ is the square matrix of order $(\mathrm{n}-\mathrm{n}_{\mathfrak{B}}) \times(\mathrm{n}-\mathrm{n}_{\mathfrak{B}}), I$ is the square matrix of order $(\mathrm{n}_{\mathfrak{B}} \times \mathrm{n}_{\mathfrak{B}})$, and the nodes at the boundary are $\mathrm{n}_{\mathfrak{B}}$. The solution of the system defined in Eqs (13)–(14) is attained by solving the system

where, the values of interior-boundary points are accumulated via $\widehat{\mathrm{T}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$ and $\widehat{\mathrm{b}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$ correspondingly. Finally, we will utilize the NILTM to obtain the solution of the problem in Eqs (1)–(3).

3.3.   Error bound of the collocation method

As, $\mathfrak{I}_{\mathrm{n}}: C(\mathrm{D}) \rightarrow \mathrm{P}_{n}$, depend on the Chebyshev nodes (15) and Lagrange polynomial (16) as follows:

For the calculation of the error bound of CCM, we utilize the work of Borm et al. ^[59]. Suppose for all $\mathbf{U} \in C[-1, 1], \; \exists \; \mathrm{M}_{\mathrm{n}} > 0$, a constant, satisfies the inequality

Furthermore, for all $\mathbf{U} \in \mathrm{P}_{\mathrm{n}}$

For interpolation based on Chebyshev points, we have

The stability constant gets larger very sluggishly ^[59], the approximation bound is given as

Theorem 3.1. ^[59] If the polynomial interpolation error bound in (22) and the approximation bound in (25) hold for all $\mathbf{U} \in C^{\mathrm{n}+1}$, with $q \in\{0, 1, 2, \ldots, \mathrm{n}\}$, then

where $\mathrm{M}_{\mathrm{n}}^{(q)}$ is given by

The error bound is formulated by utilizing the results in Eq (25) and Eq (26): for the $1-D$ case, the linear differential operator in Eqs (1)–(3) is $\mathrm{L} = c_{1} \frac{\partial^{2}}{\partial \vartheta^{2}}-c_{2} I,$

the time derivatives are computed precisely, so the bound of error of $\|\partial_{\tau}^{\alpha}(\mathbf{U}-\mathfrak{I}_{\mathrm{n}}(\mathbf{U}))\|_{\infty}$ and $\| \partial_{\tau}^{\alpha_{q}}(\mathbf{U}-\mathfrak{I}_{\mathrm{n}}(\mathbf{U})) \|_{\infty}$, have the same orders as $\|\mathbf{U}-\mathfrak{I}_{\mathrm{n}}(\mathbf{U})\|_{\infty},$ respectively. Finally, the error bound is given by

After the calculation of coefficients of $\|\mathbf{U}^{(\mathrm{n}+1)}\|_{\infty}$, the constant $c_{1}$ is obtained. For the models in 2-D, interpolation operators based on tensor products are used to calculate the similar error bound.

3.4.   NILTM

The inversion of the LT is explained in this section. The CIM is used to numerically invert the LT. In CIM, $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ is computed as follows

here $\Upsilon$ is an appropriately selected line of integration in the complex plane linking $\sigma-i \infty$ to $\sigma+i \infty$. The integral presented in Eq (27) is called the Bromwich integral. Various numerical algorithms are used in the literature to compute the integral in Eq (27). As $\exp (\mathrm{s} \tau)$ on $\Upsilon$ is a very gradually decaying complex function, the numerical integration of the Eq (27) is very tough to execute. The contour deformation ^[49] can be utilized to handle $\exp (\mathrm{s} \tau)$, in such a way that $\Omega:(\sigma+i \infty, \sigma-i \infty)$ is deformed to Hankel's contour that begins and ends in the half plane (to the left) such that $Re(\mathrm{s}) \longrightarrow - \infty$ at both ends, so $\exp (\mathrm{s} \tau)$ depreciates very fast, therefore making the integral in Eq (27) convenient for the approximation by using the mid-point or the trapezoidal rule. This type of deformation can be established by the Cauchy's theorem in conformation to reality that $\Upsilon$ resides in the neighborhood where $\widehat{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$ is analytic. The Talbot approach may fail if the transformed function $\widehat{\mathbf{U}}({\boldsymbol{\bar{\vartheta}}}, \mathrm{s})$ has some singularities in the imaginary region of the complex plane. Eq (27) can be effectively solved using an ideal parabolic or hyperbolic contour. Several CIMs have been developed to approximate the inverse Laplace transform numerically ^[50,52]. This study employs the hyperbolic contour of the ^[52], which is expressed in parametric form as

with $\zeta > 0, \varrho \geq 0, 0 < \rho < \theta-\frac{1}{2} \pi$, and $\frac{1}{2} \pi < \theta \leq \bar{\theta} < \pi$ ^[52]. Substituting Eq (28) in Eq (27), we get

Thus Eq (29) can be approximated through the trapezoidal rule with step $k$ as

3.5.   Error analysis of LTCCM

The error analysis of the LTCCM is analyzed in this section. The Laplace transform is implemented in the first step, which is free of error. The CCM is employed in the second step for approximating the solution of the transformed problem. The following theorem establishes the error of the CCM as follows:

Theorem 3.2. (Theorem 5, pp. 48, ^[7]). Let $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ be the given function, and for a sequence $\{{\boldsymbol{\bar{\vartheta}}}_{l}\}_{\mathrm{n} = 1}^{\infty}$ based on a set of interpolation nodes, the sequence ${\boldsymbol{\bar{\vartheta}}}_{l} \rightarrow \delta$, as $n \rightarrow \infty$, where $\vartheta$ is the density function, with the respective potential $\omega$ defined by

where

For every $\mathrm{n}$, construct $P_{\mathrm{n}},$ a polynomial of the degree less than or equal to $\mathrm{n}$ that interpolates the $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ at $\{{\boldsymbol{\bar{\vartheta}}}_{i}\}_{\mathrm{n} = 1}^{\infty}$. If $\; \exists \; \omega_{\mathbf{U}} > \omega_{[-1, 1]}$, where $\omega_{\mathbf{U}}$ is a constant, such that $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ is differentiable at each point within the region defined by

suppose for $\mathrm{M} > 0$ a constant, such that $\forall \; \; {\boldsymbol{\bar{\vartheta}}}$ and $\forall \; \; \mathrm{n}$,

The above estimate is valid for any order derivatives $(\mathbf{U}^{\alpha}-P_{\mathrm{n}}^{\alpha},$ where $\alpha \geq 1)$ with a new constant, which will still not be dependent on $\mathrm{n}$ and $x$.

Finally, we employ the CIM in order to approximate the Eq (27). While approximating Eq (27) the convergence of the proposed scheme depends on $\Upsilon$ and the set of optimal parameters. The error analysis of the CIM is discussed in the following theorem.

Theorem 3.3. (^[52], Theorem 2.2). Let $\mathbf{U}({\boldsymbol{\bar{\vartheta}}}, \tau)$ be the solution of Eq (1), and let $\widehat{\mathrm{S}}$ be differentiable in the set $\sum_{\theta}^{\varrho}$. Set $\gamma = 2 \vartheta, \quad \vartheta \in(0, 1]$, and let $\Upsilon \subset S_{r} \subset \sum_{\theta}^{\varrho}$ be defined as $\zeta = \frac{\gamma}{\kappa T}$, with $\kappa = 1-\sin (\rho-r)$. Suppose $\mathbf{U}_{k}({\boldsymbol{\bar{\vartheta}}}, \tau)$ be the approximate solution to Eqs (1)–(3) with $k = \sqrt{\frac{\bar{r}}{\gamma \mathrm{m}}} \leq \frac{\bar{r}}{\log 2}$. Then, if $\vartheta_{0}, \nu \geq 0$ and $\vartheta_{0}+\frac{1}{2} \nu \geq \vartheta$, we have, with $C = C_{\rho, r, \theta, \vartheta, \vartheta_{0}}$ for $0 \leq \tau \leq T$,

In the current work, we use the optimal contour of integration defined in Eq (28) with the best parameters suggested in ^[52] as follows:

therefore, the error estimate is given as

4.   Application

We consider four different test problems to validate and check the efficiency of the LTCCM. The maximum absolute error $(\varepsilon r^{\infty})$ norm is computed among the numerical solution and the analytical solution to measure the accuracy of the suggested LTCCM technique. The error norm is defined as:

where $\mathbf{U}$ is the analytical solution and $\mathbf{U}_{k}$ is the approximate solution.

4.1.   Test Example 1

Consider the 2D two-term time fractional diffusion-wave equation ^[35] with the exact solution $\tau^{2+\alpha+\alpha_{1}} e^{\vartheta+\upsilon},$ and $c_{1} = 1, c_{2} = 0$ of the form

with initial and Dirichlet's boundary conditions

where $(\vartheta, \upsilon, \tau) \in[-1, 1] \times[-1, 1] \times[0, 1]$. The results presented in Table and figures are computed at $\tau = 1$. The computed outcomes of the $er^{\infty}$ error for the varying values of $\mathrm{n}$ and $\mathrm{m}$ with $\alpha = 1.90$ and $\alpha = 1.50$, and $\alpha_1 = 1.20$ and $\alpha_{1} = 1.70$ are displayed in Table 1. The results presented in Table 1 are compared with the Legendre collocation method ^[35], which clearly demonstrates that the proposed method is highly accurate. Moreover, for a very small number of points $\mathrm{n}$ a very reasonable accuracy is obtained. The numerical solution of test example 1 determined by the proposed method with $\mathrm{n} = 25, \mathrm{\; m} = 90, \alpha = 1.45$, and $\alpha_{1} = 1.85$ is presented in Figure 1a. In Figure 1b, the graphs of the $er^{\infty}$ for various values of $\alpha$, and $\alpha_{1}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\} \; \text{and} \; \mathrm{m} = 90$ are displayed. In Figure 1c, the plots of the $er^{\infty}$ for various values of $\mathrm{m}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\}, \alpha = 1.45$, and $\alpha_{1} = 1.85$ are displayed. The graphs of the $er^{\infty}$ for various values of $\mathrm{n}$ with $\mathrm{m} = \{50, 55, 60, \ldots, 100\}, \alpha = 1.45$, and $\alpha_{1} = 1.85$ are presented in Figure 1d.

4.2.   Test Example 2

Consider the 2D two-term time fractional diffusion–wave equation with the exact solution $\tau^{2} \cos \bigg(\frac{\pi}{2} \vartheta\bigg) \cos \bigg(\frac{\pi}{2} \upsilon \bigg),$ and $c_{1} = 1, c_{2} = 0$ of the form

with initial and Dirichlet's boundary conditions

where $(\vartheta, \upsilon, \tau) \in[-1, 1] \times[-1, 1] \times[0, 1]$. The results presented in Table and figures are computed at $\tau = 1$. The computed outcomes of the $er^{\infty}$ error for the varying values of $\mathrm{n}$ and $\mathrm{m}$ with $\alpha = 1.50, \alpha_{1} = 1.70$ are displayed in Table 2. The results presented in Table 2 clearly demonstrates that the proposed method is highly accurate. Moreover, for very low values of $\mathrm{n},$ a reasonable accuracy is achieved. The numerical solution of test example 2 determined by the proposed method with $\mathrm{n} = 25, \mathrm{\; m} = 90, \alpha = 1.45$, and $\alpha_{1} = 1.85$ is presented in Figure 2a. In Figure 2b, the graphs of the $er^{\infty}$ for various values of $\alpha$, and $\alpha_{1}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\}\; \text{and} \; \mathrm{m} = 90$ are displayed. In Figure 2c, the plots of the $er^{\infty}$ for various values of $\mathrm{m}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\}, \alpha = 1.45$, and $\alpha_{1} = 1.85$ are displayed. The graphs of the $er^{\infty}$ for various values of $\mathrm{n}$ with $\mathrm{m} = \{50, 55, 60, \ldots, 100\}, \; \alpha = 1.45$, and $\alpha_{1} = 1.85$ are presented in Figure 2d.

4.3.   Test Example 3

Consider the 2D three-term time fractional diffusion–wave equation ^[32] with the exact solution $\tau^{3} \sin (\pi \vartheta) \sin (\pi \upsilon),$ and $c_{1} = c_{2} = 1$ of the form

with initial and Dirichlet's boundary conditions

where $(\vartheta, \upsilon, \tau) \in[-1, 1] \times[-1, 1] \times[0, 1]$. The results presented in Table and figures are computed at $\tau = 1$. The computed outcomes of the $er^{\infty}$ error for the varying values of $\mathrm{n}$ and $\mathrm{m}$ with $\alpha = 1.50$ and $\alpha = 1.30$, $\alpha_{1} = 1.30$ and $\alpha_{1} = 1.50$, and $\alpha_{2} = 1.15$ and $\alpha_{2} = 1.95$ are displayed in Table 3. The results presented in Table 3 are compared with the Galerkin spectral element method ^[32], which clearly demonstrates that the proposed method is highly accurate. Moreover, for very small number of nodes $\mathrm{n}$ a reasonable accuracy is obtained. The numerical solution of test example 3 determined by the proposed method with $\mathrm{n} = 25, \mathrm{\; m} = 90, \alpha = 1.45, \alpha_{1} = 1.65$, and $\alpha_{2} = 1.95$ is presented in Figure 3a. In Figure 3b, the graphs of the $er^{\infty}$ for various values of $\alpha, \alpha_{1}$, and $\alpha_{2}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\} \; \text{and} \; \mathrm{m} = 90$ are displayed. InFigure 3c, the plots of the $er^{\infty}$ for various values of $\mathrm{m}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\}, \alpha = 1.45, \alpha_{1} = 1.65$, and $\alpha_{2} = 1.95$ are displayed. The graphs of $er^{\infty}$ for various values of $\mathrm{n}$ with $\mathrm{m} = \{50, 55, 60, \ldots, 100\}, \alpha = 1.45, \alpha_{1} = 1.65$, and $\alpha_{2} = 1.95$ are presented in Figure 3d.

4.4.   Test Example 4

Consider the 2D four-term time fractional diffusion–wave equation with the exact solution $(1+\tau+\tau^{2})\bigg(\frac{\vartheta^{3}}{3}-\vartheta^{2}+\frac{\vartheta}{3}+\frac{1}{3}\bigg)+\bigg(\frac{\upsilon^{3}}{3}-\upsilon^{2}+\frac{\upsilon}{3}+\frac{1}{3}\bigg),$ and $c_{1} = 1, c_{2} = 1$ of the form

with initial and Dirichlet's boundary conditions

where $(\vartheta, \upsilon, \tau) \in[-1, 1] \times[-1, 1] \times[0, 1]$. The results presented in Table and figures are computed at $\tau = 1$. The computed outcomes of the $er^{\infty}$ error for the varying values of $\mathrm{n}$ and $\mathrm{m}$ with $\alpha = 1.30, \alpha_{1} =$ 1.50, $\alpha_{2} = 1.75$, and $\alpha_{3} = 1.95$ are displayed in Table 4. The results presented in Table 4 clearly demonstrate that the proposed method is highly accurate. Further, for very small number of points $\mathrm{n}$ accurate results are obtained. The numerical solution of test example 4 determined by the proposed method with $\mathrm{n} = 25, \mathrm{\; m} =$ 120, $\alpha = 1.25, \alpha_{1} = 1.45, \alpha_{2} = 1.75$, and $\alpha_{3} = 1.95$ is shown in Figure 4a. In Figure 4b, the plots of the $er^{\infty}$ for various values of $\alpha, \alpha_{1}, \alpha_{2}$, and $\alpha_{3}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\} \; \text{and} \; \mathrm{m} = 120$ are displayed. In Figure 4c, the plots of the $er^{\infty}$ for various values of $\mathrm{m}$ with $\mathrm{n} \in\{6, 8, 10, \ldots, 16\}, \alpha = 1.55, \alpha_{1} = 1.75, \alpha_{2} = 1.95$, and $\alpha_{3} = 1.85$ are displayed. The graphs of the $er^{\infty}$ for various values of $\mathrm{n}$ with $\mathrm{m} = \{50, 55, 60, \ldots, 100\}, \alpha = 1.25, \alpha_{1} =$ $1.45, \alpha_{2} = 1.75$, and $\alpha_{3} = 1.95$ are presented in Figure 4d.

5.   Conclusions

In this article, we have successfully implemented an efficient numerical scheme for the numerical modeling of the time-fractional diffusion-wave equation. A robust framework was developed by combining the LT with the CCM and employing the contour integration method for inverting the LT numerically. The suggested, technique proved to be capable of handling complex fractional dynamics. We tested the proposed approach on single-term, two-term, and three-term fractional problems. The outcomes of these tests suggest that the proposed method produces accurate and stable solutions. The solution's spatial and temporal features are accurately captured by combining the efficient handling of spatial derivatives using the CCM and the time fractional derivative using the LT. The contour integration method further increased the adaptability of our numerical technique by providing precise inversion of the Laplace transform, which is vital for handling the time-fractional derivatives.

Author contributions

Conceptualization, F.A.S and K; Methodology, F.A.S and K; software, F.A.S; validation, K, Z.A.K and F.A; formal analysis, F.A and N.M; investigation, K., F.A and N.M; writing{original draft preparation, F.A.S, K and Z.A.K; data curation, F.A.S, Z.A.K and F.A; visualization, K, Z.A. K, F.A and N. M; funding acquisition, F.A and N.M. All authors reviewed the manuscript.

Acknowledgments

The author Z.A. Khan expresses her gratitude to Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2024R8). Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

The author F. Azmi and N. Mlaiki would like to thank Prince Sultan University for paying the APC and for the support through the TAS research laboratory.

Conflict of interest

The authors declare no conflicts of interest.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.