An accurate approximation technique based on Schröder operational matrices for numerical treatments of time-fractional nonlinear generalized Kawahara equation

Abstract

This study presents a novel numerical method for solving the nonlinear time-fractional generalized Kawahara equations (NTFGKE) under uniform initial boundary conditions (IBCs). To satisfy the aforementioned IBCs, we construct a class of modified Schröder polynomials (MSPs). The suggested approach makes use of operational matrices (OMs) to compute MSPs’ fractional derivatives (FDs) as well as the ordinary derivatives (ODs). These OMs are integrated alongside the spectral collocation method (SCM) to achieve the computational framework. The proposed approach’s convergence and error analysis are provided. We show three numerical test cases to illustrate the accuracy and usefulness of our methodology. The efficacy of the method is demonstrated by benchmarking the obtained numerical outcomes against established solutions. Tables and graphs demonstrate that the suggested method yields very accurate approximate solutions (ASols). Our method demonstrates significant improvements in accuracy and computational efficiency, making it particularly valuable for complex nonlinear problems. This research contributes to the field by providing a robust tool for accurately solving time-fractional equations, which are crucial in modeling various physical phenomena.

Introduction

Fractional differential equations (FDEs) are utilized to describe significant events and find applications in various scientific disciplines [1–3]. Their solutions are highly imperative. Acquiring analytical solutions for these equations is frequently unavailable, hence necessitating the use of numerical methods. Various methodologies are employed to solve distinct types of FDEs. For instance, [4] used the wavelet method. Operational matrices were utilized in multiple articles, such as [5, 6]. The authors in [7] employed a unified predictor–corrector technique to address certain FDEs. The authors in [8, 9] employed a collocation method to address nonlinear variable-order FDEs. The authors in [10] examined the generalized modified Liouville–Caputo FDEs. The authors of [11] utilized a difference scheme to address specific nonlinear FDEs. The researchers in [12] successfully resolved several boundary value problems related to an FDE. The authors of [13–16] addressed several multiorder FDEs.

The numerical techniques for solving the nonlinear and linear partial differential equations (PDEs) attract many researchers [17–24]. Particularly, the numerical methods for solving time-FDEs had been the subject of numerous studies [25–31]. One of the time-FDEs that plays a significant role in physical applications is the Kawahara equation (KE), characterizing the propagation of wave packets in a dispersive medium. The well-known Korteweg–de Vries equation, which models shallow water waves, was generalized in 1972 by Toshio Kawahara. The solution of this particular equation was the focus of a few contributions. As an illustration, the authors achieved precise traveling wave solutions to the KE in [32]. The time-fractional generalized dissipative KE was explicitly solved by the authors in [33]. Various kinds of KEs were studied for their stability, uniqueness, and existence by some authors [34–37]. To address the various KEs, a number of numerical approaches are employed. The replicating kernel technique was applied in [38] to solve the time-fractional KE. In [39], the authors dealt with the generalized KE by using a conservative difference algorithm. Wave solutions for some KEs were created by the author of [40].

The term Schröder number refers to its emergence in the number theory [41]. The initial set of Schröder numbers consists of the following: 1, 2, 6, 22, 90, 394, 1806, 8558, etc. Ernst Schröder, a German mathematician, was credited for his work on this set. A number of articles focused on these numbers and attempted to prove identities for them. The authors in [42] provided some intriguing features of these numbers, such as determinant and product inequalities. Furthermore, in [43], the author presented several recursive formulas for these integers. In [44], the authors discovered certain mathematical characteristics of these numbers. Additional contributions pertaining to these numbers can be discovered in [45, 46]. The literature lacks any contributions related to Schröder polynomials (SPs) and the numerous formulas they contain, despite the existence of various contributions related to Schröder numbers. The authors in [47] encouraged to study these polynomials further and give new formulas such as the inversion formula and the moments of SPs.

We now have an additional reason to look into SPs and utilize them to solve different kinds of differential equations. In the current manuscript, we will demonstrate that using MSPs as basis functions produces precise solutions for the following NTFGKE [48, 49]:

1.1

$\begin{array}{c}({\mathbb{D}}_{\mathfrak{t}}^{\beta }-{\mathcal{D}}_{\mathfrak{z}}^5+{\mathcal{D}}_{\mathfrak{z}}^3+(\mathcal{Y}(\mathfrak{z},\mathfrak{t})+g_1(\mathfrak{z},\mathfrak{t})){\mathcal{D}}_{\mathfrak{z}}+g_2(\mathfrak{z},\mathfrak{t}))\mathcal{Y}(\mathfrak{z},\mathfrak{t})=g_3(\mathfrak{z},\mathfrak{t}),\\ (\mathfrak{z},\mathfrak{t})\in [0,1]\times [0,1],0<\beta \le 1,\end{array}$ with the following initial and boundary conditions:

1.2

$\mathcal{Y}(\mathfrak{z},0)=0,\mathcal{Y}(0,\mathfrak{t})=\mathcal{Y}(1,\mathfrak{t})=0,{\mathcal{Y}}_{\mathfrak{z}}(0,\mathfrak{t}))={\mathcal{Y}}_{\mathfrak{z}}(1,\mathfrak{t})=0,{\mathcal{Y}}_{\mathfrak{z}\mathfrak{z}}(1,\mathfrak{t})=0.$ Here, the functions $g_k(\mathfrak{z},\mathfrak{t})$ are given functions for $k=1,2,3$, and the operator ${\mathbb{D}}_{\mathfrak{t}}^{\beta }$ stands for the fractional derivative interpreted in the sense of Liouville–Caputo. The mathematical framework governed by this equation characterizes the dynamics of nonlinear wave propagation in dispersive media. This work introduces an expanded variant of the standard Kawahara equation in the form [50]

1.3

${\mathbb{D}}_{\mathfrak{t}}^{\beta }y(\mathfrak{z},\mathfrak{t})+\gamma {\mathcal{D}}_{\mathfrak{z}}^{5}y(\mathfrak{z},\mathfrak{t})+\beta {\mathcal{D}}_{\mathfrak{z}}^{3}y(\mathfrak{z},\mathfrak{t})+\beta y(\mathfrak{z},\mathfrak{t}){\mathcal{D}}_{\mathfrak{z}}y(\mathfrak{z},\mathfrak{t})=0,$ through the inclusion of a time-fractional derivative term, which captures memory-dependent dynamics and nonlocal spatial dependencies within the system. The Liouville–Caputo derivative was chosen because of its stability and consistency in numerical implementations. It often yields better convergence properties in numerical solutions (NSols) of fractional differential equations. Additionally, it can accurately model systems with nonlocal effects and time-dependent behaviors, providing a flexible framework for capturing the dynamics of complex systems.

The NTFGKE has garnered a lot of interest in nonlinear research and applied mathematics due to its rich dynamics. To understand the behavior of nonlinear waves in physical systems, it is essential to understand the NTFGKE and its solutions.

We introduce a novel Galerkin OM for FDs and ODs in the Liouville–Caputo sense. Its creation was motivated by a novel approach to solving the problem (1.1) with the conditions (1.2). Particularly designed for MSPs basis vectors are these OMs. By leveraging these operational matrices, we have developed a robust computational framework that enables accurate NSols to be systematically derived using the SCM. As far as we are aware, this study pioneers the use of MSPs in addressing a diverse spectrum of NTFGKE (1.1), establishing a novel computational framework for such systems.

To our knowledge, we are unaware of any previous work that makes use of the OD and FD OMs presented in this study. This gap motivates us to develop and apply these matrices, as they offer a promising alternative that could enhance both stability and accuracy in approximations. Furthermore, we believe these OMs may have broader applications in the future, enabling researchers to tackle a wider range of problems effectively, which provides an additional motivation. On top of that, we will prove that using them yields approximations that match the exact ones. The following are some of the presented unique aspects:

• We create novel OMs based on the SPs for the FDs as well as ODs.

• We utilize these OMs for solving NTFGKE in an efficient and accurate manner.

• The provided algorithms’ great precision, which is on par with certain other algorithms, highlights the importance of the new OMs.

Existing numerical methods, such as finite difference and spectral methods, offer advantages but also face limitations in terms of accuracy, stability, and computational efficiency. Finite difference methods can lack precision in nonlinear scenarios, while spectral methods may require extensive computational resources. This study addresses these challenges by introducing MSPs, which provide a flexible framework for approximating solutions to complex nonlinear equations, including NTFGKE. In comparison to previous approaches, the suggested numerical technique for solving NTFGKE has many notable advantages. It successfully deals with complicated nonlinear problems and achieves excellent accuracy by utilizing MSPs and OMs. Better computing efficiency is another advantage of OM integration with the SCM, which allows for faster computing without losing high accuracy. In addition to being flexible in approximating the solutions of a wide range of FDEs, our technique shows flexibility when dealing with different boundary conditions. Our method is better than others for solving time-FDEs because it gives accurate solutions faster and stays reliable even when conditions change, based on thorough performance reviews.

The paper’s outline comes as given next. In Sect. 2, we state the fundamental facts and characteristics of FDs in detail. Section 3 provides a concise overview of the fundamental properties of SPs. In Sects. 4.1 and 4.2, the primary focus is on developing novel OMs for ODs and FDs of MSPs. The aforesaid OMs are generated in order to solve equation (1.1) when IBCs (1.2) are provided. In Sect. 5, we investigate the potential numerical technique for treating the NTFGKE in the SCM using the built OMs. Section 6 evaluates the error estimate of this new approximation. The suggested method is contrasted with several previously published approaches in Sect. 7, and three examples are given to demonstrate its efficacy. The chief findings and conclusions of our investigation are summarized in Sect. 8.

Basic facts of Liouville–Caputo FDs

In order to build upon our suggested method, this section defines key terms and provides essential resources. These resources are essential for constructing a solid framework for dealing with the NTFGKE.

Definition 2.1

