A Novel and Accurate Algorithm for Solving Fractional Diffusion-Wave Equations

Abstract

The main objective of this work is to apply a novel and accurate algorithm for solving the second-order and fourth-order fractional diffusion-wave equations (FDWEs). First, the desired equation is reduced to the corresponding Volterra integral equation (VIE). Then, the collocation method is applied, for which the Chebyshev cardinal functions (CCFs) have been considered as the bases. In this paper, the CCFs based on a Lobatto grid are introduced and used for the first time to solve these kinds of equations. To this end, the derivative and fractional integral operators are represented in CCFs. The main features of the method are simplicity, compliance with boundary conditions, and good accuracy. An exact analysis to show the convergence of the scheme is presented, and illustrative examples confirm our investigation.

1. Introduction

Recently, the modeling of various physical phenomena greatly relies on the significance of fractional differential equations (FDEs), making this subject a focal point for many researchers [1,2,3,4]. Numerous mathematical schemes have been introduced and applied to study and solve FDEs, including the adaptive grid technique [5], the B-spline collocation method [6], the wavelet method [7,8], multi-step methods [9], fast second-order accurate difference schemes [10], the Petrov–Galerkin finite element-meshfree formulation [11], the finite element method [12], an implicit integration factor method [13], Adomian decomposition [14], the Kuratowski measure of the noncompactness technique [15], etc.

The FDWE serves as a valuable mathematical model for representing various significant physical phenomena. The foundation of this model lies in a linear integro-partial fractional differential equation derived from the classical diffusion-wave equation. Notably, in its formulation, the conventional second-order time derivative is substituted with a fractional derivative of order $\kappa$ ($\kappa \in (1,2]$). The FWDE accurately describes numerous universal mechanical, electromagnetic, and acoustic responses; see, e.g., [16]. The fourth-order spatial derivative is accrued in the context of wave propagation in beams and the modeling of grooves on a flat surface. As a result, substantial attention has been focused on the fourth-order FDWE and its practical applications [17]. The second-order and fourth-order FDWEs have been considered in this paper as follows.

Our objective in this study was to solve the second-order FDWE,

(1)${}_0{}^C\mathcal{D}_t^{\kappa }\left(y\right)(x,t)={\mathcal{D}}_x^2\left(y\right)(x,t)+g(x,t),t\in [0,T],x\in [0,b],$

with