Following the Liouville–Caputo framework, the fractional derivative of $v(\mathfrak{X})$ is expressed by [51]:

2.1

The following characteristics are possessed by the Liouville–Caputo operator [51]. The first is the linearity:

2.2

${\mathbb{D}}_{\mathfrak{X}}^{\beta }(d_1{\ell }_1(\mathfrak{X})+d_2{\ell }_2(\mathfrak{X}))=d_1{\mathbb{D}}_{\mathfrak{X}}^{\beta }{\ell }_1(\mathfrak{X})+d_2{\mathcal{D}}_{\mathfrak{X}}^{\beta }{\ell }_2(\mathfrak{X}).$ In addition, for a constant $\mathcal{C}$ and power-like functions, we have

2.3a

${\mathbb{D}}_{\mathfrak{X}}^{\beta }(\mathcal{C})=0,$

2.3b

${\mathbb{D}}_{\mathfrak{X}}^{\beta }{\mathfrak{X}}^p=\{\begin{array}{cc}0,\hfill & \text{if }p<\lceil \beta \rceil ,\text{and}p\in {\mathbb{N}}_0,\hfill \\ \frac{\mathrm{\Gamma }(p+1)}{\mathrm{\Gamma }(p+1-\beta )}{\mathfrak{X}}^{p-\beta },\hfill & \text{if }p\ge \lceil \beta \rceil .\hfill \end{array}$

Remark 2.1

For those who are interested, FDs can be defined in various ways and have extra characteristics listed in [51].

An overview on Schröder polynomials

Presenting the SPs’ essential characteristics is the primary goal of this section. An SP is defined as follows [41]:

3.1

$S_i(\mathfrak{z})=\sum _{k=0}^i{c}_k^{(i)}{\mathfrak{z}}^k,$ where $c_k^{(i)}=\frac{\left(\genfrac{}{}{0ex}{}{2k}{k}\right)\left(\genfrac{}{}{0ex}{}{i+k}{i-k}\right)}{k+1}$.

Also, SPs can be created by using the following three-term recursive formula [47]:

3.2

$\mathfrak{z}{S}_i(\mathfrak{z})=\frac{i-1}{2(2i+1)}{S}_{i-1}(\mathfrak{z})-\frac{1}{2}{S}_i(\mathfrak{z})+\frac{i+2}{2(2i+1)}{S}_{i+1}(\mathfrak{z}),i\ge 1,$ with $S_0(\mathfrak{z})=1$, $S_1(\mathfrak{z})=1+\mathfrak{z}$. Additionally, the inversion formula takes the form [47]:

3.3

${\mathfrak{z}}^k=\sum _{r=0}^k{b}_r^{(k)}{S}_r(\mathfrak{z}),$ where $b_r^{(k)}=\frac{{(-1)}^{k-r}(2r+1)}{(k+r+1)!(k-r)!}k!(k+1)!$.

To create the intended OMs, the next two expansions are required:

Lemma 3.1

In the expansions

3.4

$S_i(\mathfrak{z})=\sum _{k=0}^i{\overline{c}}_k^{(i)}{(1-\mathfrak{z})}^k$and

3.5

${(1-\mathfrak{z})}^m=\sum _{r=0}^m{\overline{b}}_r^{(m)}{S}_r(\mathfrak{z}),m=0,1,2,\dots ,$we can express the coefficients${\overline{c}}_k^{(i)}$and${\overline{b}}_k^{(i)}$$(i\ge 0)$as follows:

3.6

${\overline{c}}_k^{(i)}=\frac{{(-1)}^k}{k+1}\left(\genfrac{}{}{0ex}{}{2k}{k}\right)\left(\genfrac{}{}{0ex}{}{i+k}{i-k}\right){}_{2}F_1\left[\begin{array}{c}k-i,k+i+1\\ k+2\end{array}\right|-1]$and

3.7

${\overline{b}}_r^{(m)}=\frac{\sqrt{\pi }{(-\frac{1}{4})}^r\mathrm{\Gamma }(r+2)\left(\genfrac{}{}{0ex}{}{m}{r}\right)}{\mathrm{\Gamma }(r+\frac{1}{2})}{}_{2}F_1\left[\begin{array}{c}r-m,r+2\\ 2r+2\end{array}\right|-1].$

Proof

Formula (3.4) is proven first. Here is the form of the polynomial $S_i(\mathfrak{z})$:

3.8

$S_i(\mathfrak{z})=\sum _{k=0}^i{\overline{c}}_k^{(i)}{(1-\mathfrak{z})}^k,{\overline{c}}_k^{(i)}=\frac{{(-1)}^k}{k!}{\mathcal{D}}^k{S}_i(1).$ According to formula (3.1), the form of the coefficients ${\overline{c}}_k^{(i)}$ can be written as (3.6), thanks to Mathematica. We now demonstrate (3.5). Using the structure

3.9

${(1-\mathfrak{z})}^m=\sum _{r=0}^m{(-1)}^r\left(\genfrac{}{}{0ex}{}{m}{r}\right){\mathfrak{z}}^r,$ one gets

3.10

${(1-\mathfrak{z})}^m=\sum _{r=0}^m{(-1)}^r\left(\genfrac{}{}{0ex}{}{m}{r}\right)\sum _{j=0}^r{b}_j^{(r)}{S}_j(\mathfrak{z}).$ After algebraic expansion and grouping of similar terms, we get the formula (3.10) is expressed as (3.5). □

It is now useful to think about the two sets of polynomials ${\{{\phi }_j(\mathfrak{t})\}}_{j\ge 0}$ and ${\{{\mathfrak{U}}_j(\mathfrak{z})\}}_{j\ge 0}$:

3.11

${\phi }_j(\mathfrak{t})=\mathfrak{t}{S}_j(\mathfrak{t}),$

3.12

${\mathfrak{U}}_j(\mathfrak{z})={\mathfrak{z}}^2{(1-\mathfrak{z})}^3{S}_j(\mathfrak{z}),$ to satisfy the specified IBCs (1.2).

Note 3.1

Here, it is important to remember that the generalized hypergeometric function is defined as [52]

3.13

${}_{p}F_q\left[\begin{array}{c}{a}_1,a_2,\dots ,a_p\\ b_1,b_2,\dots ,b_q\end{array}\right|\mathfrak{z}]=\sum _{k=0}^{\mathrm{\infty }}\frac{{(a_1)}_k\dots {(a_p)}_k}{{(b_1)}_k\dots {(b_q)}_k}\frac{{\mathfrak{z}}^k}{k!},$ where $b_j\ne 0$, for all $1\le j\le q$.

The OM for FD and ODs

The construction of OM for FD of ${\phi }_j(\mathfrak{t})$

In this part, we will derive the operational matrices for fractional derivatives associated to ${\phi }_j(\mathfrak{t})$. Consequently, the subsequent theorem is established:

Theorem 4.1

The expression of${\mathbb{D}}^{\beta }{\phi }_i(\mathfrak{t})$may be written as

4.1

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\phi }_i(\mathfrak{t})={\mathfrak{t}}^{-\beta }\sum _{j=0}^i{\mathrm{\Delta }}_{i,j}(\beta ){\phi }_j(\mathfrak{t}),$and as a result

4.2

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\mathit{\phi }}_Q(\mathfrak{t})={\mathfrak{t}}^{-\beta }{\mathbf{D}}^{(\beta )}{\mathit{\phi }}_Q(\mathfrak{t}),$where${\mathit{\phi }}_Q(\mathfrak{t})={[{\phi }_0(\mathfrak{t}),{\phi }_1(\mathfrak{t}),\dots ,{\phi }_Q(\mathfrak{t})]}^T$, and${\mathbf{D}}^{(\beta )}=(d_{i,j}(\beta ))$is a matrix of size$(Q+1)\times (Q+1)$given by

4.3

$\left(\begin{array}{cccccc}{\mathrm{\Delta }}_{0,0}(\beta )& 0& \dots & \dots & \dots & 0\\ {\mathrm{\Delta }}_{1,0}(\beta )& {\mathrm{\Delta }}_{1,1}(\beta )& 0& \dots & \dots & 0\\ ⋮& & \ddots & & & ⋮\\ {\mathrm{\Delta }}_{i,0}(\beta )& \dots & {\mathrm{\Delta }}_{i,i}(\beta )& 0& \dots & 0\\ ⋮& & & \ddots & & ⋮\\ ⋮& & & & \ddots & 0\\ {\mathrm{\Delta }}_{Q,0}(\beta )& \dots & \dots & \dots & \dots & {\mathrm{\Delta }}_{Q,Q}(\beta )\end{array}\right),$where

4.4

$d_{i,j}(\beta )=\{\begin{array}{cc}{\mathrm{\Delta }}_{i,j}(\beta ),\hfill & i\ge j,\hfill \\ 0,\hfill & \mathit{\text{otherwise,}}\hfill \end{array}$and

4.5

${\mathrm{\Delta }}_{i,j}(\beta )=\frac{\sqrt{\pi }{4}^{-j}(j+1)\mathrm{\Gamma }(i+j+1)}{\mathrm{\Gamma }(j+\frac{1}{2})\mathrm{\Gamma }(i-j+1)\mathrm{\Gamma }(j-\beta +2)}{}_{3}F_2\left[\begin{array}{c}j+2,j-i,i+j+1\\ 2j+2,j-\beta +2\end{array}\right|1].$

Proof

By applying the two equations (3.1) and (2.3b), it may be inferred that

4.6

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\phi }_i(\mathfrak{t})={\mathfrak{t}}^{1-\beta }\sum _{k=0}^i{\hat{c}}_k^{(i)}{\mathfrak{t}}^k,{\hat{c}}_k^{(i)}=\frac{\mathrm{\Gamma }(k+2)}{\mathrm{\Gamma }(k-\beta +2)}{c}_k^{(i)},$ and with the help of (3.3) we get

4.7

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\phi }_i(\mathfrak{t})={\mathfrak{t}}^{-\beta }\sum _{j=0}^i{\mathrm{\Delta }}_{i,j}(\beta ){\phi }_j(\mathfrak{t}),$ where

4.8

${\mathrm{\Delta }}_{i,j}(\beta )=\sum _{r=0}^{i-j}\frac{\mathrm{\Gamma }(2+j+r)}{\mathrm{\Gamma }(2+j+r-\beta )}{c}_{r+j}^{(i)}{b}_j^{(r+j)},$ which can be written as (4.5), thanks to Mathematica. Additionally, formula (4.1) can be expressed as

4.9

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\phi }_i(\mathfrak{t})={\mathfrak{t}}^{-\beta }[{\mathrm{\Delta }}_{i,0}(\beta ),{\mathrm{\Delta }}_{i,1}(\beta ),\dots ,{\mathrm{\Delta }}_{i,i}(\beta ),0,\dots ,0]{\mathit{\phi }}_Q(\mathfrak{t}),$ and this form leads to (4.2). □

For instance, if $Q=4$, we get

4.10

${\mathbf{D}}^{(\beta )}={\left(\begin{array}{ccccc}\frac{1}{\mathrm{\Gamma }(2-\beta )}& 0& 0& 0& 0\\ \frac{3\beta }{\mathrm{\Gamma }(3-\beta )}& \frac{2}{\mathrm{\Gamma }(3-\beta )}& 0& 0& 0\\ \frac{6\beta (7\beta -3)}{7\mathrm{\Gamma }(4-\beta )}& \frac{64\beta }{7\mathrm{\Gamma }(4-\beta )}& \frac{6}{\mathrm{\Gamma }(4-\beta )}& 0& 0\\ \frac{5\beta (2\beta (7\beta -9)+13)}{7\mathrm{\Gamma }(5-\beta )}& \frac{60\beta (3\beta -1)}{7\mathrm{\Gamma }(5-\beta )}& \frac{75\beta }{2\mathrm{\Gamma }(5-\beta )}& \frac{24}{\mathrm{\Gamma }(5-\beta )}& 0\\ \frac{15\beta (\beta (\beta (7\beta -18)+32)-13)}{7\mathrm{\Gamma }(6-\beta )}& \frac{80\beta (5(\beta -1)\beta +4)}{7\mathrm{\Gamma }(6-\beta )}& \frac{25\beta (11\beta -3)}{2\mathrm{\Gamma }(6-\beta )}& \frac{192\beta }{\mathrm{\Gamma }(6-\beta )}& \frac{120}{\mathrm{\Gamma }(6-\beta )}\end{array}\right)}_{5\times 5}.$

Remark 4.1

We provide a detailed breakdown of the derivation process for the OMs associated with MSPs. Theorem 4.1 establishes the FD of $\phi (\mathfrak{t})$ as given in (4.1). Here, ${\mathrm{\Delta }}_{i,j}(\beta )$ represent the elements of the OM, which capture the influence of different polynomial degrees on the FD. The matrix ${\mathbf{D}}^{(\beta )}$ is structured to be lower triangular, facilitating efficient computations and reflecting the interactions among polynomial terms. Mathematically, ${\mathrm{\Delta }}_{i,j}(\beta )$ is derived from the integral representation of the FD, indicating how changes in one term influence another. Physically, this matrix models dynamic behaviors in systems described by FDEs, capturing essential characteristics of memory and nonlocality. Thus, ${\mathbf{D}}^{(\beta )}$ systematically accounts for the influence of each polynomial term on the behavior of solutions to NTFGKE, crucial for accurately modeling their dynamics.

The construction of OM for ODs of ${\mathfrak{U}}_i(\mathfrak{z})$

To proceed, we determine the OM for the ODs of ${\mathfrak{U}}_i(\mathfrak{z})$. The following two lemmas are necessary to accomplish this:

Lemma 4.1

The following is a representation of the polynomials$S_i(\mathfrak{z}),i\ge 0$:

4.11

$S_i(\mathfrak{z})=(1-\mathfrak{z})\sum _{j=0}^{i-1}{\mu }_{i,j}{S}_j(\mathfrak{z})+S_i(1),$where

4.12

${\mu }_{i,j}=\sum _{r=0}^{i-1-j}{\overline{c}}_{r+j+1}^{(i)}{\overline{b}}_j^{(r+j)}.$

Proof

We have

4.13

$S_i(\mathfrak{z})=\left(\frac{S_i(\mathfrak{z})-S_i(1)}{\mathfrak{z}-1}\right)(\mathfrak{z}-1)+S_i(1).$ Using (3.4) and (3.5) leads to

4.14

$\frac{S_i(\mathfrak{z})-S_i(1)}{\mathfrak{z}-1}=-\sum _{k=0}^{i-1}{\overline{c}}_{k+1}^{(i)}{(1-\mathfrak{z})}^k=-\sum _{k=0}^{i-1}{\overline{c}}_{k+1}^{(i)}(\sum _{r=0}^k{\overline{b}}_r^{(k)}{S}_i(\mathfrak{z})).$ After expanding, collecting similar terms, and performing some algebraic manipulations,

4.15

$\frac{S_i(\mathfrak{z})-S_i(1)}{\mathfrak{z}-1}=-\sum _{j=0}^{i-1}{\mu }_{i,j}{S}_j(\mathfrak{z}),$ in which ${\mu }_{i,j}$ appears as (4.12). Then (4.11) is obtained by substituting (4.15) into (4.13). □

Lemma 4.2

One possible way to express the derivative$\mathcal{D}({\mathfrak{z}}^2{S}_i(\mathfrak{z})),i\ge 0$is as

4.16

$\mathcal{D}({\mathfrak{z}}^2{S}_i(\mathfrak{z}))=\sum _{j=0}^{i-1}{\gamma }_{i,j}({\mathfrak{z}}^2{S}_j(\mathfrak{z}))+2\mathfrak{z},$where

4.17

${\gamma }_{i,j}=\sum _{r=0}^{i-j-1}(r+j+3)c_{r+j+1}^{(i)}{b}_j^{(r+j)}.$

Proof

We can infer from formula (3.1) that

4.18

$\mathcal{D}({\mathfrak{z}}^2{S}_i(\mathfrak{z}))={\mathfrak{z}}^2\sum _{k=0}^i{\tilde{c}}_k^{(i)}{\mathfrak{z}}^{k-1},{\tilde{c}}_k^{(i)}=(k+2)c_k^{(i)},$ then

4.19

$\mathcal{D}({\mathfrak{z}}^2{S}_i(\mathfrak{z}))={\mathfrak{z}}^2\sum _{k=0}^{i-1}{\tilde{c}}_{k+1}^{(i)}{\mathfrak{z}}^k+{\tilde{c}}_0^{(i)}\mathfrak{z}.$ The structure in (4.16) will result from substituting (3.3) to (4.19) and following some algebraic computation. □

Theorem 4.2

The following representation holds for$\mathcal{D}{\mathfrak{U}}_i(\mathfrak{z})$:

4.20

$\mathcal{D}{\mathfrak{U}}_i(\mathfrak{z})=\sum _{j=0}^{i-1}{\nu }_{i,j}{\mathfrak{U}}_j(\mathfrak{z})+{\varrho }_i(\mathfrak{z}),$where

4.21

${\varrho }_i(\mathfrak{z})={(1-\mathfrak{z})}^2\mathfrak{z}(2-(2+3S_i(1))\mathfrak{z}),$

4.22

${\nu }_{i,j}={\gamma }_{i,j}-3{\mu }_{i,j}.$

Proof

We have

4.23

$\mathcal{D}{\mathfrak{U}}_i(\mathfrak{z})={(1-\mathfrak{z})}^3\mathcal{D}({\mathfrak{z}}^2{S}_i(\mathfrak{z}))-3{(1-\mathfrak{z})}^2{\mathfrak{z}}^2{S}_i(\mathfrak{z}).$ Equation (4.20) is obtained by using Lemmas 4.1 and 4.2. □

We can now calculate the needed OM for the vector

4.24

${\mathfrak{U}}_Q(\mathfrak{z})={[{\mathfrak{U}}_0(\mathfrak{z}),{\mathfrak{U}}_1(\mathfrak{z}),\dots ,{\mathfrak{U}}_Q(\mathfrak{z})]}^T,$ using the following corollary:

Corollary 4.1

We can express themth derivative of${\mathfrak{U}}_Q(\mathfrak{z})$in the form

4.25

$\frac{d^m{\mathfrak{U}}_Q(\mathfrak{z})}{d{\mathfrak{z}}^m}={\mathit{H}}^m{\mathfrak{U}}_Q(\mathfrak{z})+{\mathit{\zeta }}_Q^{(m)}(\mathfrak{z}),$with${\mathit{\zeta }}_Q^{(m)}(\mathfrak{z})=\sum _{k=0}^{m-1}{\mathit{H}}^k{\mathit{\varrho }}_Q^{(m-k-1)}(\mathfrak{z})$, where${\mathit{\varrho }}_Q(\mathfrak{z})={[{\varrho }_0(\mathfrak{z}),{\varrho }_1(\mathfrak{z}),\dots ,{\varrho }_Q(\mathfrak{z})]}^T$Also, we have$\mathit{H}={\left(h_{r,s}\right)}_{0\le r,s\le Q}$, where its components are$h_{r,s}={\nu }_{r,s}$if$r>s$, and$h_{r,s}=0$, otherwise.

For instance, if $Q=4$, we get

4.26

$H={\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ 5& 0& 0& 0& 0\\ -\frac{5}{2}& 16& 0& 0& 0\\ 10& -\frac{5}{4}& \frac{52}{3}& 0& 0\\ -\frac{35}{16}& \frac{214}{12}& -\frac{37}{16}& 21& 0\end{array}\right)}_{5\times 5}.$