$\left\{\begin{array}{ccc}y(x,0)={\varphi }_0\left(x\right),\hfill & {\mathcal{D}}_t\left(y\right)(x,0)={\varphi }_1\left(x\right),\hfill & x\in [0,b],\left(\mathrm{initial}\mathrm{conditions}\right),\hfill \\ y(0,t)={\psi }_0\left(t\right),\hfill & y(1,t)={\psi }_1\left(t\right),\hfill & t\in [0,T],\left(\mathrm{boundary}\mathrm{conditions}\right),\hfill \end{array}\right.$

and the fourth-order FDWE,

(2)${}_0{}^C\mathcal{D}_t^{\kappa }\left(y\right)(x,t)+{\mathcal{D}}_x^4\left(y\right)(x,t)=g(x,t),t\in [0,T],x\in [0,b],$

with

$\left\{\begin{array}{ccc}y(x,0)={\varphi }_0\left(x\right),\hfill & {\mathcal{D}}_t\left(y\right)(x,0)={\varphi }_1\left(x\right),\hfill & x\in [0,b],\left(\mathrm{initial}\mathrm{conditions}\right),\hfill \\ y(0,t)={\psi }_0\left(t\right),\hfill & y(1,t)={\psi }_1\left(t\right),\hfill & \hfill \\ {\mathcal{D}}_x^2\left(y\right)(0,t)={\psi }_2\left(t\right),\hfill & {\mathcal{D}}_x^2\left(y\right)(1,t)={\psi }_3\left(t\right),\hfill & \left(\mathrm{boundary}\mathrm{conditions}\right),\hfill \end{array}\right.$

where $\kappa \in (1,2]$, and $\mathcal{D}$ indicates the derivative operator.

Several schemes have been introduced and implemented to solve these kinds of equations, including a compact difference scheme [18,19], the finite difference scheme [20], the Sumudu transform method [21], high-order compact finite difference methods [22], the implicit difference scheme [23], and the spectral tau method [24].

We strive to present a novel and accurate algorithm for solving Equations (1) and (2). After converting the desired equation into a Volterra integral equation, the next step is to solve it using the collocation method that relies on the CCFs. The CCFs are the powerful bases for solving a variety of equations [6,25,26]. In this work, the CCFs obtained from a Lobatto grid are used for the first time. Notably, the Lobatto grid is popular for solving boundary value problems because the boundary conditions instantly furnish two grid point values for the unknown. On the other hand, when the differential equation is singular at the endpoints, then the Lobatto grid is preferable.

The upcoming portions of the document are arranged as follows: Section 2 contains a review and introduction of CCFs and their properties. The collocation method is applied to solve second-order and fourth-order FDWEs using CCFs in Section 3. This section includes an investigation of the convergence analysis as well. Section 4 is dedicated to illustrating the accuracy and applicability of the method. To sum up our work, we included a conclusion in Section 5.

2. Chebyshev Cardinal Polynomials Based on Lobatto Grid

Let us denote the Chebyshev polynomials of order N by $T_N$. As we know, $T_N$ has N real, simple, and different roots, sometimes called Chebyshev nodes. The Chebyshev nodes can be set as ${\left\{x_n\right\}}_{n=1}^N$. Defining the Chebyshev polynomials onto the interval $[a,b]$, called shifted Chebyshev polynomials and characterized by $T_N^{*}$, is easy, and can be carried out through an affine transformation $x\mapsto {\displaystyle \frac{2(x-a)}{b-a}}-1$ that maps $[-1,1]$ into $[a,b]$, viz.,

(3)$T_N^{*}\left(x\right):=T_N\left({\displaystyle \frac{2(x-a)}{b-a}}-1\right).$

Consequently, the roots of $T_N^{*}$ are obtained through

(4)$x_n^{*}={\displaystyle \frac{1}{2}}(x_n+1)(b-a)+a,\forall 1\le n\le N.$

The CCFs defined rely on two types of grids of the Chebyshev nodes. The first definition is associated with the Gauss–Chebyshev nodes of $T_{N+1}^{*}$, and the second involves the Lobatto grid, which includes the extrema of Chebyshev polynomial $T_N^{*}$ along with the endpoints. Considering ${\mathcal{D}}_x$ as the derivative operator with respect to x, the CCFs corresponding to the Gauss–Chebyshev grid [25,26,27] are stated by the following:

(5)$C_n\left(x\right)={\displaystyle \frac{T_{N+1}^{*}\left(x\right)}{{\mathcal{D}}_x\left(T_{N+1}^{*}\right)\left(x_n^{*}\right)(x-x_n^{*})}},x\in [a,b],1\le n\le N+1.$

As mentioned above, there is another definition of CCFs with a substitute of nodes, namely the Lobatto grid, which are determined as follows:

(6)$C_n\left(x\right)={\displaystyle \frac{X\left(x\right){\mathcal{D}}_x\left(T_N^{*}\right)\left(x\right)}{{\mathcal{D}}_x\left(X{T}_N^{*}\right)\left(x_n^{*}\right)(x-x_n^{*})}},x\in [a,b],1\le n\le N+1,$

in which $X\left(x\right)=1-{({\displaystyle \frac{2(x-a)}{b-a}}-1)}^2$.

For $\nu \in \mathbb{N}$, the Sobolev space ${\mathcal{H}}^{\nu }\left([a,b]\right)$ is denoted by

${\mathcal{H}}^{\nu }\left([a,b]\right)=\left\{y\in C^{\nu }\left([a,b]\right):{\mathcal{D}}_x^{\sigma }y\in L^2\left([a,b]\right),\sigma \le \nu ,\sigma \in {\mathbb{N}}_0\right\},$

endowed with the inner product

(7)${(y_1,y_2)}_{{\mathcal{H}}^{\nu }\left([a,b]\right)}=\sum _{\sigma =0}^{\nu }{({\mathcal{D}}_x^{\sigma }\left(y_1\right),{\mathcal{D}}_x^{\sigma }\left(y_2\right))}_{L^2\left([a,b]\right)},$

which induces the norm

(8)${\parallel y\parallel }_{{\mathcal{H}}^{\nu }\left([a,b]\right)}^2=\sum _{\sigma =0}^{\nu }{\parallel {\mathcal{D}}_x^{\sigma }\left(y\right)\left(x\right)\parallel }_2^2,$

and the semi-norm

(9)${\left|y\right|}_{{\mathcal{H}}^{\nu ,N}\left([a,b]\right)}^2=\sum _{\sigma =min\{\nu ,N\}}^N{\parallel {\mathcal{D}}_x^{\sigma }\left(y\right)\left(x\right)\parallel }_2^2.$

The remarkable characteristic of the CCFs is their cardinality, i.e.,

(10)$C_n\left(x_{n^{\prime }}\right)={\delta }_{n{n}^{\prime }},1\le n,n^{\prime }\le N+1,$

where ${\delta }_{n{n}^{\prime }}$ denotes the Kronecker delta. The significance of this property lies in its ability to help determine coefficients without the need for calculating integrals when approximating any function $y\in {\mathcal{H}}^{\nu }\left([a,b]\right)$ with them. Introducing a projection operator ${\mathcal{Q}}_N$ is essential to map $y\in C[a,b]$ into $P_{N+1}$ (the space of all polynomials of degree less than $N+1$), i.e.,

(11)$y\left(x\right)\approx {\mathcal{Q}}_N\left(y\right)\left(x\right)=\sum _{n=1}^{N+1}y\left(x_n^{*}\right)C_n\left(x\right)\in P_{N+1}.$

Lemma 1 

(cf [28]). Given $N\ge 0$, if $y\in {\mathcal{H}}^{\nu }\left([a,b]\right)$, we have

(12)$\parallel y-{\mathcal{Q}}_N{\left(y\right)\parallel }_2\le C{(b-a)}^{\nu }{N}^{-\nu }{\left|y\right|}_{{\mathcal{H}}^{\nu ,N}\left([a,b]\right)},$

where C is a constant independent of N.

2.1. Matrix Representation of the Derivative Operator in CCFs Obtained by Lobatto Grids

Consider $\mathbf{C}\left(x\right)$ as a vector function whose components are ${\left\{C_n\right\}}_{1\le n\le N+1}$. The focus of this part is to specify a square matrix D, denoted as

(13)${\mathcal{D}}_x\left(\mathbf{C}\right)\left(x\right)=D\mathbf{C}\left(x\right).$

The process below outlines how to specify the entries for D. Thanks to Equation (11), it is easy to confirm that

(14)${\left[D\right]}_{n,n^{\prime }}={\mathcal{D}}_x\left(C_n\right)\left(x_{n^{\prime }}^{*}\right).$

According to the definition of the CCFs using the Lobbato grid, an alternate formulation can be presented (see, e.g., [29]):

(15)${\displaystyle C_n\left(x\right)={\rho }_n\prod _{\begin{array}{c}{n}^{\prime }=1\\ n^{\prime }\ne n\end{array}}^{N+1}(x-x_{n^{\prime }}^{*}),}$

where

${\rho }_n={\displaystyle \frac{-4\Gamma (2N-2)}{{(b-a)}^{N+1}{\Gamma }^2(N-1)\left({\mathcal{D}}_x\left(X{T}_N^{*}\right)\left(x_n^{*}\right)\right)}}.$

Motivated by (6), taking the derivative from both sides of (15) leads to the following:

(16)${\displaystyle {\mathcal{D}}_x\left(C_n\right)\left(x\right)={\rho }_n\prod _{\begin{array}{c}{n}^{\prime }=1\\ n^{\prime }\ne n\end{array}}^{N+1}{\mathcal{D}}_x(x-x_{n^{\prime }}^{*}).}$

• If $n=n^{\prime }$, then we obtain

  (17)$\begin{array}{cc}\hfill {\mathcal{D}}_x\left(C_n\right)\left(x\right)& ={\rho }_n\sum _{\begin{array}{c}{n}^{\prime \prime }=1\\ n^{\prime \prime }\ne n\end{array}}^{N+1}\prod _{\begin{array}{c}{n}^{\prime }=1\\ n^{\prime }\ne n,n^{\prime \prime }\end{array}}^{N+1}(x-x_{n^{\prime }}^{*})\hfill \\ \hfill & =\sum _{\begin{array}{c}{n}^{\prime \prime }=1\\ n^{\prime \prime }\ne n\end{array}}^{N+1}{\displaystyle \frac{C_n\left(x\right)}{(x-x_{n^{\prime \prime }}^{*})}}.\hfill \end{array}$

  From Equation (10), we conclude that

  (18)${\mathcal{D}}_x\left(C_n\right)\left(x_{n^{\prime }}^{*}\right)=\sum _{\begin{array}{c}{n}^{\prime \prime }=1\\ n^{\prime \prime }\ne n\end{array}}^{N+1}{\displaystyle \frac{1}{(x_{n^{\prime }}-x_{n^{\prime \prime }})}}.$

• If $n\ne n^{\prime }$, then one can obtain

  (19)${\displaystyle {\mathcal{D}}_x\left(C_n\right)\left(x_{n^{\prime }}^{*}\right)={\rho }_n\prod _{\begin{array}{c}{n}^{\prime \prime }=1\\ n^{\prime \prime }\ne n,n^{\prime }\end{array}}^{N+1}(x_{n^{\prime }}^{*}-x_{n^{\prime \prime }}^{*}).}$

2.2. Matrix Representation of the Fractional Integral Operator in CCFs Obtained by Lobatto Grids

In recalling the fractional integral operator (FIO), see, e.g., [30,31],

(20)${\mathcal{I}}_0^{\kappa }\left(y\right)\left(x\right)={\displaystyle \frac{1}{\Gamma \left(\kappa \right)}}{\int }_0^x{(x-s)}^{\kappa -1}y\left(s\right)ds,x\in [0,b],\kappa \in {\mathbb{R}}^{+},$

the square matrix $I_{\kappa }$ fulfills

(21)${\mathcal{I}}_0^{\kappa }\left(\mathbf{C}\left(x\right)\right)\approx I_{\kappa }\mathbf{C}\left(x\right),x\in (0,b).$

Thanks to (11), one can calculate the elements of $I_{\kappa }$ as

(22)${\left[I_{\kappa }\right]}_{n,n^{\prime }}={\mathcal{I}}_0^{\kappa }\left(C_n\left(x_{n^{\prime }}^{*}\right)\right).$

For simplicity, the alternate formula of (15) can be considered, introduced in [25], as follows:

(23)${\displaystyle C_n\left(x\right)={\rho }_n\prod _{\begin{array}{c}{n}^{\prime }=1\\ n^{\prime }\ne n\end{array}}^{N+1}(x-x_{n^{\prime }}^{*})={\rho }_n\sum _{n^{\prime }=0}^N{\varrho }_{n,n^{\prime }}{x}^{N-n^{\prime }}}$

where

${\varrho }_{n,0}=1,{\varrho }_{n,n^{\prime }}={\displaystyle \frac{1}{n^{\prime }}}\sum _{n^{\prime \prime }=0}^{n^{\prime }}{\varsigma }_{n,n^{\prime \prime }}{\varrho }_{n,n^{\prime }-n^{\prime \prime }},n^{\prime }=1,\dots ,N,n=1,\dots ,N+1,$

with

${\varsigma }_{n,n^{\prime \prime }}=\sum _{\begin{array}{c}i=1\\ i\ne n\end{array}}^{N+1}{\left(x_i^{*}\right)}^{n^{\prime \prime }},n^{\prime \prime }=1,\dots ,N,n=1,\dots ,N+1.$

Taking into account the Property 2.1 [30] and using (23), one can write

(24)$\begin{array}{cc}\hfill {\mathcal{I}}_0^{\kappa }\left(C_n\left(x\right)\right)& ={\rho }_n{\mathcal{I}}_0^{\mu }(\sum _{n^{\prime }=0}^N{\varrho }_{n,n^{\prime }}{x}^{N-n^{\prime }})\hfill \\ \hfill & ={\rho }_n\sum _{n^{\prime }=0}^N{\varrho }_{n,n^{\prime }}{\mathcal{I}}_0^{\kappa }\left(x^{N-n^{\prime }}\right)\hfill \\ \hfill & ={\rho }_n\sum _{n^{\prime }=0}^N{\varrho }_{n,n^{\prime }}{\displaystyle \frac{(N-n^{\prime })!}{\Gamma (N-n^{\prime }+\kappa +1)}}{x}^{N-n^{\prime }+\kappa }.\hfill \end{array}$

Consequently, (22) and (24) lead to determining

(25)${\left[I_{\kappa }\right]}_{n,n^{\prime }}={\rho }_n\sum _{n^{\prime }=0}^N{\varrho }_{n,n^{\prime }}{\displaystyle \frac{(N-n^{\prime })!}{\Gamma (N-n^{\prime }+\kappa +1)}}{\left(x_{n^{\prime }}^{*}\right)}^{N-n^{\prime }+\kappa }.$

3. Proposed Algorithm

Let us combine Equations (1) and (2) and refer to them as

(26)${}_0{}^C\mathcal{D}_t^{\kappa }\left(y\right)(x,t)+{(-1)}^{m/2}{\mathcal{D}}_x^m\left(y\right)(x,t)=g(x,t),\mathrm{with}m=2,4.$

It is clear that Equation (26) corresponds to (1) when $m=2$, and it corresponds to (2) when $m=4$.

Taking into account the Lemma 2.22 [30] and taking the fractional integral from both sides of (26), one can obtain the corresponding integral equation, viz.,

(27)$y(x,t)-v(x,t)+{\mathcal{I}}_0^{\kappa }\left({(-1)}^{m/2}{\mathcal{D}}_x^m\left(y\right)-g\right)(x,t)=0,t\in [0,T],x\in [0,b],$

where $1<\kappa \le 2$, and $v(x,t)=t{\mathcal{D}}_t\left(y\right)(x,0)+y(x,0)$. To give rise to the proposed algorithm, the solution y must be approximated by CCFs, i.e.,

(28)$y(x,t)\approx {\mathcal{Q}}_{N,N^{\prime }}\left(y\right)(x,t)=\sum _{n=1}^{N+1}\sum _{n^{\prime }=1}^{N^{\prime }+1}{y}_{n,n^{\prime }}{C}_n\left(x\right)C_{n^{\prime }}\left(t\right):=\tilde{y}(x,t),$

where N and $N^{\prime }$ can be different. Strictly speaking, they depend on our selection of the number of bases in each axis. Substituting $\tilde{y}$ in integral Equation (27) leads to

(29)$\tilde{y}(x,t)-\tilde{v}(x,t)+{\mathcal{I}}_0^{\kappa }\left({(-1)}^{m/2}{\mathcal{D}}_x^m\left(\tilde{y}\right)-\tilde{g}\right)(x,t)=0,$

where $\tilde{v}(x,t)={\mathcal{Q}}_{N,N^{\prime }}\left(v\right)(x,t)$, and $\tilde{g}(x,t)={\mathcal{Q}}_{N,N^{\prime }}\left(g\right)(x,t)$. In the matrix form, one can write

(30)${\mathbf{C}}^T\left(x\right)\left(Y-V+D^2{I}_{\kappa }Y-G\right){\mathbf{C}}^T\left(t\right)=0$

where $Y,V,G\in {\mathbb{R}}^{(N+1)\times (N^{\prime }+1)}$, and

$\begin{array}{cc}\hfill \tilde{y}(x,t)& ={\mathbf{C}}^T\left(x\right)Y{\mathbf{C}}^T\left(t\right),\hfill \\ \hfill \tilde{v}(x,t)& ={\mathbf{C}}^T\left(x\right)V{\mathbf{C}}^T\left(t\right),\hfill \\ \hfill \tilde{g}(x,t)& ={\mathbf{C}}^T\left(x\right)G{\mathbf{C}}^T\left(t\right).\hfill \end{array}$

The suggested algorithm utilizes the collocation approach as its foundation. Note that the collocation method selects the solution that satisfies the equation at the collocation points. Equivalently, in this approach, the residual function $r(x,t)$ must approach zero at the collocation points. So, we obtain the residual function

(31)$r(x,t)=\tilde{y}(x,t)-\tilde{v}(x,t)+{\mathcal{I}}_0^{\kappa }\left({(-1)}^{m/2}{\mathcal{D}}_x^m\left(\tilde{y}\right)-g\right)(x,t),$

and set the Lobatto nodes as collocation points.

Thanks to the matrix representation of FIO and derivative operator, and using the Lobatto grid, the following linear system can be obtained:

(32)$\tilde{y}(x_n^{*},x_{n^{\prime }}^{*})-\tilde{v}(x_n^{*},,x_{n^{\prime }}^{*})+{\mathcal{I}}_0^{\kappa }\left({(-1)}^{m/2}{\mathcal{D}}_x^m\left(\tilde{y}\right)-g\right)(x_n^{*},x_{n^{\prime }}^{*})=0,1\le n^{\prime }\le N^{\prime },1\le n\le N.$

Consequently, we have the linear system $AY=G$.

In considering the endpoints as collocation points in the Lobatto grid, only two conditions in the fourth-order FDWE remain unapplied. To use these two conditions, we make the following substitutions:

$\begin{array}{cc}\hfill {\left[AY\right]}_{2,n^{\prime }}={\left[{\mathbf{C}}^T\left(0\right){D^2}^{T}Y\right]}_{n^{\prime }},& {\left[G\right]}_{2,n^{\prime }}={\left[{\Psi }_2\right]}_{n^{\prime }},\hfill \\ \hfill {\left[AY\right]}_{N^{\prime }-1,n^{\prime }}={\left[{\mathbf{C}}^T\left(b\right){D^2}^{T}Y\right]}_{n^{\prime }},& {\left[G\right]}_{N^{\prime }-1,n^{\prime }}={\left[{\Psi }_3\right]}_{n^{\prime }},\hfill \end{array}$

where ${\Psi }_2,{\Psi }_3\in {\mathbb{R}}^{N^{\prime }+1}$ are obtained by

${\psi }_i\approx {\mathcal{Q}}_{N^{\prime }}\left({\psi }_i\right)\left(t\right)=\sum _{n^{\prime }=1}^{N^{\prime }+1}{[\Psi ]}_{n^{\prime }}{C}_{n^{\prime }}\left(t\right):={\Psi }_i^T\mathbf{C}\left(t\right),fori=2,3.$

By converting the matrices A, Y, and G into vectors $\overline{A}$, $\overline{Y}$, and $\overline{G}$, respectively, the following system can be obtained:

(33)$\overline{A}\overline{Y}=\overline{G}.$

Solving this system using the Linsolve command of Matlab (version 2022) leads to finding the unknowns $\left\{y_{n,n^{\prime }}\right\}$, for $1\le n^{\prime }\le N^{\prime }+1$, $1\le n\le N+1$.

Convergence Analysis

The boundedness of the FIO in $L^p[0,b]$ was presented in [30] as follows:

(34)$\parallel {\mathcal{I}}_0^{\kappa }{\left(y\right)\parallel }_p\le K{\parallel y\parallel }_p,\mathrm{with}K={\displaystyle \frac{b^{\kappa }}{\Gamma (\kappa +1)}}.$

Theorem 1. 

Let $1\le N\in \mathbb{N}$. The proposed method converges, and the error satisfies

(35) $\parallel e_N{\parallel }_2\le C\Lambda h^{\nu }{N}^{m-\nu }{\left|y\right|}_{{\mathcal{H}}^{\nu ,N}},$

where C and Λ are constants, $\Lambda =max\left\{\right|v{|}_{{\mathcal{H}}^{\nu ,N}},\left|g{|}_{{\mathcal{H}}^{\nu ,N}}\right\}$, and $h=max\{b,T\}$.

Proof. 

Subtracting Equation (29) from (27) leads to

(36)$e_N-(v-\tilde{v})+{\mathcal{I}}_0^{\kappa }\left({(-1)}^{m/2}{\mathcal{D}}_x^m\left(e_N\right)-(g-\tilde{g})\right)=0,$

where $e_N=y-\tilde{y}$. Taking the norm of (36) and using the triangular inequality gives rise to

(37)$\parallel e_N{\parallel }_2\le \parallel v-\tilde{v}{\parallel }_2+\parallel {\mathcal{I}}_0^{\kappa }\left({\mathcal{D}}_x^m\left(e_N\right)\right){\parallel }_2+{\parallel {\mathcal{I}}_0^{\kappa }(g-\tilde{g})\parallel }_2.$

Motivated by Lemma 1 and in using $h=max\{b,T\}$, the first norm on the left side can be bounded as follows:

(38)$\parallel v-\tilde{v}{\parallel }_2\le C_0{h}^{\nu }{N}^{-\nu }{\left|v\right|}_{{\mathcal{H}}^{\nu ,N}},$

where $C_0$ is a constant. Similarly, thanks to Lemma 1 and using Equation (34), we have

(39)$\parallel {\mathcal{I}}_0^{\kappa }(g-\tilde{g}){\parallel }_2\le K\parallel g-\tilde{g}{\parallel }_2\le C_1{h}^{\nu }{N}^{-\nu }{\left|g\right|}_{{\mathcal{H}}^{\nu ,N}},$

in which $C_1$ is a constant, and K is determined by (34). In addition to Lemma 1 and Equation (34), we need to use Theorem 8.6 [28] to find the bound of the second norm on the left side of Equation (37), viz.,

(40)$|{\mathcal{I}}_0^{\kappa }\left({\mathcal{D}}_x^m\left(e_N\right)\right){\parallel }_2\le K\parallel {\mathcal{D}}_x^m\left(e_N\right){\parallel }_2\le C_2{h}^{\nu -m}{N}^{m-\nu }{\left|y\right|}_{{\mathcal{H}}^{\nu ,N}},0\le m\le \nu \le N+1.$

Consequently, it is easy to confirm that

(41)$\parallel e_N{\parallel }_2\le C_0{h}^{\nu }{N}^{-\nu }\left(max\left\{\right|v{|}_{{\mathcal{H}}^{\nu ,N}},\left|g{|}_{{\mathcal{H}}^{\nu ,N}}\right\}\right)+C_2{h}^{\nu -m}{N}^{m-\nu }{\left|y\right|}_{{\mathcal{H}}^{\nu ,N}}.$

Considering $\Lambda =max\left\{\right|v{|}_{{\mathcal{H}}^{\nu ,N}},\left|g{|}_{{\mathcal{H}}^{\nu ,N}}\right\}$, one can obtain

(42)$\parallel e_N{\parallel }_2\le C\Lambda h^{\nu }{N}^{m-\nu }{\left|y\right|}_{{\mathcal{H}}^{\nu ,N}}.$

□

4. Numerical Experiments

To show the accuracy of the proposed scheme, we calculated the $L^2$-error and $L^{\infty }$-error. These errors are given by the following equations:

$L^2-error={\left({\int }_0^b{\int }_0^T{\left|e_N\right|}^{2}dtdx\right)}^{1/2},$

and

$L^{\infty }-error=\underset{x\in [0,b],t\in [0,T]}{max}\left|e_N(x,t)\right|.$

Example 1. 

Consider the second-order FDWE as follows [32]:

${}_0{}^C\mathcal{D}_t^{\kappa }\left(y\right)(x,t)={\displaystyle \frac{x^2}{2}}{\mathcal{D}}_x^2\left(y\right)(x,t),0\le t\le 1,0\le x\le 1,1<\kappa \le 2,$

with conditions

$\left\{\begin{array}{cc}y(0,t)=0,\hfill & y(1,t)=sinh\left(t\right),\hfill \\ y(x,0)=x,\hfill & {\mathcal{D}}_t\left(y\right)(x,0)=x^2.\hfill \end{array}\right.$

The analytic solution for this equation is reported in [33], obtained by VIM, as

$y(x,t)=x+x^{2}t{E}_{\kappa ,2}\left(t^{\kappa }\right),$

in which $E_{\kappa ,2}$ denotes the Mittag–Leffler function [30].

To demonstrate the presented method’s ability to solve second-order FDWE, Table 1 is provided. As you can see, the method solves the problem with high accuracy. Table 2 tabulates the method’s convergence. Increasing the value of N (number of bases) reduces error. Also, in this table, we report the CPU time. Figure 1 illustrates the approximate solution and corresponding absolute error, taking $\kappa =2$ and $N=10$. Figure 2 plots the absolute errors obtained with different values of N at $t=1$, where $\kappa$ was selected as $15/8$. It is worth mentioning that when $\kappa \to 2$, we expected the solutions to tend toward the corresponding solution obtained with $\kappa =2$. Figure 3 demonstrates this fact.

Example 2. 

The second example delves into the fourth-order FDWE [20]:

${}_0{}^C\mathcal{D}_t^{\kappa }\left(y\right)(x,t)+{\mathcal{D}}_x^4\left(y\right)(x,t)=e^x\left({\displaystyle \frac{\Gamma (\kappa +4)}{6}}{t}^3+t^{\kappa +3}+t\right),0<t\le 1,0<x\le 1,1<\kappa \le 2,$

with

$y(x,0)=0,{\mathcal{D}}_t\left(y\right)(x,0)=e^x,(initial conditions),$

and

$\left\{\begin{array}{ccc}y(0,t)=t^{\kappa +3}+t,\hfill & y(1,t)=e(t^{\kappa +3}+t),\hfill & \hfill \\ {\mathcal{D}}_x^2\left(y\right)(0,t)=t^{\kappa +3}+t,\hfill & {\mathcal{D}}_x^2\left(y\right)(1,t)=e(t^{\kappa +3}+t),\hfill & (boundary conditions).\hfill \end{array}\right.$

As reported in [20], the exact solution is denoted as $y(x,t)=e^x(t^{\kappa +3}+t)$.

The capability of solving fourth-order FDWEs is showcased in Table 3 through the presented method. We tabulated Table 4 to confirm that the method is convergent. Increasing the number of bases (N) reduces error. Also, in this table, we report the CPU time. Table 5 is reported to compare the presented method with CCNFDM [20]. Our scheme gives better accuracy than CCNFDM, and it requires fewer calculations. Figure 4 illustrates the approximate solution and corresponding absolute error, taking $\kappa =1.8$ and $N=10$. We plot the absolute errors at $t=1$, obtained for different values of N, in Figure 5.

5. Conclusions

Considering the importance and application of the FDWE, mentioned in the introduction section, solving this type of equation can be extremely valuable. In this study, we introduce an innovative and precise algorithm for solving the given equation. The algorithm is derived by reducing the equation to the Volterra integral equation and solving this integral equation with the collocation method. To confirm the convergence of the method, we performed a rigorous analysis, and the numerical examples support our findings.

Based on the experimental observations, the following conclusions can be drawn:

• The presented algorithm demonstrates efficacy in solving the FDWE.

• Our algorithm demonstrates convergence when tackling these equations.

• The method presented offers reduced computational cost by avoiding integration to find coefficients.

• The proposed algorithm is simple to implement, and it also provides good accuracy.

Author Contributions

Conceptualization, H.B.J. and I.D.; methodology, H.B.J.; software, H.B.J. and I.D.; validation, H.B.J. and I.D.; formal analysis, H.B.J. and I.D.; investigation, H.B.J. and I.D.; writing—original draft preparation, H.B.J. and I.D.; writing—review and editing, H.B.J. and I.D.; visualization, H.B.J. and I.D.; supervision, H.B.J.; project administration, H.B.J. and I.D.; funding acquisition, H.B.J. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The writers state that they have no known personal relationships or competing financial interests that could have appeared to affect the work reported in this work.

Footnotes

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

Figures and Tables

Figure 1. Plots of the approximate solution and corresponding error, where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (Example 1).

Figure 2. Errors reported at [Forumla omitted. See PDF.] for different values of N, taking [Forumla omitted. See PDF.] for Example 1.

Figure 3. Approximate solution obtained for different values of [Forumla omitted. See PDF.] at [Forumla omitted. See PDF.] for Example 1.

Figure 4. Plots of the approximate solution and corresponding error, where [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (Example 2).

Figure 5. Errors reported at [Forumla omitted. See PDF.] for different values of N, taking [Forumla omitted. See PDF.] for Example 2.

Table 1

The errors obtained using the presented method for Example 1.

$\mathsf{\kappa }$	9/8	11/8	13/8	15/8	2
$L^2$-error	$1.799\times {10}^{-07}$	$8.693\times {10}^{-08}$	$1.460\times {10}^{-08}$	$1.828\times {10}^{-09}$	$4.331\times {10}^{-15}$
$L^{\infty }$-error	$2.493\times {10}^{-06}$	$2.438\times {10}^{-06}$	$8.292\times {10}^{-07}$	$1.730\times {10}^{-07}$	$5.507\times {10}^{-13}$

Table 2

The $L^2$-error, taking $\kappa =15/8$ for Example 1.

$\mathit{t}\setminus \mathit{N}$	6	8	10	12
$0.1$	$7.189\times {10}^{-07}$	$4.229\times {10}^{-07}$	$4.617\times {10}^{-08}$	$2.081\times {10}^{-08}$
$0.2$	$3.153\times {10}^{-06}$	$8.469\times {10}^{-08}$	$7.727\times {10}^{-08}$	$1.740\times {10}^{-08}$
$0.3$	$1.440\times {10}^{-06}$	$3.683\times {10}^{-07}$	$6.500\times {10}^{-08}$	$3.862\times {10}^{-09}$
$0.4$	$1.455\times {10}^{-06}$	$5.512\times {10}^{-08}$	$1.610\times {10}^{-08}$	$1.078\times {10}^{-08}$
$0.5$	$2.465\times {10}^{-06}$	$3.009\times {10}^{-07}$	$5.719\times {10}^{-08}$	$1.527\times {10}^{-08}$
$0.6$	$1.150\times {10}^{-06}$	$4.268\times {10}^{-08}$	$1.175\times {10}^{-08}$	$7.676\times {10}^{-09}$
$0.7$	$8.585\times {10}^{-07}$	$2.090\times {10}^{-07}$	$3.415\times {10}^{-08}$	$1.919\times {10}^{-09}$
$0.8$	$1.443\times {10}^{-06}$	$3.508\times {10}^{-08}$	$2.764\times {10}^{-08}$	$5.731\times {10}^{-09}$
$0.9$	$8.747\times {10}^{-08}$	$1.162\times {10}^{-08}$	$1.016\times {10}^{-08}$	$3.922\times {10}^{-09}$
$1.0$	$2.740\times {10}^{-08}$	$2.053\times {10}^{-09}$	$2.269\times {10}^{-10}$	$3.882\times {10}^{-11}$
CPU Time	$0.281$	$1.422$	$7.266$	$19.484$

Table 3

The errors obtained using the presented method for Example 2.

$\mathsf{\kappa }$	1.2	1.3	1.5	1.7	1.8
$L^2$-error	$1.563\times {10}^{-08}$	$1.552\times {10}^{-08}$	$1.528\times {10}^{-08}$	$1.499\times {10}^{-08}$	$1.485\times {10}^{-08}$
$L^{\infty }$-error	$2.532\times {10}^{-07}$	$3.035\times {10}^{-07}$	$2.889\times {10}^{-07}$	$1.838\times {10}^{-07}$	$1.193\times {10}^{-07}$

Table 4

The $L^2$-error, taking $\kappa =1.5$ for Example 2.

$\mathit{t}\setminus \mathit{N}$	6	8	10	12
$0.1$	$6.698\times {10}^{-05}$	$8.771\times {10}^{-07}$	$1.898\times {10}^{-07}$	$3.667\times {10}^{-09}$
$0.2$	$2.578\times {10}^{-05}$	$6.383\times {10}^{-07}$	$1.503\times {10}^{-07}$	$2.720\times {10}^{-09}$
$0.3$	$3.728\times {10}^{-05}$	$1.890\times {10}^{-06}$	$7.078\times {10}^{-08}$	$3.717\times {10}^{-09}$
$0.4$	$1.158\times {10}^{-04}$	$4.502\times {10}^{-07}$	$4.619\times {10}^{-08}$	$1.043\times {10}^{-08}$
$0.5$	$1.754\times {10}^{-04}$	$9.287\times {10}^{-07}$	$9.714\times {10}^{-08}$	$1.201\times {10}^{-08}$
$0.6$	$1.731\times {10}^{-04}$	$7.166\times {10}^{-07}$	$3.592\times {10}^{-08}$	$7.929\times {10}^{-09}$
$0.7$	$1.676\times {10}^{-04}$	$1.923\times {10}^{-06}$	$4.265\times {10}^{-08}$	$2.120\times {10}^{-09}$
$0.8$	$2.330\times {10}^{-04}$	$7.509\times {10}^{-07}$	$7.270\times {10}^{-08}$	$1.126\times {10}^{-09}$
$0.9$	$3.776\times {10}^{-04}$	$1.104\times {10}^{-06}$	$6.305\times {10}^{-08}$	$1.028\times {10}^{-09}$
$1.0$	$3.967\times {10}^{-04}$	$1.173\times {10}^{-06}$	$5.836\times {10}^{-08}$	$7.207\times {10}^{-09}$
CPU Time	$0.391$	$0.469$	$0.953$	$3.281$

Table 5

A comparison between the presented scheme and CCNFDM [20] for Example 2.

Our Scheme				CCNFDM [20] with $\mathit{h}=1/500$
	$\mathsf{\kappa }$				$\mathsf{\kappa }$
$\mathit{N}={\mathit{N}}^{\prime }$	$\mathbf{1}.\mathbf{3}$	$\mathbf{1}.\mathbf{5}$	$\mathbf{1}.\mathbf{8}$	$\mathsf{\tau }$	$\mathbf{1}.\mathbf{3}$	$\mathbf{1}.\mathbf{5}$	$\mathbf{1}.\mathbf{8}$
6	$5.859\times {10}^{-04}$	$5.563\times {10}^{-04}$	$5.604\times {10}^{-04}$	$1/10$	$2.25\times {10}^{-03}$	$6.21\times {10}^{-03}$	$2.57\times {10}^{-02}$
8	$2.254\times {10}^{-06}$	$2.480\times {10}^{-06}$	$2.060\times {10}^{-06}$	$1/20$	$6.96\times {10}^{-04}$	$2.21\times {10}^{-03}$	$1.13\times {10}^{-02}$
10	$3.035\times {10}^{-07}$	$2.889\times {10}^{-07}$	$1.193\times {10}^{-07}$	$1/40$	$2.15\times {10}^{-04}$	$7.88\times {10}^{-04}$	$4.99\times {10}^{-03}$

References

1. Arif, M.; Ali, F.; Khan, I.; Nisar, K.S. A time fractional model with non-singular kernel the generalized couette flow of couple stress nanofluid. IEEE Access; 2020; 8, pp. 77378-77395. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.2982028]

2. Chang, A.; Sun, H.; Zheng, C.; Bingqing, L.; Chengpeng, L.; Rui, M.; Yong, Z. A Time Fractional Convection-Diffusion Equation to Model Gas Transport through Heterogeneous Soil and Gas Reservoirs. Physica A; 2018; 502, pp. 356-369. [DOI: https://dx.doi.org/10.1016/j.physa.2018.02.080]

3. Tenreiro, M.; Silva, M.F.; Barbosa, R.S.; Jesus, I.S.; Reis, C.M.; Marcos, M.G.; Galhano, A.F. Some Applications of Fractional Calculus in Engineering. Math. Probl. Eng.; 2010; 2010, 639801. [DOI: https://dx.doi.org/10.1155/2010/639801]

4. Mainardi, F. Fractional Calculus and Waves in Linear Viscoelasticity; Imperial College Press: London, UK, 2010.

5. Maji, S.; Natesan, S. Adaptive-grid technique for the numerical solution of a class of fractional boundary-value-problems. Comput. Methods Differ. Equ.; 2010; 12, pp. 338-349.

6. Lakestani, M.; Dehghan, M. The use of Chebyshev cardinal functions for the solution of a partial differential equation with an unknown time-dependent coefficient subject to an extra measurement. J. Comput. Appl. Math.; 2010; 235, pp. 669-678. [DOI: https://dx.doi.org/10.1016/j.cam.2010.06.020]

7. Asadzadeh, M.; Saray, B.N. On a multiwavelet spectral element method for integral equation of a generalized Cauchy problem. BIT Numer. Math.; 2022; 62, pp. 383-1416. [DOI: https://dx.doi.org/10.1007/s10543-022-00915-1]

8. Shi, L.; Saray, B.N.; Soleymani, F. Sparse wavelet Galerkin method: Application for fractional Pantograph problem. J. Comput. Appl. Math.; 2024; 451, 116081. [DOI: https://dx.doi.org/10.1016/j.cam.2024.116081]

9. Garrappa, R. On some explicit Adams multistep methods for fractional differential equations. J. Comput. Appl. Math.; 2009; 229, pp. 392-399. [DOI: https://dx.doi.org/10.1016/j.cam.2008.04.004]

10. Jian, H.Y.; Huang, T.Z.; Zhao, X.; Zhao, Y.L. A fast second-order accurate difference schemes for time distributed-order and Riesz space fractional diffusion equations. J. Appl. Anal. Comput.; 2019; 9, pp. 1359-1392. [DOI: https://dx.doi.org/10.11948/2156-907X.20180247]

11. Lin, Z.; Wang, D.; Qi, D.; Deng, L. A Petrov–Galerkin finite element-meshfree formulation for multi-dimensional fractional diffusion equations. Comput. Mech.; 2020; 66, pp. 323-350. [DOI: https://dx.doi.org/10.1007/s00466-020-01853-x]

12. Fix, G.J.; Roop, J.P. Least squares finite element solution of a fractional order two-point boundary value problem. Comput. Math. Appl.; 2004; 48, pp. 1017-1033. [DOI: https://dx.doi.org/10.1016/j.camwa.2004.10.003]

13. Zhao, Y.L.; Zhu, P.; Gu, X.M.; Zhao, X.; Jian, H.Y. An implicit integration factor method for a kind of spatial fractional diffusion equations. J. Phys. Conf. Ser.; 2019; 1324, 012030. [DOI: https://dx.doi.org/10.1088/1742-6596/1324/1/012030]

14. Daftardar-Gejji, V.; Jafari, A. Adomian decomposition: A tool for solving a system of fractional differential equations. J. Math. Anal. Appl.; 2005; 301, pp. 508-518. [DOI: https://dx.doi.org/10.1016/j.jmaa.2004.07.039]

15. Benkerrouche, B.; Baleanu, D.; Souid, M.S.; Hakem, A. Boundary value problem for nonlinear fractional differential equations of variable order via Kuratowski MNC technique. Adv. Differ. Equ.; 2021; 2021, 365. [DOI: https://dx.doi.org/10.1186/s13662-021-03520-8]

16. Nigmatullin, R.R. To the theoretical explanation of the universal response. Phys. Status Solidi B Basic Res.; 1984; 123, pp. 739-745. [DOI: https://dx.doi.org/10.1002/pssb.2221230241]

17. Oldham, K.B.; Spanier, J. The Fractional Calculus; Academic Press: New York, NY, USA, 1974.

18. Du, R.; Cao, W.R.; Sun, Z.Z. A compact difference scheme for the fractional diffusion-wave equation. Appl. Math. Model.; 2010; 34, pp. 2998-3007. [DOI: https://dx.doi.org/10.1016/j.apm.2010.01.008]

19. Hu, X.; Zhang, L. A compact finite difference scheme for the fourth-order fractional diffusion-wave system. Comput. Phys. Commun.; 2011; 182, pp. 1645-1650. [DOI: https://dx.doi.org/10.1016/j.cpc.2011.04.013]

20. Hu, X.; Zhang, L. On finite difference methods for fourth-order fractional diffusion-wave and subdiffusion systems. Appl. Math. Comput.; 2012; 218, pp. 5019-5034. [DOI: https://dx.doi.org/10.1016/j.amc.2011.10.069]

21. Darzi, R.; Mohammadzade, B.; Mousavi, S.; Beheshti, R. Sumudu transform method for solving fractional differential equations and fractional diffusionwave equation. J. Math. Comput. Sci.; 2013; 6, pp. 79-84. [DOI: https://dx.doi.org/10.22436/jmcs.06.01.08]

22. Cui, M.R. Convergence analysis of high-order compact alternating direction implicit schemes for the two-dimensional time fractional diffusion equation. Numer. Algorithms; 2013; 62, pp. 383-409. [DOI: https://dx.doi.org/10.1007/s11075-012-9589-3]

23. Chen, J.; Liu, F.; Anh, V.; Shen, S.; Liu, Q.; Liao, C. The analytical solution and numerical solution of the fractional diffusion-wave equation with damping. Appl. Math. Comput.; 2012; 219, pp. 1737-1748. [DOI: https://dx.doi.org/10.1016/j.amc.2012.08.014]

24. Bhrawy, A.H.; Doha, E.H.; Baleanu, D.; Ezz-Eldien, S.S. A spectral tau algorithm based on Jacobi operational matrix for numerical solution of time fractional diffusion-wave equations. J. Comput. Phys.; 2015; 293, pp. 142-156. [DOI: https://dx.doi.org/10.1016/j.jcp.2014.03.039]

25. Afarideh, A.; Dastmalchi Saei, F.; Lakestani, M.; Saray, B.N. Pseudospectral method for solving fractional Sturm-Liouville problem using Chebyshev cardinal functions. Phys. Scr.; 2021; 96, 125267. [DOI: https://dx.doi.org/10.1088/1402-4896/ac3c59]

26. Shahriari, M.; Saray, B.N.; Mohammadalipour, B.; Saeidian, S. Pseudospectral method for solving the fractional one-dimensional Dirac operator using Chebyshev cardinal functions. Phys. Scr.; 2023; 98, 055205. [DOI: https://dx.doi.org/10.1088/1402-4896/acc7d3]

27. Boyd, J.P. Chebyshev and Fourier Spectral Methods; 2nd ed. Dover Publications: Garden City, NY, USA, 2001.

28. Shen, J.; Tang, T.; Wang, L.L. Spectral Methods: Algorithms. Analysis, Applications; Springer: Berlin, Germany, 2011.

29. Afarideh, A.; Dastmalchi Saei, F.; Saray, B.N. Eigenvalue problem with fractional differential operator: Chebyshev cardinal spectral method. J. Math. Model.; 2021; 11, pp. 343-355.

30. Kilbas, A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; Elsevier B. V.: Amsterdam, The Netherlands, 2006.

31. Saray, B.N. Abel’s integral operator: Sparse representation based on multiwavelets. BIT Numer. Math.; 2021; 61, pp. 587-606. [DOI: https://dx.doi.org/10.1007/s10543-020-00832-1]

32. Jafari, H.; Momani, S. Solving fractional diffusion and wave equations by modified homotopy perturbation method. Phys. Lett. A; 2007; 370, pp. 388-396. [DOI: https://dx.doi.org/10.1016/j.physleta.2007.05.118]

33. Momani, S.; Odibat, Z. Comparison between the homotopy perturbation method and the variational iteration method for linear fractional partial differential equations. Comput. Math. Appl.; 2007; 54, pp. 910-919. [DOI: https://dx.doi.org/10.1016/j.camwa.2006.12.037]