Note 4.1

It is important to note that ${\mathit{\zeta }}_Q^{(m)}(\mathfrak{z})$ captures contributions from previous derivatives, reflecting the memory aspect inherent in FDEs. This term highlights how the current state depends not only on the present inputs but also on historical data. For example, in processes like in viscoelastic materials or diffusion, the response at any time is influenced by the entire history of deformation or concentration. The OM H effectively models how these past influences propagate through the system. By incorporating ${\mathit{\zeta }}_Q^{(m)}(\mathfrak{z})$, we enhance our model’s ability to simulate complex dynamics governed by the NTFGKE, providing a more accurate representation of phenomena where the behavior is shaped by both current conditions and historical context.

The computational treatment of NTFGKE under IBCs (1.2)

For the purpose of treating the NTFGKE, we are going to present the MSPs collocation operational matrices method (MSPCOMM) below.

In this section, the NSols for the NTFGKE (1.1) in the presence of IBCs (1.2) are computed using the OMs from Sects. 4.1 and 4.2 which are shown in the following steps:

Step 1. Approximation of the solution

To find an approximation for the solution $\mathcal{Y}(\mathfrak{z},\mathfrak{t})$ of NTFGKE, we express it in terms of OMs:

5.1

$\mathcal{Y}(\mathfrak{z},\mathfrak{t})\simeq {\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})=\sum _{r=0}^Q\sum _{s=0}^Q{a}_{r,s}{\mathfrak{U}}_r(\mathfrak{z}){\phi }_s(\mathfrak{t})={\mathfrak{U}}_Q(\mathfrak{z})\mathit{A}{\mathit{\phi }}_Q(\mathfrak{t}).$ Here the matrix $\mathit{A}={(a_{r,s})}_{r,s=0}^Q$ has size $(Q+1)\times (Q+1)$.

Step 2. Application of OMs

The subsequent equations are established using Theorem 4.1 and Corollary 4.1:

1. Time Derivative

5.2

${\mathbb{D}}_{\mathfrak{t}}^{\beta }{\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})=\sum _{r=0}^Q\sum _{s=0}^Q{a}_{r,s}{\mathfrak{U}}_r(\mathfrak{z}){\mathbb{D}}_{\mathfrak{t}}^{\beta }{\phi }_s(\mathfrak{t})={\mathfrak{t}}^{-\beta }{\mathfrak{U}}_Q(\mathfrak{z})\mathit{A}{\mathit{D}}^{(\beta )}{\mathit{\phi }}_Q(\mathfrak{t}).$

2. Spatial Derivative

5.3

${\mathcal{D}}_{\mathfrak{z}}^m{\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})=\sum _{r=0}^Q\sum _{s=0}^Q{a}_{r,s}{\mathcal{D}}_{\mathfrak{z}}^m{\mathfrak{U}}_r(\mathfrak{z}){\phi }_s(\mathfrak{t})=({\mathit{H}}^m{\mathfrak{U}}_Q(\mathfrak{z})+{\mathit{\zeta }}_Q^{(m)}(\mathfrak{z}))\mathit{A}{\mathit{\phi }}_Q(\mathfrak{t}),$ where $m=1,3,5$.

Step 3. Formulation of the residual

By employing (5.1)–(5.3), it becomes possible to write the residual of (1.1) in the form

5.4

${\mathcal{R}}_Q(\mathfrak{z},\mathfrak{t})=({\mathbb{D}}_{\mathfrak{t}}^{\beta }-{\mathcal{D}}_{\mathfrak{z}}^5+{\mathcal{D}}_{\mathfrak{z}}^3+({\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})+g_1(\mathfrak{z},\mathfrak{t})){\mathcal{D}}_{\mathfrak{z}}+g_2(\mathfrak{z},\mathfrak{t})){\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})-g_3(\mathfrak{z},\mathfrak{t}).$

Step 4. Implementation of the SCM

This is how we find the NSol of (1.1) under IBCs (1.2): we use SCM and the set of collocation nodes $\{({\mathfrak{z}}_i,{\mathfrak{t}}_i):(0\le i\le Q)$. These points can be taken as roots $S_{Q+1}(\mathfrak{t})$. The other option is the uniform points given by ${\mathfrak{z}}_i,{\mathfrak{t}}_i=\frac{(i+1)}{Q+2}$.

Step 5. Establishing the system of equations

Our proposed algorithm, MSPCOMM, is based on the aforementioned points, which serve as the foundation for the spectral approximation. This gives the system

5.5

${\mathcal{R}}_Q({\mathfrak{z}}_r,{\mathfrak{t}}_s)=0,\text{for all}r,s=0,1,\dots ,Q.$

Step 6. Solving for the coefficients

In order to determine $a_{r,s}(r,s=0,1,\dots ,Q)$, the system (5.5) must be solved via a suitable numerical solver. These coefficients critically govern the accuracy of the resulting NSols.

Remark 5.1

In Sect. 7, various numerical examples are solved using MSPCOMM. All computations were performed on an Intel(R) Core(TM) i9-10850 CPU (3.60 GHz, 10 cores, 20 logical processors) using Mathematica 13.3.1.0. Algorithm 1 outlines the systematic workflow for solving the system (1.1)–(1.2), emphasizing the procedural steps essential for deriving the NSols.

[See PDF for image]

Algorithm 1

MSPCOMM to find an approximate solution of (1.1)–(1.2)

Error bound and convergence results for the developed algorithm

Convergence analysis and error estimates for the proposed expansion method are the focus of this section. To continue, we introduce the associated space

6.1

${\mathcal{S}}_{\mathcal{Q}}=\text{Span}⟨{\mathfrak{U}}_r(\mathfrak{z}){\phi }_s(\mathfrak{t})|r,s=0,1,\dots ,Q⟩,$ and the error term

6.2

${\mathcal{E}}_Q(\mathfrak{z},\mathfrak{t})=|{\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})-\mathcal{Y}(\mathfrak{z},\mathfrak{t})|.$ The research uses the $L_{\mathrm{\infty }}$ norm error estimation to assess the numerical scheme’s error,

6.3

${\parallel {\mathcal{E}}_Q\parallel }_{\mathrm{\infty }}={\parallel \mathcal{Y}-{\mathcal{Y}}_Q\parallel }_{\mathrm{\infty }}=\underset{(\mathfrak{z},\mathfrak{t})\in J}{max}|\mathcal{Y}(\mathfrak{z},\mathfrak{t})-{\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})|.$ Here, we have set $J=[0,1]\times [0,1]$. The error estimate is addressed in Theorem 6.1.

Theorem 6.1

Let assume that we have$\mathcal{Y}(\mathfrak{z},\mathfrak{t})={\mathfrak{z}}^2{(1-\mathfrak{z})}^3\mathfrak{t}\upsilon (\mathfrak{z},\mathfrak{t})$. Also, suppose${\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})$is given by (5.1) that provides the closest (best) approximation for the actual solution$\mathcal{Y}(\mathfrak{z},\mathfrak{t})$belonging to${\mathcal{S}}_{\mathcal{Q}}$. The following upper bound is then valid:

6.4

${\parallel {\mathcal{E}}_Q\parallel }_{\mathrm{\infty }}\le \frac{C}{4^Q\mathrm{\Gamma }(Q+1/2)},$where$C=C_1+C_2+C_3$and

6.5

$C_1=\underset{(\mathfrak{z},\mathfrak{t})\in J}{max}\left|\frac{{\partial }^{Q+1}\upsilon (\mathfrak{z},\mathfrak{t})}{\partial {\mathfrak{z}}^{Q+1}}\right|,C_2=\underset{(\mathfrak{z},\mathfrak{t})\in J}{max}\left|\frac{{\partial }^{Q+1}\upsilon (\mathfrak{z},\mathfrak{t})}{\partial {\mathfrak{t}}^{Q+1}}\right|,C_3=\underset{(\mathfrak{z},\mathfrak{t})\in J}{max}\left|\frac{{\partial }^{2Q+2}\upsilon (\mathfrak{z},\mathfrak{t})}{\partial {\mathfrak{z}}^{Q+1}\partial {\mathfrak{t}}^{Q+1}}\right|.$

Proof

Let ${\upsilon }_Q(\mathfrak{z},\mathfrak{t})$ be defined as the polynomial that interpolates $\upsilon (\mathfrak{z},\mathfrak{t})$ at the points $({\mathfrak{z}}_r,{\mathfrak{t}}_s)$ where $S_{Q+1}({\mathfrak{t}}_s)=S_{Q+1}({\mathfrak{z}}_r)=0$, for $r,s=0,1,\dots ,Q$. Then, $\upsilon (\mathfrak{z},\mathfrak{t})$ takes the form [53]

6.6

$\begin{array}{rl}\upsilon (\mathfrak{z},\mathfrak{t})& ={\upsilon }_Q(\mathfrak{z},\mathfrak{t})+\frac{{\partial }^{Q+1}\upsilon ({\eta }_{\mathfrak{z}},\mathfrak{t})}{\partial {\mathfrak{z}}^{Q+1}(Q+1)!}{\mathcal{Q}}_1(\mathfrak{z})+\frac{{\partial }^{Q+1}\upsilon (\mathfrak{z},{\eta }_{\mathfrak{t}})}{\partial {\mathfrak{t}}^{Q+1}(Q+1)!}{\mathcal{Q}}_2(\mathfrak{t})\\ & -\frac{{\partial }^{2Q+2}\upsilon ({\eta }_{\mathfrak{z}}^{\prime },{\eta }_{\mathfrak{t}}^{\prime })}{\partial {\mathfrak{z}}^{Q+1}\partial {\mathfrak{t}}^{Q+1}{((Q+1)!)}^2}{\mathcal{Q}}_1(\mathfrak{z}){\mathcal{Q}}_2(\mathfrak{t}),\end{array}$ where ${\mathcal{Q}}_1(\mathfrak{z})={\prod }_{i=0}^Q(\mathfrak{z}-{\mathfrak{z}}_i)$, ${\mathcal{Q}}_2(\mathfrak{t})={\prod }_{j=0}^Q(\mathfrak{t}-{\mathfrak{t}}_j)$, ${\eta }_{\mathfrak{z}},{\eta }_{\mathfrak{t}},{\eta }_{\mathfrak{z}}^{\prime },{\eta }_{\mathfrak{t}}^{\prime }\in [0,1]$.

Consequently, this results in

6.7

$\begin{array}{rl}{\parallel \upsilon -{\upsilon }_Q\parallel }_{\mathrm{\infty }}& \le \frac{C_1}{(Q+1)!c_{Q+1}^{Q+1}}{\parallel S_{Q+1}(\mathfrak{z})\parallel }_{\mathrm{\infty }}+\frac{C_2}{(Q+1)!c_{Q+1}^{Q+1}}{\parallel S_{Q+1}(\mathfrak{t})\parallel }_{\mathrm{\infty }}\\ & +\frac{C_3}{{((Q+1)!)}^2{(c_{Q+1}^{Q+1})}^2}{\parallel S_{Q+1}(\mathfrak{z})\parallel }_{\mathrm{\infty }}{\parallel S_{Q+1}(\mathfrak{t})\parallel }_{\mathrm{\infty }},\end{array}$ where $c_{Q+1}^{Q+1}=\frac{1}{Q+2}\left(\genfrac{}{}{0ex}{}{2Q+2}{Q+1}\right)$ is the leading coefficient of $S_{Q+1}(\mathfrak{z})$. One then obtains

6.8

$c_{Q+1}^{Q+1}=\frac{4^{Q+1}\mathrm{\Gamma }(Q+3/2)}{(Q+1)!(Q+2)\sqrt{\pi }},Q\ge 0.$ Hence

6.9

${\parallel \upsilon -{\upsilon }_Q\parallel }_{\mathrm{\infty }}\le (C_1+C_2)\frac{1}{4^Q\mathrm{\Gamma }(Q+1/2)}+C_3{\left(\frac{1}{4^Q\mathrm{\Gamma }(Q+1/2)}\right)}^2.$ The inequality (6.9) leads to

6.10

${\parallel \upsilon -{\upsilon }_Q\parallel }_{\mathrm{\infty }}\le \frac{C}{4^Q\mathrm{\Gamma }(Q+1/2)},C=C_1+C_2+C_3.$ Now, take into account the estimate $\mathcal{Y}(\mathfrak{z},\mathfrak{t})\simeq Y_Q(\mathfrak{z},\mathfrak{t})=\mathfrak{t}{\mathfrak{z}}^2{(1-\mathfrak{z})}^3{\upsilon }_Q(\mathfrak{z},\mathfrak{t})$, so

6.11

${\parallel \mathcal{Y}-Y_Q\parallel }_{\mathrm{\infty }}\le \frac{108}{3125}{\parallel \mathcal{Y}-{\upsilon }_Q\parallel }_{\mathrm{\infty }}\le \frac{C}{4^Q\mathrm{\Gamma }(Q+1/2)}.$ But ${\mathcal{Y}}_Q(\mathfrak{z},\mathfrak{t})\in {\mathcal{S}}_Q$ is the best approximation to $\mathcal{Y}(\mathfrak{z},\mathfrak{t})$, so

6.12

${\parallel \mathcal{Y}-{\mathcal{Y}}_Q\parallel }_{\mathrm{\infty }}\le {\parallel \mathcal{Y}-w\parallel }_{\mathrm{\infty }},\mathrm{\forall }w\in {\mathcal{S}}_Q,$ therefore,

6.13

${\parallel \mathcal{Y}-{\mathcal{Y}}_Q\parallel }_{\mathrm{\infty }}\le {\parallel \mathcal{Y}-Y_Q\parallel }_{\mathrm{\infty }}\le \frac{1}{4^Q\mathrm{\Gamma }(Q+1/2)}.$ □

The next corollary illustrates the convergence of the resultant errors.

Corollary 6.1

The following estimate is valid:

6.14

${\parallel {\mathcal{E}}_Q\parallel }_{\mathrm{\infty }}=\mathcal{O}\left(\frac{1}{4^Q\mathrm{\Gamma }(Q+1/2)}\right),Q\ge 1.$

The next theorem confirms the stability of errors, specifically in estimating error propagation.

Theorem 6.2

The following inequality is true:

6.15

$|{\mathcal{Y}}_{Q+1}-{\mathcal{Y}}_Q|\lesssim \mathcal{O}\left(\frac{1}{4^Q\mathrm{\Gamma }(Q+1/2)}\right),Q\ge 1.$

Proof

It is not hard to achieve (6.15) given (6.14). □

Computational outcomes and discussions

Here, we intend to evaluate how well our suggested collocation algorithm (MSPCOMM) handles the NTFGKE (1.1) when subjected to IBCs (1.2). To gauge the practicality and precision of our proposed strategy, we solve a number of test examples and perform comparisons. The evaluation framework relies on the absolute error (AE) as well as the maximum absolute error (MAE) derived from the differences between the exact solution (ESol) and the ASol.

The NSols $y_Q(\mathfrak{z},\mathfrak{t})$ were found using MSPCOMM with $Q=0,1,\dots ,12$, and the estimated errors are visualized in Tables 1, 3, and 7. The computational outcomes presented in these tables exhibit remarkable precision and robustness, underscoring the method’s reliability. Comparisons of MSPCOMM with other methods in [38, 49] are provided in Tables 2, 5, 6, 8, and 9. According to these tables, MSPCOMM produces the most precise outcomes when contrasted with the alternative methods. Moreover, as illustrated through Figs. 1, 2, 3, 4a, 5, 6, 7a, and 8a, there is an excellent level of agreement between the exact and ASols in Examples 7.1–7.3. We demonstrate the convergence, stability, and accuracy of the NSols in Examples 7.1 and 7.2 while applying MSPCOMM with a variety of Q. The analysis includes assessing the stability by examining the behavior of the solutions under varying conditions and perturbations. Error estimates in an adaptive log-scale format are presented in Figs. 4b, 7b, and 8b, illustrating the method’s robustness and consistency across different scenarios.

[See PDF for image]

Figure 1

Plots of exact and approximate solutions and error for Problem 7.1 by utilizing $\beta =0.7$, $Q=15$

[See PDF for image]

Figure 2

Plots of exact and approximate solutions and error for Problem 7.1 by utilizing $Q=15$, $\beta =0.8$

[See PDF for image]

Figure 3

Plots of exact and approximate solutions and error for Example 7.1 by utilizing $Q=15$, $\beta =0.9$

[See PDF for image]

Figure 4

Graphics of errors by utilizing a variety of Q and β for Example 7.1.

[See PDF for image]

Figure 5

Plots of exact and approximate solutions and error for Example 7.2 (Case 1) by utilizing $Q=16$, $\beta =0.7$

[See PDF for image]

Figure 6

Plots of exact and approximate solutions and error for Example 7.2 (Case 2) by utilizing $Q=16$, $\beta =0.8$

[See PDF for image]

Figure 7

Graphics of errors by utilizing a variety of Q and β for Example 7.2 (Case 1).

[See PDF for image]

Figure 8

Graphics of errors by utilizing a variety of Q and β for Example 7.2 (Case 2)

Table 1. MAEs for Example 7.1 using a variety of β

β	Q = 1	Q = 3	Q = 6	Q = 9	Q = 12	Q = 15
0.7	2.1E–01	1.2E–04	2.5E–06	3.1E–08	2.2E–10	4.2E–12
CPU time	0.221	0.340	0.471	0.623	0.801	0.990
0.8	2.2E–01	3.4E–04	2.1E–06	1.1E–08	1.5E–10	2.3E–12
CPU time	0.271	0.352	0.491	0.617	0.812	0.992
0.9	2.2E–02	3.1E–04	2.1E–06	2.2E–08	3.1E–10	4.8E–13
CPU time	0.280	0.361	0.488	0.601	0.825	0.999

Table 2. A comparison of AE determined by MSPCOMM and [38, 49] for Example 7.1$(0.1\le \mathfrak{t}\le 0.9)$

$(\mathfrak{z},\mathfrak{t})$	β = 0.7			β = 0.9
	MSPCOMM (Q = 15)	[49]	[38]	MSPCOMM (Q = 15)	[49]	[38]
$(0.1,\mathfrak{t})$	2.51E–13	1.80E–07	6.36E–06	1.19E–14	2.78E–08	4.92E–06
$(0.2,\mathfrak{t})$	3.11E–13	5.05E–07	1.79E–05	1.32E–13	7.81E–08	1.38E–05
$(0.3,\mathfrak{t})$	3.45E–13	7.61E–07	2.70E–05	1.87E–13	1.18E–07	2.09E–05
$(0.4,\mathfrak{t})$	4.32E–13	8.53E–07	3.03E–05	1.23E–13	1.32E–07	2.34E–05
$(0.5,\mathfrak{t})$	2.17E–13	7.71E–07	2.74E–05	1.17E–13	1.19E–07	2.12E–05
$(0.6,\mathfrak{t})$	4.42E–13	5.69E–07	2.02E–05	1.27E–13	8.80E–08	1.56E–05
$(0.7,\mathfrak{t})$	3.65E–13	3.27E–07	1.16E–05	1.13E–13	5.06E–08	8.99E–06
$(0.8,\mathfrak{t})$	7.23E–13	1.26E–07	4.50E–05	1.92E–14	1.96E–08	3.48E–06
$(0.9,\mathfrak{t})$	6.11E–14	2.00E–08	7.12E–07	4.57E–15	3.10E–09	5.51E–07

Example 7.1

Take into account NTFGKE [38, 49] given by

7.1

$\begin{array}{c}({\mathbb{D}}_{\mathfrak{t}}^{\beta }+{\mathcal{D}}_{\mathfrak{z}}^5+{\mathcal{D}}_{\mathfrak{z}}^3+(\mathcal{Y}(\mathfrak{z},\mathfrak{t})+{\mathfrak{z}}^2\mathfrak{t}+1){\mathcal{D}}_{\mathfrak{z}}+(\mathfrak{z}-\mathfrak{t}))\mathcal{Y}(\mathfrak{z},\mathfrak{t})=g(\mathfrak{z},\mathfrak{t}),\\ (\mathfrak{z},\mathfrak{t})\in [0,1]\times [0,1],\end{array}$ which is controlled by (1.2). The function $g(\mathfrak{z},\mathfrak{t})$ is selected in such a way that $\mathcal{Y}(\mathfrak{z},\mathfrak{t})={\mathfrak{t}}^{\beta +1}{\mathfrak{z}}^2(\frac{{\mathfrak{z}}^3}{6}-\frac{{\mathfrak{z}}^2}{2}+\frac{\mathfrak{z}}{2}-\frac{1}{6})$.

Table 1 and Figs. 1–4 demonstrate that reasonable accuracy is achieved when the suggested technique MSPCOMM is used to get NSols for a variety of Q and β. At $Q=15$, these solutions are in complete agreement with the NSols with an accuracy of 10⁻¹³. AE for MSPCOMM with $Q=15$ compared to the approaches in [38, 49] at $\beta =0.7,0.9$ for $0.1\le \mathfrak{t}\le 0.9$ is also presented in Table 2. This demonstrates that our approach is accurate. Our technique provides efficient performance, considering the CPU time (measured in seconds) given in Table 1.

Remark 7.1

In view of the presented CPU time, our approach has efficient performance. The calculations show that the memory consumption was excellent. For example, the calculated CPU time using $N=15$ is 23.6% slower than $N=12$ and, moreover, requires 22% more RAM memory compared to the $N=12$ calculation. The numerical examples and comparisons provided in our paper highlight the superior accuracy and efficiency of our algorithm, solidifying its potential for solving NTFGKE effectively. When we compared the resource use of our method to those described in [38, 49], we saw that those papers did not give CPU time and memory usage. However, based on our analysis, our approach demonstrates better performance compared to the referenced methods.

Example 7.2

Take into account NTFGKE [38, 49] given by

7.2

$\begin{array}{c}({\mathbb{D}}_{\mathfrak{t}}^{\beta }+{\mathcal{D}}_{\mathfrak{z}}^5+{\mathcal{D}}_{\mathfrak{z}}^3+(\mathcal{Y}(\mathfrak{z},\mathfrak{t})+g_1(\mathfrak{z},\mathfrak{t})){\mathcal{D}}_{\mathfrak{z}}+g_2(\mathfrak{z},\mathfrak{t}))\mathcal{Y}(\mathfrak{z},\mathfrak{t})=g_3(\mathfrak{z},\mathfrak{t}),\\ (\mathfrak{z},\mathfrak{t})\in [0,1]\times [0,1],\end{array}$ which is controlled by (1.2). The function $g_3(\mathfrak{z},\mathfrak{t})$ is selected in such a way that $\mathcal{Y}(\mathfrak{z},\mathfrak{t})=\frac{1}{12}{\mathfrak{t}}^{\beta +1}{\mathfrak{z}}^2({\mathfrak{z}}^4-2{\mathfrak{z}}^3+2\mathfrak{z}-1)$.

The solution to Eq. (7.2) is considered in two cases:

Case 1: $g_1(\mathfrak{z},\mathfrak{t})=1$, $g_2(\mathfrak{z},\mathfrak{t})=0$

The numerical outcomes tabulated in Table 3 for a variety of Q and β are obtained by implementing the suggested method MSPCOMM. Furthermore, the NSols derived with an accuracy of 10⁻¹⁴ at $Q=16$ coincide exactly with the solutions shown in Figs. 5 and 7. When $Q=16$ and the methods in [38, 49] are used with $\beta =0.7,0.9$ and $\mathfrak{t}$ is between 0.1 and 0.9, Table 5 shows how the AE of MSPCOMM compares to those methods. This demonstrates how accurate our approach is.

Table 3. MAEs determined for Example 7.2 (Case 1) and different β

β	Q = 2	Q = 5	Q = 8	Q = 11	Q = 14	Q = 16
0.7	2.2E–03	3.2E–05	3.5E–07	3.5E–10	3.1E–12	4.6E–14
0.8	2.3E–03	3.1E–05	2.7E–07	3.1E–10	5.5E–12	5.5E–14
0.9	4.1E–04	2.8E–06	2.5E–08	2.3E–11	4.5E–13	8.1E–15

Case 2: $g_1(\mathfrak{z},\mathfrak{t})=0$, $g_2(\mathfrak{z},\mathfrak{t})=1$

For a variety of Q and β, the numerical results shown in Table 4 are obtained by applying the proposed method MSPCOMM. Figures 6 and 8 further demonstrate that NSols acquired with an accuracy of 10⁻¹⁵ at $Q=16$ are in complete agreement with the Esol. Table 6 also compares the AE of MSPCOMM with $Q=16$ with the techniques from [38, 49] with $\beta =0.7,0.9$ for $\mathfrak{t}$ values between 0.1 and 0.9. This demonstrates how accurate our approach is.

Table 4. MAEs determined for Example 7.2 (Case 2) and different β

β	Q = 2	Q = 5	Q = 8	Q = 11	Q = 14	Q = 16
0.7	4.2E–03	1.3E–05	2.5E–07	5.3E–10	4.2E–12	5.9E–14
0.8	4.2E–03	4.1E–05	5.5E–07	7.1E–10	3.5E–12	6.2E–14
0.9	2.1E–04	3.2E–06	5.2E–08	4.2E–11	3.1E–13	7.1E–15

Example 7.3

Take into account NTFGKE [38, 49] given by

7.3

$({\mathbb{D}}_{\mathfrak{t}}^{\beta }+{\mathcal{D}}_{\mathfrak{z}}^5+{\mathcal{D}}_{\mathfrak{z}}^3+\mathcal{Y}(\mathfrak{z},\mathfrak{t}){\mathcal{D}}_{\mathfrak{z}}+1)\mathcal{Y}(\mathfrak{z},\mathfrak{t})=g(\mathfrak{z},\mathfrak{t}),(\mathfrak{z},\mathfrak{t})\in [0,1]\times [0,1],$ which is controlled by (1.2). The function $g(\mathfrak{z},\mathfrak{t})$ is selected so that the ESol is $\mathcal{Y}(\mathfrak{z},\mathfrak{t})={\mathfrak{z}}^2({\mathfrak{z}}^4-2{\mathfrak{z}}^3+2\mathfrak{z}-1)sin(2\pi \beta \mathfrak{t})$.

As illustrated in Table 7 and Fig. 9, the suggested approach MSPCOMM yields sufficient accuracy when applied to get NSols for a variety of Q and β. These NSols with an accuracy of 10⁻¹⁶ at $Q=16$ coincide completely with this solution. Furthermore, for $\beta =0.7,0.8,0.9$ and $0.1\le \mathfrak{t}\le 0.9$, Tables 8 and 9 compare the AE of MSPCOMM with $Q=16$ and the approaches described in [49]. This demonstrates how accurate our approach is.

[See PDF for image]

Figure 9

Errors for Example 7.3 at $Q=16$, $\beta =0.7,0.8,0.9$

Remark 7.2

According to Tables 2, 5, 6, 8, and 9, our algorithms outperform those in [38, 49] in terms of accuracy.

Table 5. A comparison of AE determined by MSPCOMM and [38, 49] for Example 7.2 (Case 1) $(0.1\le \mathfrak{t}\le 0.9)$

$(\mathfrak{z},\mathfrak{t})$	β = 0.7			β = 0.9
	MSPCOMM (Q = 16)	[49]	[38]	MSPCOMM (Q = 16)	[49]	[38]
$(0.1,\mathfrak{t})$	1.12E–15	9.87E–08	7.73E–06	1.21E–15	1.52E–08	8.26E–06
$(0.2,\mathfrak{t})$	2.14E–14	3.03E–07	1.57E–05	3.02E–14	4.67E–08	1.75E–05
$(0.3,\mathfrak{t})$	2.12E–14	4.95E–07	1.95E–05	2.51E–14	7.64E–08	2.26E–05
$(0.4,\mathfrak{t})$	1.11E–14	5.97E–07	2.17E–05	4.23E–14	9.22E–08	2.32E–05
$(0.5,\mathfrak{t})$	3.18E–14	5.78E–07	2.09E–05	5.21E–14	8.94E–08	2.00E–05
$(0.6,\mathfrak{t})$	2.12E–14	4.55E–07	1.64E–05	3.12E–14	7.04E–08	1.44E–05
$(0.7,\mathfrak{t})$	3.21E–14	2.78E–07	1.00E–05	2.12E–14	4.30E–08	8.20E–06
$(0.8,\mathfrak{t})$	2.51E–14	1.14E–07	4.10E–06	4.11E–14	1.76E–08	3.18E–06
$(0.9,\mathfrak{t})$	3.31E–15	1.90E–08	6.83E–07	5.10E–15	2.95E–09	5.27E–07

Table 6. A comparison of AE determined by MSPCOMM and [49] for Example 7.2 (Case 2) $(0.1\le \mathfrak{t}\le 0.9)$

$(\mathfrak{z},\mathfrak{t})$	β = 0.7		β = 0.8		β = 0.9
	MSPCOMM (Q = 16)	[49]	MSPCOMM (Q = 16)	[49]	MSPCOMM (Q = 16)	[49]
(0.1,0.1)	3.11E–17	6.70E–10	1.18E–17	3.00E–10	1.16E–17	9.75E–11
(0.2,0.2)	1.14E–16	8.64E–09	2.21E–16	3.97E–09	2.15E–16	1.33E–09
(0.3,0.3)	2.25E–16	5.08E–08	5.22E–16	2.33E–08	1.36E–16	7.76E–09
(0.4,0.4)	4.21E–16	3.90E–07	4.13E–16	1.79E–07	3.22E–16	5.99E–08
(0.5,0.5)	2.21E–16	1.02E–07	5.06E–16	6.87E–10	1.10E–16	3.31E–10
(0.6,0.6)	3.18E–16	4.55E–07	2.26E–16	2.09E–07	1.54E–16	7.04E–08
(0.7,0.7)	3.12E–16	6.95E–08	6.11E–16	3.20E–08	7.10E–16	1.08E–08
(0.8,0.8)	4.11E–16	1.39E–08	5.21E–16	6.41E–09	2.21E–16	2.15E–09
(0.9,0.9)	4.13E–17	1.43E–09	4.27E–17	6.63E–10	4.37E–17	2.25E–10

Table 7. MAEs determined for Example 7.3 using a variety of β and Q

β	Q = 2	Q = 5	Q = 8	Q = 11	Q = 14	Q = 16
0.7	1.2E–04	3.1E–06	1.2E–08	2.5E–10	5.1E–13	5.5E–16
0.8	2.7E–04	1.2E–06	2.1E–08	4.0E–10	3.8E–13	3.4E–16
0.9	1.2E–05	2.1E–07	3.0E–09	3.4E–11	3.5E–14	1.7E–16

Table 8. A comparison of AE determined by MSPCOMM and [49] for Example 7.3 by utilizing $\beta =0.5$

$\mathfrak{z}$	$\mathfrak{t}=0.1$		$\mathfrak{t}=0.3$		$\mathfrak{t}=0.5$
	MSPCOMM (Q = 16)	[49]	MSPCOMM (Q = 16)	[49]	MSPCOMM (Q = 16)	[49]
0.1	2.05E–19	1.04E–15	2.16E–19	1.88E–15	1.16E–19	9.43E–15
0.2	1.11E–18	3.26E–15	1.22E–18	5.91E–15	2.21E–19	2.94E–14
0.3	1.01E–18	5.53E–15	1.54E–18	1.01E–14	1.23E–18	4.99E–14
0.4	2.11E–18	7.13E–15	1.41E–18	1.30E–14	2.13E–18	6.43E–14
0.5	1.12E–18	7.57E–15	2.22E–18	1.38E–14	2.54E–18	6.83E–14
0.6	1.51E–18	6.74E–15	1.13E–18	1.25E–14	2.12E–19	6.08E–14
0.7	1.71E–18	4.88E–15	3.28E–18	9.15E–15	1.45E–18	4.40E–14
0.8	1.55E–18	2.48E–15	2.61E–18	5.04E–15	2.53E–18	2.24E–14
0.9	1.38E–19	3.86E–16	5.53E–19	1.76E–15	2.66E–19	4.00E–15

Table 9. A comparison of AE determined by MSPCOMM and [49] for Example 7.3 by utilizing $\beta =0.5$

$\mathfrak{z}$	$\mathfrak{t}=0.7$		$\mathfrak{t}=0.9$
	MSPCOMM (Q = 16)	[49]	MSPCOMM (Q = 16)	[49]
0.1	2.24E–19	3.49E–15	7.52E–19	3.21E–15
0.2	2.63E–19	1.09E–14	2.51E–19	1.01E–14
0.3	2.52E–18	1.84E–14	2.52E–18	1.73E–14
0.4	2.42E–18	2.37E–14	2.41E–18	2.28E–14
0.5	2.78E–18	2.50E–14	1.34E–18	2.70E–14
0.6	2.63E–18	2.12E–14	2.41E–18	3.71E–14
0.7	2.74E–18	1.12E–14	3.27E–18	9.80E–14
0.8	2.52E–19	2.37E–15	3.11E–19	2.09E–13
0.9	2.37E–19	2.11E–14	4.82E–19	4.79E–13

Conclusion

This study emphasized that the selection of SPs is not only theoretically justified but also practically advantageous for achieving accurate and efficient NSols to NTFGKE. With homogeneous IBCs, this research numerically solved NTFGKE in a new way. The proposed approach appropriately displayed the solutions by making use of MSPs. We employed the SCM for accurate and fast NSol finding. This was accomplished by using OMs for the MSPs’ ODs and FDs. We ensured the reliability of the proposed approach by establishing the convergence of MSPCOMM and carrying out error analysis. Alluding to numerical examples, this study proved that the suggested strategy is accurate and applicable. Comparisons with results from previous studies further demonstrated the algorithm’s accuracy and efficiency. The graphs and tables displayed the tight agreement between the ESol and ASol, demonstrating the excellent accuracy of the proposed algorithm. This method’s efficacy demonstrated its potential to address NTFGKE issues using homogeneous IBCs. This further advances the study of numerical approaches for FDEs. Further research is needed to investigate its efficacy in other contexts and expand its applicability to more complex problems, for example, the numerical approaches to solve multiterm variable-order time-fractional diffusion-wave equations, variable-order fractional nonlinear cable equations, and the two-dimensional variable-order fractional Rayleigh–Stokes model.

Author contributions

Conceptualization: H.M.A. and M.I.; methodology: H.M.A. and M.I.; software: H.M.A.; validation: H.M.A. and M.I.; formal analysis: H.M.A. and M.I.; investigation: H.M.A. and M.I.; writing-original draft preparation: H.M.A. and M.I.; writing-review and editing: H.M.A. and M.I. All authors have and agreed to the published version of the manuscript.

Data Availability

All data generated or analyzed during this study are included in this article.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Caponetto, R.; Giovanni, D.; Luigi, F.; Petras, I. Fractional Order Systems: Modeling and Control Applications; 2010; World Scientific 72

2. Oldham, K.B. Fractional differential equations in electrochemistry. Adv. Eng. Softw.; 2010; 41, 1 pp. 9-12. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2008.12.012]

3. Martines-Arano, H.; Palacios-Barreto, S.; Castillo-Cruz, J.; Meda-Campaña, J.A.; García-Pérez, B.E.; Torres-Torres, C. Fractional photodamage triggered by chaotic attractors in human lung epithelial cancer cells. Int. J. Therm. Sci.; 2022; 181, [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2022.107734] 107734.

4. Wu, M.; Zhang, J.; Huang, Z.; Li, X.; Dong, Y. Numerical solutions of wavelet neural networks for fractional differential equations. Math. Methods Appl. Sci.; 2023; 46, 3 pp. 3031-3044.4538241 [DOI: https://dx.doi.org/10.1002/mma.7449]

5. Ahmed, H.M. A new first finite class of classical orthogonal polynomials operational matrices: an application for solving fractional differential equations. Contemp. Math.; 2023; 4, 4 pp. 974-994. [DOI: https://dx.doi.org/10.37256/cm.4420232716]

6. Pourbabaee, M.; Saadatmandi, A. A new operational matrix based on Müntz–Legendre polynomials for solving distributed order fractional differential equations. Math. Comput. Simul.; 2022; 194, pp. 210-235. [DOI: https://dx.doi.org/10.1016/j.matcom.2021.11.023]

7. Wu, G.C.; Kong, H.; Luo, M.; Fu, H.; Huang, L.L. Unified predictor–corrector method for fractional differential equations with general kernel functions. Fract. Calc. Appl. Anal.; 2022; 25, 2 pp. 648-667.4437296 [DOI: https://dx.doi.org/10.1007/s13540-022-00029-z]

8. Ahmed, H.M. Enhanced shifted Jacobi operational matrices of derivatives: spectral algorithm for solving multiterm variable-order fractional differential equations. Bound. Value Probl.; 2023; 2023, 4667280 [DOI: https://dx.doi.org/10.1186/s13661-023-01796-1] 108.

9. Zheng, X. Numerical approximation for a nonlinear variable-order fractional differential equation via a collocation method. Math. Comput. Simul.; 2022; 195, pp. 107-118.4372811 [DOI: https://dx.doi.org/10.1016/j.matcom.2022.01.005]

10. Talib, I.; Bohner, M. Numerical study of generalized modified Caputo fractional differential equations. Int. J. Comput. Math.; 2023; 100, 1 pp. 153-176.4527208 [DOI: https://dx.doi.org/10.1080/00207160.2022.2090836]

11. Jiang, X.; Wang, J.; Wang, W.; Zhang, H. A predictor–corrector compact difference scheme for a nonlinear fractional differential equation. Fractal Fract.; 2023; 7, 7 521. [DOI: https://dx.doi.org/10.3390/fractalfract7070521]

12. Sevinik Adigüzel, R., Aksoy, Ü., Karapinar, E., Erhan, İ.M.: On the solution of a boundary value problem associated with a fractional differential equation. Math. Methods Appl. Sci. (2020)

13. Izadi, M.; Sene, N.; Adel, W.; El-Mesady, A. The Layla and Majnun mathematical model of fractional order: stability analysis and numerical study. Results Phys.; 2023; 51, [DOI: https://dx.doi.org/10.1016/j.rinp.2023.106650] 106650.

14. Izadi, M.; Afshar, M. Solving the Basset equation via Chebyshev collocation and LDG methods. J. Math. Model.; 2021; 9, pp. 61-79.4185231

15. Ahmed, H.M. Enhanced shifted Jacobi operational matrices of integrals: spectral algorithm for solving some types of ordinary and fractional differential equations. Bound. Value Probl.; 2024; 2024, 14757553 [DOI: https://dx.doi.org/10.1186/s13661-024-01880-0] 75.

16. Youssri, Y.H.; Abd-Elhameed, W.M.; Ahmed, H.M. New fractional derivative expression of the shifted third-kind Chebyshev polynomials: application to a type of nonlinear fractional pantograph differential equations. J. Funct. Spaces; 2022; 2022, 14410782 3966135.

17. Izadi, M.; Srivastava, H.M. Numerical treatments of nonlinear Burgers–Fisher equation via a combined approximation technique. Kuwait J. Sci.; 2024; 51, 4721601 [DOI: https://dx.doi.org/10.1016/j.kjs.2023.12.003] 100163.

18. Izadi, M.; Ahmad, H.; Srivastava, H.M. Numerical computations of time-dependent auto-catalytic glycolysis chemical reaction–diffusion system. MATCH Commun. Math. Comput. Chem.; 2025; 93, pp. 69-97.4914358 [DOI: https://dx.doi.org/10.46793/match.93-1.069I]

19. Ansari, K.J.; Izadi, M.; Noeiaghdam, S. Enhancing the accuracy and efficiency of two uniformly convergent numerical solvers for singularly perturbed parabolic convection–diffusion–reaction problems with two small parameters. Demonstr. Math.; 2024; 57, 4718211 [DOI: https://dx.doi.org/10.1515/dema-2023-0144] 20230144.

20. Izadi, M.; Roul, P. Spectral semi-discretization algorithm for a class of nonlinear parabolic PDEs with applications. Appl. Math. Comput.; 2022; 429, 4421635 127226.

21. Ahmed, H.M. New generalized Jacobi Galerkin operational matrices of derivatives: an algorithm for solving multi-term variable-order time-fractional diffusion-wave equations. Fractal Fract.; 2024; 8, 1 68. [DOI: https://dx.doi.org/10.3390/fractalfract8010068]

22. Dehestani, H.; Ordokhani, Y.; Razzaghi, M. Piecewise numerical approach for solving piecewise fractional optimal control and variational problems. Numer. Algorithms; 2025; [DOI: https://dx.doi.org/10.1007/s11075-025-02024-w]

23. Dehestani, H. Performance of Ritz-piecewise Gegenbauer approach for two types of fractional pantograph equations including piecewise fractional derivative. Math. Methods Appl. Sci.; 2025; 48, pp. 6889-6903.4879277 [DOI: https://dx.doi.org/10.1002/mma.10724]

24. Dehestani, H.; Ordokhani, Y.; Razzaghi, M. Ritz-least squares support vector regression technique for the system of fractional Fredholm–Volterra integro-differential equations. J. Appl. Math. Comput.; 2025; 71, pp. 3477-3508.4910859 [DOI: https://dx.doi.org/10.1007/s12190-025-02365-8]

25. Santra, S.; Mohapatra, J. A novel finite difference technique with error estimate for time fractional partial integro-differential equation of Volterra type. J. Comput. Appl. Math.; 2022; 400, 4293712 [DOI: https://dx.doi.org/10.1016/j.cam.2021.113746] 113746.

26. Lin, J.; Bai, J.; Reutskiy, S.; Lu, J. A novel RBF-based meshless method for solving time-fractional transport equations in 2D and 3D arbitrary domains. Eng. Comput.; 2023; 39, 3 pp. 1905-1922. [DOI: https://dx.doi.org/10.1007/s00366-022-01601-0]

27. Ali, U.; Iqbal, A.; Sohail, M.; Abdullah, F.A.; Khan, Z. Compact implicit difference approximation for time-fractional diffusion-wave equation. Alex. Eng. J.; 2022; 61, 5 pp. 4119-4126. [DOI: https://dx.doi.org/10.1016/j.aej.2021.09.005]

28. Saifullah, S.; Ali, A.; Khan, Z.A. Analysis of nonlinear time-fractional Klein–Gordon equation with power law kernel. AIMS Math.; 2022; 7, 4 pp. 5275-5290.4370227 [DOI: https://dx.doi.org/10.3934/math.2022293]

29. Abd-Elhameed, W.M.; Ahmed, H.M. Spectral solutions for the time-fractional heat differential equation through a novel unified sequence of Chebyshev polynomials. AIMS Math.; 2024; 9, 1 pp. 2137-2166.4685397 [DOI: https://dx.doi.org/10.3934/math.2024107]

30. Javed, S.; Malik, S.A. On the solution of multi-term time fractional diffusion-wave equation involving ultra-hyperbolic operator. Phys. Scr.; 2024; 99, 3 [DOI: https://dx.doi.org/10.1088/1402-4896/ad2250] 035004.

31. Khristenko, U.; Wohlmuth, B. Solving time-fractional differential equations via rational approximation. IMA J. Numer. Anal.; 2023; 43, 3 pp. 1263-1290.4597184 [DOI: https://dx.doi.org/10.1093/imanum/drac022]

32. Gasmi, B.; Moussa, A.A.; Mati, Y.; Alhakim, L.A.; Akgül, A. New exact traveling wave solutions to the Kawahara equation using the $tanh(\xi )$ expansion method. Int. J. Appl. Comput. Math.; 2023; 9, 5 98. [DOI: https://dx.doi.org/10.1007/s40819-023-01568-6]

33. Alquran, M., Ali, M., Alshboul, O.: Explicit solutions to the time-fractional generalized dissipative Kawahara equation. J. Ocean Eng. Sci. (2022)

34. Shuangping, T.; Shuangbin, C. Existence and uniqueness of solution to nonlinear Kawahara equation. Chin. Ann. Math., Ser. A; 2002; 23, 2 pp. 221-228.1900639

35. Cui, S.B.; Deng, D.G.; Tao, S.P. Global existence of solutions for the Cauchy problem of the Kawahara equation with $L^2$ initial data. Acta Math. Sin. Engl. Ser.; 2006; 22, 5 pp. 1457-1466.2251404 [DOI: https://dx.doi.org/10.1007/s10114-005-0710-6]

36. Haragus, M.; Lombardi, E.; Scheel, A. Spectral stability of wave trains in the Kawahara equation. J. Math. Fluid Mech.; 2006; 8, pp. 482-509.2286731 [DOI: https://dx.doi.org/10.1007/s00021-005-0185-3]

37. Kabakouala, A.; Molinet, L. On the stability of the solitary waves to the (generalized) Kawahara equation. J. Math. Anal. Appl.; 2018; 457, 1 pp. 478-497.3702717 [DOI: https://dx.doi.org/10.1016/j.jmaa.2017.08.021]

38. Saldır, O.; Sakar, M.G.; Erdogan, F. Numerical solution of time-fractional Kawahara equation using reproducing kernel method with error estimate. Comput. Appl. Math.; 2019; 38, pp. 1-23.4029559 [DOI: https://dx.doi.org/10.1007/s40314-019-0979-1]

39. Chen, H.; He, Y. Conservative compact difference scheme based on the scalar auxiliary variable method for the generalized Kawahara equation. Math. Methods Appl. Sci.; 2023; 46, 4 pp. 4546-4562.4551381 [DOI: https://dx.doi.org/10.1002/mma.8778]

40. Varol, D. Solitary and periodic wave solutions of the space-time fractional extended Kawahara equation. Fractal Fract.; 2023; 7, 7 539. [DOI: https://dx.doi.org/10.3390/fractalfract7070539]

41. Liu, J.C. A supercongruence involving Delannoy numbers and Schröder numbers. J. Number Theory; 2016; 168, pp. 117-127.3515810 [DOI: https://dx.doi.org/10.1016/j.jnt.2016.04.019]

42. Qi, F.; Shi, X.T.; Guo, B.N. Some properties of the Schröder numbers. Indian J. Pure Appl. Math.; 2016; 47, pp. 717-732.3589757 [DOI: https://dx.doi.org/10.1007/s13226-016-0211-6]

43. Qi, F.; Guo, B.N. Some explicit and recursive formulas of the large and little Schröder numbers. Arab J. Math. Sci.; 2017; 23, 2 pp. 141-147.3662973

44. Sun, Z.W. Arithmetic properties of Delannoy numbers and Schröder numbers. J. Number Theory; 2018; 183, pp. 146-171.3715232 [DOI: https://dx.doi.org/10.1016/j.jnt.2017.07.011]

45. Qi, F.; Guo, B.N. Explicit and recursive formulas, integral representations, and properties of the large Schröder numbers. Kragujev. J. Math.; 2017; 41, 1 pp. 121-141. [DOI: https://dx.doi.org/10.5937/KgJMath1701121F]

46. Qi, F.; Shi, X.; Guo, B. Two explicit formulas of the Schröder numbers. Integers; 2016; 16, A23.

47. Abd-Elhameed, W.M.; Amin, A.K. Novel formulas of Schröder polynomials and their related numbers. Mathematics; 2023; 11, 2 468. [DOI: https://dx.doi.org/10.3390/math11020468]

48. Saldır, O.; Sakar, M.; Erdogan, F. Numerical solution of time-fractional Kawahara equation using reproducing kernel method with error estimate. Comput. Appl. Math.; 2019; 38, 198.4029559 [DOI: https://dx.doi.org/10.1007/s40314-019-0979-1]

49. Abd-Elhameed, W.M.; Youssri, Y.H.; Amin, A.K.; Atta, A.G. Eighth-kind Chebyshev polynomials collocation algorithm for the nonlinear time-fractional generalized Kawahara equation. Fractal Fract.; 2023; 7, 9 652. [DOI: https://dx.doi.org/10.3390/fractalfract7090652]

50. Kawahara, T. Oscillatory solitary waves in dispersive media. J. Phys. Soc. Jpn.; 1972; 33, 1 pp. 260-264. [DOI: https://dx.doi.org/10.1143/JPSJ.33.260]

51. Podlubny, I. Fractional Differential Equations; 1999; San Diego, Academic Press:

52. Luke, Y.L. Special Functions and Their Approximations; 1969; New York, Academic Press:

53. Narumi, S. Some formulas in the theory of interpolation of many independent variables. Tohoku Math. J.; 1920; 18, pp. 309-321.