Novel Approach by Shifted Fibonacci Polynomials for Solving the Fractional Burgers Equation

Abstract

This paper analyzes a novel use of the shifted Fibonacci polynomials (SFPs) to treat the time-fractional Burgers equation (TFBE). We first develop the fundamental formulas of these polynomials, which include their power series representation and the inversion formula. We establish other new formulas for the SFPs, including integer and fractional derivatives, in order to design the collocation approach for treating the TFBE. These derivative formulas serve as tools that aid in constructing the operational metrics for the integer and fractional derivatives of the SFPs. We use these matrices to transform the problem and its underlying conditions into a system of nonlinear equations that can be treated numerically. An error analysis is analyzed in detail. We also present three illustrative numerical examples and comparisons to test our proposed algorithm. These results showed that the proposed algorithm is advantageous since highly accurate approximate solutions can be obtained by choosing a few terms of retained modes of SFPs.

1. Introduction

Special functions in general and special sequences of polynomials, in particular, have vital parts in many applied sciences disciplines, including engineering and physics; see, for example, [1,2]. These sequences are essential in solving various differential equations (DEs). For instance, the authors in [3] studied some fractional differential equations (FDEs) using Morgan-Voyce polynomials. The same sequences of polynomials were utilized in [4] to treat an epidemic SIR model. The authors of [5] studied a fractional-order logistic equation using a type of Dickson polynomials. In [6], the authors handled some sinh-Gordon equations using Lucas polynomials. Some other models that certain FDEs describe were presented in [7]. Some authors used orthogonal sequences of the fifth-kind Chebyshev polynomials (CPs) to solve some differential equations (DEs) in [8,9,10]. In contrast, the sixth-kind Chebyshev polynomials were utilized in [11,12]. Some other generalized CPs were used in [13] to treat some FDEs.

Fibonacci polynomials are among the most significant sequences in both theory and practice. Numerous disciplines make use of Fibonacci polynomials and their corresponding numbers. These disciplines range from mathematics and computer science to physics and biology; see, for example, [14]. In numerical analysis, the role of Fibonacci polynomials has increased. For instance, in [15], the authors used these polynomials to solve the two-dimensional Sobolev equation. Some FDEs are treated using Fibonacci polynomials in [16]. Some delay DEs are solved via the Fibonacci wavelets in [17]. For other applications of Fibonacci polynomials, one can refer to [18,19,20].

Several fields within the applied sciences rely heavily on FDEs. They explain many phenomena that standard DEs cannot. Being able to simulate memory and genetic processes is a significant factor. FDEs may model several physiological and biological processes, including tumor development and neuronal behavior (see [21]). Additionally, these equations model electromagnetic phenomena, wave propagation in complicated mediums, and anomalous diffusion (see [22]). In addition, FDEs have been used to represent the complex mechanical response of viscoelastic materials when subjected to stress or strain (see [23]). Many recent articles were interested in solving some fractional differential models that arise in different applied sciences. For example, in [24], some fractional-order epidemic models were investigated. Some solutions for a certain fractional financial model were developed in [25]. Another fractional differential equation that arises in fluid mechanics was investigated in [26]. For some other fractional models and their treatments, one can refer to [27,28,29].

It is often not possible to obtain analytical solutions for FDEs. As a result, the use of different numerical algorithms becomes necessary. Various FDEs were numerically treated using a wide variety of techniques. Some of these methods are as follows: a predictor–corrector difference method in [30], the Adomian decomposition method in [31,32], the operational collocation method in [33,34,35], matrix methods in [36], and the splines method in [37].

Among the essential FDEs is the classical Burgers equation, which uses fractional derivatives to explain various physical phenomena. The fractional Burgers equation is a substantial expansion of this equation. Because standard integer-order differential equations fail to account for systems with unusual diffusion or memory effects, this equation takes on added significance when representing such processes. A traveling wave with a sharpening front may be obtained by solving the Burgers equations, which are nonlinear PDEs. These equations depict the traffic flow models, including nonlinear propagation and diffusion effects. Due to the importance of the Burgers equation, several approaches, including analytical, numerical, and semi-analytical methods, have been utilized to handle the fractional Burgers equation. For example, the authors of [38] followed an approach to the time and time-space fractional Burgers equations. The authors of [39] obtained numerical solutions for such an equation using the homotopy method. Several approaches are followed to treat numerically coupled systems of Burgers equations. The authors in [40,41] applied collocation procedures to treat some types of coupled Burgers equations. Another wavelet approach is followed in [42] using Gegenbauer polynomials. The Haar-Sinc spectral method is followed in [43] to treat the time-fractional Burgers equation. In [44], a spectral approach is followed to obtain approximate solutions of fractional Burgers equations. For some other contributions regarding the types of Burgers equations, one can refer to [45,46,47].

Spectral methods become vital due to their ability to introduce an additive value in the scope of numerical solutions for DEs of all types. These methods have various advantages comparable to those of other standard numerical methods; see [48]. One of these advantages is the high accuracy of the solutions they produce. There is a strong link between spectral methods and special functions because the solutions can be written as suitable special functions, which could include special polynomials. There are three main approaches to spectral methods. There are some limitations to applying the Galerkin method from the choice of the basis functions; see, for example, [49,50,51]. The tau and collocation methods are also extensively used. Some contributions that utilize the tau method can be found in [52,53], while others regarding the collocation method can be found in [54,55,56]. Some contributions for different spectral methods can be found in [57,58,59].

This paper concentrates on introducing types of polynomials related to Fibonacci polynomials, namely, shifted Fibonacci polynomials. We will introduce some of the basic properties of these polynomials. These formulas will help us find expressions for fractional derivatives of these polynomials, allowing us to tackle the fractional Burgers equation. A spectral collocation method is utilized to obtain the desired numerical solutions. To our knowledge, these polynomials were not previously utilized as basis functions for numerical solutions to different DEs. This gives a strong motivation for investigating them theoretically and employing them numerically. We can enumerate the objectives of this article as follows:

• Introducing the shifted Fibonacci polynomials.
• Establishing some new formulas for these polynomials. These formulas will be pivotal in proposing our numerical algorithm.
• Designing a spectral algorithm based on the typical collocation method to obtain new solutions for fractional Burgers equations.
• Investigating the convergence analysis by developing new inequalities regarding the SFPs.
• We will provide numerical examples and comparisons to test our method.

We point out here that the novelty of our contribution in this paper can be summarized as follows:

• Developing new shifted Fibonacci polynomials.
• Constructing theoretical background concerning these polynomials, more precisely, the fundamental formulas of these polynomials, such as their analytic and inversion formulas. In addition, the integer and fractional derivatives of these polynomials are established. These formulas will be the backbone of applying various numerical methods to different DEs.
• Establishing new operational matrices of integer and fractional derivatives for these polynomials. These matrices are considered important tools for treating DEs.

To the best of our knowledge, using these polynomials in numerical analysis is new and has not been utilized. In addition, we expect that the introduced polynomials will open new horizons in using non-orthogonal polynomials in numerical analysis. Of the advantages of the presented approach is that by choosing modified sets of FPs as basis functions, a few terms of the retained modes make it possible to produce approximations with excellent precision. Less calculation is required, and the resulting errors are small. We also note that the presented collocation algorithm for treating the TFBE is new, which motivates us to analyze it.

The paper follows the following structure: Section 2 presents some preliminaries and essential formulas. Some new useful formulas regarding the SFPs are established in Section 3. Section 4 introduces a collocation algorithm to treat the fractional Burgers equation by employing the SFPs’ integer and fractional derivatives. The convergence and error analysis of the shifted Fibonacci expansion are examined in Section 5. Section 6 presents some illustrative examples and comparisons to test the applicability and efficiency of the proposed method. Section 7 presents some concluding remarks and discussions.

2. Fundamentals and Key Formulas

This section concerns some basic properties of Fibonacci polynomials and new types of polynomials called shifted Fibonacci polynomials. In addition, a brief account of the Caputo fractional derivative is displayed.

2.1. An Overview on Fibonacci Polynomials and Their Shifted Ones

This part provides some characteristics of the Fibonacci polynomials and their shifted polynomials.

The Fibonacci polynomials may be generated by the following recursive formula [14]:

$$F_{\mu +2}\left(\zeta \right)=\zeta F_{\mu +1}\left(\zeta \right)+F_{\mu }\left(\zeta \right),F_0\left(\zeta \right)=0,F_1\left(\zeta \right)=1,\mu =0,1,\dots ,$$

(1)

and they can be written as

$$F_{\mu }\left(\zeta \right)=\sum _{r=0}^{⌊{\displaystyle \frac{\mu -1}{2}}⌋}\left({\displaystyle \genfrac{}{}{0pt}{}{\mu -r-1}{r}}\right){\zeta }^{\mu -2r-1},$$

(2)

where $⌊z⌋$ is the greatest integer less than or equal z.

The first few Fibonacci polynomials are:

$$\begin{array}{cc}\hfill & F_0\left(\zeta \right)=0,F_1\left(\zeta \right)=1,F_2\left(\zeta \right)=\zeta ,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & F_3\left(\zeta \right)={\zeta }^2+1,F_4\left(\zeta \right)={\zeta }^3+2\zeta .\hfill \end{array}$$

(4)

Remark 1. 

It is clear from the recursive Formula (1) that $F_{\mu +1}\left(\zeta \right)$ is of degree μ for any positive integer number μ.

Formula (2) can be inverted to give the following expression for any non-negative integer $\rho$:

$${\zeta }^{\rho }=\sum _{r=0}^{⌊{\displaystyle \frac{\rho }{2}}⌋}{\displaystyle \frac{{(-1)}^r}{r!}}\left(\rho -2r+1\right){\left(\rho -r+2\right)}_{r-1}{F}_{\rho -2r+1}\left(\zeta \right),$$

(5)

where ${\left(a\right)}_{\mu }$ represents the Pochhammer function, which is defined as

$${\displaystyle {\left(a\right)}_{\mu }=\frac{\mathsf{\Gamma }(a+\mu )}{\mathsf{\Gamma }\left(a\right)}.}$$

We will introduce polynomials related to Fibonacci polynomials, which we will utilize in this paper. We will define them as

$$F_{\nu }^{*}\left(\zeta \right)=F_{\nu }(2\zeta -1);$$

(6)

therefore, it is convenient to call them the shifted Fibonacci polynomials (SFPs).

It is to be noted that $F_{\nu +1}^{*}\left(\zeta \right)$ is a polynomial of degree $\nu$.

2.2. The Caputo Fractional Derivative

Definition 1. 

The Caputo fractional derivative is defined as [60]:

$${\displaystyle D_z^{\eta }Y\left(z\right)=\frac{1}{\mathsf{\Gamma }(r-\eta )}{\int }_0^z{(z-t)}^{r-\eta -1}{Y}^{\left(r\right)}\left(t\right)dt,\eta >0,z>0,}$$

(7)

$$r-1<\eta <r,r\in \mathbb{N}.$$

In addition, we have

$$\begin{array}{cc}\hfill D_z^{\eta }C& =0,\left(C\mathrm{is}\mathrm{a}\mathrm{constant}\right),\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill D_z^{\eta }{z}^{ϵ}& =\left\{\begin{array}{cc}0,\hfill & ifϵ\in {\mathbb{N}}_{0}andϵ<\lceil \eta \rceil ,\hfill \\ {\displaystyle \frac{ϵ!}{\mathsf{\Gamma }(ϵ+1-\eta )}}{z}^{ϵ-\eta },\hfill & ifϵ\in {\mathbb{N}}_{0}andϵ\ge \lceil \eta \rceil ,\hfill \end{array}\right.\hfill \end{array}$$

(9)

where $\mathbb{N}=\{1,2,\dots \}$ and  ${\mathbb{N}}_0=\{0,1,2,\dots \}$; $\lceil \eta \rceil$ denotes the ceiling function.

In the following section, some formulas regarding the SFPs will be developed. As far as we know, they are novel. In addition, they will be the backbone to derive our proposed theoretical study and the numerical treatment for the fractional Burgers differential equations.

3. Some Novel Formulas Regarding the SFPs

In this section, we prove several new formulae for the SFPs. First, we will present and prove the power form, inversion representations, and the expression for the SFPs’ derivatives. These formulas will be the keys to deriving the integer and fractional operational matrices of derivatives of the SFPs.

The following theorem gives an explicit power form representation of the shifted polynomials $F_{\nu +1}^{*}\left(\zeta \right)$.

Theorem 1. 

Let $\nu \in {\mathbb{Z}}^{\ge 0}$. The series representation for $F_{\nu +1}^{*}\left(\zeta \right)$ has the form

$$F_{\nu +1}^{*}\left(\zeta \right)=\sum _{\theta =0}^{\nu }{B}_{\theta ,\nu }{\zeta }^{\theta },$$

(10)

where

$$B_{\theta ,\nu }={(-1)}^{\nu +\theta }{\left(2\right)}^{\theta }\left({\displaystyle \genfrac{}{}{0pt}{}{\nu }{\theta }}\right){}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1}{2}}(-\nu +\theta ),{\displaystyle \frac{1}{2}}(1-\nu +\theta )\\ -\nu \end{array}|-4\right).$$

(11)

Proof. 

Starting with (2), we can write

$$F_{\nu +1}^{*}\left(\zeta \right)=\sum _{ϵ=0}^{⌊{\displaystyle \frac{\nu }{2}}⌋}\left({\displaystyle \genfrac{}{}{0pt}{}{\nu -ϵ}{ϵ}}\right){(2\zeta -1)}^{\nu -2ϵ},$$

(12)

which can be expressed as follows using the binomial theorem

$$F_{\nu +1}^{*}\left(\zeta \right)=\sum _{ϵ=0}^{⌊{\displaystyle \frac{\nu }{2}}⌋}{(-1)}^{\nu }\left({\displaystyle \genfrac{}{}{0pt}{}{\nu -ϵ}{ϵ}}\right)\sum _{r=0}^{\nu -2ϵ}{(-2)}^r\left({\displaystyle \genfrac{}{}{0pt}{}{\nu -2ϵ}{r}}\right){\zeta }^r.$$

(13)

Formula (13) can be rearranged to take the following alternative form

$$F_{\nu +1}^{*}\left(\zeta \right)=\sum _{\theta =0}^{\nu }\sum _{\rho =0}^{⌊{\displaystyle \frac{\nu }{2}}⌋}{(-1)}^{\nu +\theta }{2}^{\theta }\left({\displaystyle \genfrac{}{}{0pt}{}{\nu -2\rho }{\theta }}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\nu -\rho }{\rho }}\right){\zeta }^{\theta }.$$

(14)

The last formula is equivalent to the form in (10).    □

Now, we will state and prove the inversion formula to the power form representation of the SFPs given in (10). First, the following lemma is needed to derive the inversion formula.

Lemma 1. 

Let $\theta ,\mu \in {\mathbb{Z}}^{\ge 0}$. We have the following two identities:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)={\displaystyle \frac{2(1+\mu -\theta )\theta }{(1+2\mu )(1+2\mu -2\theta )}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\hfill \\ \hfill & +{\displaystyle \frac{1}{1+2\mu }}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta )}{(1+2\mu )(1+2\mu -2\theta )}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right),\hfill \end{array}\end{array}$$

(15)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {}_{2}F_1\left(\begin{array}{c}-\theta ,-2-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)={\displaystyle \frac{(3+2\mu -2\theta )\theta }{2(1+2\mu )(1+\mu -\theta )}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-2-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\hfill \\ \hfill & +{\displaystyle \frac{\theta (-2-2\mu +\theta )}{1+2\mu }}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\hfill \\ \hfill & +{\displaystyle \frac{(1+2\mu -2\theta )(2+2\mu -\theta )}{2(1+2\mu )(1+\mu -\theta )}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right).\hfill \end{array}\end{array}$$

(16)

Proof. 

The proofs of identities (15) and (16) are similar. We will prove identity (15).

Now, let

$$\begin{array}{cc}\hfill S_{\theta ,\mu }=& {}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)-{\displaystyle \frac{2(1+\mu -\theta )\theta }{(1+2\mu )(1+2\mu -2\theta )}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\hfill \\ \hfill & -{\displaystyle \frac{1}{1+2\mu }}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)-{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta )}{(1+2\mu )(1+2\mu -2\theta )}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right).\hfill \end{array}$$

(17)

To prove Formula (15), it is required to show that $S_{\theta ,\mu }=0$.

Based on the well-known definition of the ${}_{2}F_1\left(z\right)$ ([61]), it is easy to note that the following identity holds for any non-negative integer n and any real numbers $A,B,B\ne 0$:

$${}_{2}F_1\left(\begin{array}{c}-n,A\\ B\end{array}|z\right)=\sum _{r=0}^n{\displaystyle \frac{{(-n)}_r{\left(A\right)}_r}{{\left(B\right)}_{r}r!}}{z}^r.$$

The above identity enables one to write $S_{\theta ,\mu }$ in the following form:

$$\begin{array}{cc}\hfill S_{\theta ,\mu }=& \sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}-{\displaystyle \frac{2(1+\mu -\theta )\theta }{(1+2\mu )(1+2\mu -2\theta )}}\sum _{r=0}^{\theta -1}{\displaystyle \frac{{(1-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}\hfill \\ \hfill & -{\displaystyle \frac{1}{1+2\mu }}\sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{1}{2}\right)}_{r}r!{(-4)}^r}}-{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta )}{(1+2\mu )(1+2\mu -2\theta )}}\sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}.\hfill \end{array}$$

(18)

Based on the following simple identities:

$$\begin{array}{cc}\hfill {(-2\mu +\theta )}_r=& {\displaystyle \frac{1+2\mu -\theta -r}{1+2\mu -\theta }}{(-2\mu +\theta -1)}_r,\hfill \\ \hfill {\left({\displaystyle \frac{1}{2}}\right)}_r=& {\displaystyle \frac{1}{1+2r}}{\left({\displaystyle \frac{3}{2}}\right)}_r,\hfill \\ \hfill \theta {(1-\theta )}_r=& (\theta -r){(-\theta )}_r,\hfill \end{array}$$

$S_{\theta ,\mu }$ can be converted into the following expression

$$\begin{array}{cc}\hfill S_{\theta ,\mu }=& \sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}-{\displaystyle \frac{2(1+\mu -\theta )}{(1+2\mu )(1+2\mu -2\theta )}}\sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r(\theta -r){(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}\hfill \\ \hfill & -{\displaystyle \frac{2r+1}{2\mu +1}}\sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}-{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta )}{(1+2\mu )(1+2\mu -2\theta )}}\sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-2\mu +\theta -1)}_r(1+2\mu -\theta -r)}{(1+2\mu -\theta )\left({\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r\right)}},\hfill \end{array}$$

(19)

which can be simplified again into the form

$$\begin{array}{cc}\hfill S_{\theta ,\mu }=& \sum _{r=0}^{\theta }{\displaystyle \frac{{(-\theta )}_r{(-1-2\mu +\theta )}_r}{{\left(\frac{3}{2}\right)}_{r}r!{(-4)}^r}}\left(1-{\displaystyle \frac{1+2r}{1+2\mu }}-{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta -r)}{(1+2\mu )(1+2\mu -2\theta )}}-{\displaystyle \frac{2(1+\mu -\theta )(\theta -r)}{(1+2\mu )(1+2\mu -2\theta )}}\right).\hfill \end{array}$$

(20)

Noting that

$$1-{\displaystyle \frac{1+2r}{1+2\mu }}-{\displaystyle \frac{2(\mu -\theta )(1+2\mu -\theta -r)}{(1+2\mu )(1+2\mu -2\theta )}}-{\displaystyle \frac{2(1+\mu -\theta )(\theta -r)}{(1+2\mu )(1+2\mu -2\theta )}}=0,$$

we therefore conclude that

$$S_{\theta ,\mu }=0.$$

This proves (15).    □

Now, we are in a position to state and prove the following inversion formula of SFPs.

Theorem 2. 

For a positive integer μ, we can write

$${\zeta }^{\mu }=\sum _{\rho =0}^{\mu }{H}_{\rho ,\mu }{F}_{\rho +1}^{*}\left(\zeta \right),$$

(21)

where

$$H_{\rho ,\mu }=\left\{\begin{array}{cc}{\displaystyle \frac{{(-1)}^{\frac{\mu -\rho }{2}}{2}^{-\mu }(1+\rho )\mu !}{\left(\frac{\mu -\rho }{2}\right)!\left(\frac{\mu +\rho }{2}+1\right)!}{}_{2}F_1\left(\begin{array}{c}\frac{1}{2}(-2-\mu -\rho ),\frac{1}{2}(-\mu +\rho )\\ \frac{1}{2}\end{array}|-\frac{1}{4}\right),}\hfill & (\rho +\mu )\mathit{even},\hfill \\ {\displaystyle \frac{{(-1)}^{\frac{1}{2}(-1+\mu -\rho )}{2}^{-\mu }(1+\rho )\mu !}{\left(\frac{1}{2}(\mu -\rho -1)\right)!\left(\frac{1}{2}(\mu +\rho +1)\right)!}{}_{2}F_1\left(\begin{array}{c}\frac{1}{2}(-1-\mu -\rho ),\frac{1}{2}(1-\mu +\rho )\\ \frac{3}{2}\end{array}|-\frac{1}{4}\right),}\hfill & (\rho +\mu )\mathit{odd}.\hfill \end{array}\right.$$

(22)

Proof. 

To prove Formula (21), we are going to prove its alternative formula:

$${\zeta }^{\mu }=\sum _{\theta =0}^{⌊{\displaystyle \frac{\mu }{2}}⌋}{G}_{\theta ,\mu }{F}_{\mu -2\theta }\left(\zeta \right)+\sum _{\theta =0}^{⌊{\displaystyle \frac{\mu -1}{2}}⌋}{R}_{\theta ,\mu }{F}_{\mu -2\theta -1}\left(\zeta \right),$$

(23)

where

$$\begin{array}{cc}\hfill G_{\theta ,\mu }=& {\displaystyle \frac{{(-1)}^{\theta +1}{2}^{-\mu }(-1-\mu +2\theta )\mu !}{\theta !(\mu -\theta +1)!}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill R_{\theta ,\mu }=& {\displaystyle \frac{{(-1)}^{\theta +1}{2}^{-\mu }(-\mu +2\theta )\mu !}{\theta !(\mu -\theta )!}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right).\hfill \end{array}$$

(25)

We are going to prove (23) by induction. It is easy to see that (23) is valid for $\mu =0$. Assuming that Formula (23) holds, to prove the inductive step, we need to demonstrate that the following formula is valid:

$${\zeta }^{\mu +1}=\sum _{\theta =0}^{⌊{\displaystyle \frac{\mu +1}{2}}⌋}{G}_{\theta ,\mu +1}{F}_{\mu -2\theta +1}^{*}\left(\zeta \right)+\sum _{\theta =0}^{⌊{\displaystyle \frac{\mu }{2}}⌋}{R}_{\theta ,\mu +1}{F}_{\mu -2\theta }^{*}\left(\zeta \right).$$

(26)

It is possible to divide the previous formula into the following two formulas:

$$\begin{array}{cc}\hfill {\zeta }^{2\mu }=& \sum _{\theta =0}^{\mu }{G}_{\theta ,2\mu }{F}_{2\mu -2\theta }^{*}\left(\zeta \right)+\sum _{\theta =0}^{\mu }{R}_{\theta ,2\mu }{F}_{2\mu -2\theta -1}^{*}\left(\zeta \right),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\zeta }^{2\mu +1}=& \sum _{\theta =0}^{\mu }{R}_{\theta ,2\mu +1}{F}_{2\mu -2\theta }^{*}\left(\zeta \right)+\sum _{\theta =0}^{\mu }{G}_{\theta ,2\mu +1}{F}_{2\mu -2\theta +1}^{*}\left(\zeta \right).\hfill \end{array}$$

(28)

Similarities exist between the proofs of (27) and (28). Equation (28) will be proved.

Starting from the valid Formula (23), but replacing $\mu$ with $2\mu +1$, we can write

$$\begin{array}{cc}\hfill {\zeta }^{2\mu +1}=& {\displaystyle \frac{1}{2}}\sum _{\theta =0}^{\mu }{G}_{\theta ,2\mu }\left(F_{2\mu -2\theta +1}^{*}\left(\zeta \right)+F_{2\mu -2\theta }^{*}\left(\zeta \right)-F_{2\mu -2\theta -1}^{*}\left(\zeta \right)\right)\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\sum _{\theta =0}^{\mu -1}{R}_{\theta ,2\mu }\left(F_{2\mu -2\theta }^{*}\left(\zeta \right)+F_{2\mu -2\theta -1}^{*}\left(\zeta \right)-F_{2\mu -2\theta -2}^{*}\left(\zeta \right)\right).\hfill \end{array}$$

(29)

Through a series of algebraic calculations, the final formula is transformed into

$${\zeta }^{2\mu +1}=\sum _{\theta =0}^{\mu }{H}_{\theta ,\mu }{F}_{2\mu -2\theta }^{*}\left(\zeta \right)+\sum _{\theta =0}^{\mu }{\overline{H}}_{\theta ,\mu }{F}_{2\mu -2\theta +1}^{*}\left(\zeta \right),$$

(30)

where

$$\begin{array}{cc}\hfill H_{\theta ,\mu }=& {\displaystyle \frac{1}{2}}\left(G_{\theta ,2\mu }+R_{\theta ,2\mu }-R_{\theta -1,2\mu }\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\overline{H}}_{\theta ,\mu }=& {\displaystyle \frac{1}{2}}\left(G_{\theta ,2\mu }-G_{\theta -1,2\mu }+R_{\theta -1,2\mu }\right).\hfill \end{array}$$

(32)

The coefficients $H_{\theta ,\mu }$ and ${\overline{H}}_{\theta ,\mu }$ can be explicitly computed to give

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_{\theta ,\mu }=& \frac{1}{2}{(-1)}^{\theta +1}{2}^{-2\mu }\left(2\mu \right)!\left({\displaystyle \frac{-2\mu +2(\theta -1)}{(1+2\mu -\theta )!(\theta -1)!}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{-1-2\mu +2\theta }{(1+2\mu -\theta )!\theta !}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{-2\mu +2\theta }{(2\mu -\theta )!\theta !}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\right),\hfill \end{array}\end{array}$$

(33)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{H}}_{\theta ,\mu }=& {\displaystyle \frac{1}{2}}{(-1)}^{\theta +1}{2}^{-2\mu }\left(2\mu \right)!\left({\displaystyle \frac{-1-2\mu +2(\theta -1)}{(2+2\mu -\theta )!(\theta -1)!}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-2-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\right.\hfill \\ \hfill & \left.-{\displaystyle \frac{-2\mu +2(\theta -1)}{(1+2\mu -\theta )!(\theta -1)!}}{}_{2}F_1\left(\begin{array}{c}1-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)+{\displaystyle \frac{-1-2\mu +2\theta }{(1+2\mu -\theta )!\theta !}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\right).\hfill \end{array}\end{array}$$

(34)

The application of the two identities of Lemma 1 reduces the coefficients $H_{\theta ,\mu }$ and ${\overline{H}}_{\theta ,\mu }$, given by (33) and (34), respectively, to the following forms:

$$\begin{array}{cc}\hfill H_{\theta ,\mu }=& {\displaystyle \frac{{(-1)}^{\theta +1}{2}^{-1-2\mu }(-1-2\mu +2\theta )(2\mu +1)!}{\theta !(2\mu -\theta +1)!}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-1-2\mu +\theta \\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\overline{H}}_{\theta ,\mu }=& {\displaystyle \frac{{(-1)}^{\theta +1}{4}^{-\mu }(-1-\mu +\theta )(2\mu +1)!}{\theta !(2\mu -\theta +2)!}}{}_{2}F_1\left(\begin{array}{c}-\theta ,-2-2\mu +\theta \\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right).\hfill \end{array}$$

(36)

and this leads to the two identities

$$H_{\theta ,\mu }=R_{\theta ,2\mu +1},{\overline{H}}_{\theta ,\mu }=R_{\theta ,2\mu +1},$$

and this proves Formula (28).    □

Now, the following theorem expresses the derivatives of the SFPs in terms of SFPs. This expression will be the key to establishing the SFPs’ operational matrix of integer derivatives.

Theorem 3. 

Let $\nu ,q\in {\mathbb{Z}}^{+}$ with $\nu \ge q$. We have

$${\displaystyle D^q{F}_{\nu +1}^{*}\left(\zeta \right)=\sum _{ϵ=0}^{\nu -q}{M}_{ϵ,\nu ,q}{F}_{ϵ}^{*}\left(\zeta \right),}$$

(37)

where

$$M_{ϵ,\nu ,q}={\theta }_{ϵ,\nu ,q}{2}^q{(-1)}^{{\displaystyle \frac{\nu -ϵ-q}{2}}}(ϵ+1)\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}(-2+\nu -ϵ+q)}{\frac{1}{2}(\nu -ϵ-q)}}\right){\left({\displaystyle \frac{1}{2}}(4+\nu +ϵ-q)\right)}_{q-1},$$

with

$${\theta }_{ϵ,\nu ,q}=\left\{\begin{array}{cc}1,\hfill & (\nu -ϵ-q)\mathit{even},\hfill \\ 0,\hfill & \mathit{otherwise}.\hfill \end{array}\right.$$

Proof. 

From [62], the $qth$ derivative of $F_{\nu +1}\left(\zeta \right)$ was expressed as combinations of their original polynomials as

$$D^q{F}_{\nu +1}\left(\zeta \right)=\sum _{m=0}^{⌊{\displaystyle \frac{\nu -q}{2}}⌋}{(-1)}^m(1+\nu -2m-q)\left({\displaystyle \genfrac{}{}{0pt}{}{m+q-1}{m}}\right){(\nu -m-q+2)}_{q-1}{F}_{\nu -2m-q+1}\left(\zeta \right).$$

(38)

Formula (38) can be transformed into the following formula:

$$D^q{F}_{\nu +1}\left(\zeta \right)=\sum _{ϵ=0}^{\nu -q}{\theta }_{\nu ,ϵ,q}{(-1)}^{{\displaystyle \frac{\nu -ϵ-q}{2}}}(ϵ+1)\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}(-2+\nu -ϵ+q)}{\frac{1}{2}(\nu -ϵ-q)}}\right){\left({\displaystyle \frac{1}{2}}(4+\nu +ϵ-q)\right)}_{q-1}{F}_{ϵ+1}\left(\zeta \right),$$

(39)

where

$${\theta }_{ϵ,\nu ,q}=\left\{\begin{array}{cc}1,\hfill & (\nu -ϵ-q)\mathrm{even},\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(40)

If we replace $\zeta$ by $(2\zeta -1)$ in (39), then we can obtain the following formula:

$$D^q{F}_{\nu +1}^{*}\left(\zeta \right)=\sum _{ϵ=0}^{\nu -q}{\theta }_{\nu ,ϵ,q}{2}^q{(-1)}^{{\displaystyle \frac{\nu -ϵ-q}{2}}}(ϵ+1)\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}(-2+\nu -ϵ+q)}{\frac{1}{2}(\nu -ϵ-q)}}\right){\left({\displaystyle \frac{1}{2}}(4+\nu +ϵ-q)\right)}_{q-1}{F}_{ϵ+1}^{*}\left(\zeta \right),$$

(41)

which can be written as in (37). This completes the proof of Theorem 3.    □

Remark 2. 

The first and second derivatives of the SFPs can be deduced explicitly as special cases of Formula (37).

Corollary 1. 

${\displaystyle \frac{d{F}_{\nu +1}^{*}\left(\zeta \right)}{d\zeta }}$ has the form

$${\displaystyle \frac{d{F}_{\nu +1}^{*}\left(\zeta \right)}{d\zeta }}=\sum _{ϵ=0}^{\nu -1}{\lambda }_{\nu ,ϵ}{F}_{ϵ+1}^{*}\left(\zeta \right),\nu \ge 1,$$

(42)

where

$${\lambda }_{\nu ,ϵ}=2(ϵ+1){(-1)}^{{\displaystyle \frac{1}{2}}(\nu -ϵ-1)}{\theta }_{ϵ,\nu ,1}.$$

(43)

Proof. 

Formula (42) can be easily obtained by setting $q=1$ in Theorem 3.    □

Corollary 2. 

${\displaystyle \frac{d^2{F}_{\nu +1}^{*}\left(\zeta \right)}{d{\zeta }^2}}$ has the form

$${\displaystyle \frac{d^2{F}_{\nu +1}^{*}\left(\zeta \right)}{d{\zeta }^2}}=\sum _{ϵ=0}^{\nu -2}{\beta }_{\nu ,ϵ}{F}_{ϵ+1}^{*}\left(\zeta \right),\nu \ge 2,$$

(44)

where

$${\beta }_{\nu ,ϵ}=-(ϵ+1){(-1)}^{{\displaystyle \frac{\nu -ϵ}{2}}}(\nu -ϵ)(\nu +ϵ+2){\theta }_{ϵ,\nu ,2}.$$

(45)

Proof. 

Formula (44) can be easily obtained by setting $q=2$ in Theorem 3.    □

The following corollary gives the operational matrices of the first and second derivatives of the SFPs.

Corollary 3. 

If we consider the vector ${\mathcal{F}}^{*}\left(\zeta \right)={[F_1^{*}\left(\zeta \right),F_2^{*}\left(\zeta \right),\dots ,F_{\mathcal{N}+1}^{*}\left(\zeta \right)]}^T$, then the first and second derivatives of the vector ${\mathcal{F}}^{*}\left(\zeta \right)$ can be written in matrix form as

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{\mathcal{F}}^{*}\left(\zeta \right)}{d\zeta }}=\mathit{\lambda }{\mathcal{F}}^{*}\left(\zeta \right),\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d^2{\mathcal{F}}^{*}\left(\zeta \right)}{d^2\zeta }}=\mathit{\beta }{\mathcal{F}}^{*}\left(\zeta \right),\hfill \end{array}$$

(47)

where $\mathit{\lambda }=\left({\lambda }_{\nu ,ϵ}\right)$ and $\mathit{\beta }=\left({\beta }_{\nu ,ϵ}\right)$ are the operational matrices of derivatives of order $(\mathcal{N}+1)\times (\mathcal{N}+1)$ whose entries are given, respectively, as

$$\begin{array}{cc}& {\lambda }_{\nu ,ϵ}=\left\{\begin{array}{cc}2(ϵ+1){(-1)}^{{\displaystyle \frac{1}{2}}(\nu -ϵ-1)},\hfill & if\nu >ϵ,(\nu -ϵ)odd,\hfill \\ 0,\hfill & otherwise,\hfill \end{array}\right.\hfill \\ & {\beta }_{\nu ,ϵ}=\left\{\begin{array}{cc}-(ϵ+1){(-1)}^{{\displaystyle \frac{\nu -ϵ}{2}}}(\nu -ϵ)(\nu +ϵ+2),\hfill & if\nu >ϵ+1,(\nu -ϵ)even,\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.\hfill \end{array}$$

(48)

Proof. 

Formulas (46) and (47) are immediate results of Formulas (42) and (44), respectively.    □

The following theorem exhibits an explicit expression for the fractional derivatives of the SFPs.

Theorem 4. 

For $\eta \in (0,1)$, the fractional derivative of $F_{\mu +1}^{*}\left(t\right)$ can be represented as

$$D_t^{\eta }{F}_{\mu +1}^{*}\left(t\right)=t^{-\eta }\left(\sum _{\theta =0}^{\mu }{\gamma }_{\theta ,\mu }^{\eta }{F}_{\mu +1}^{*}\left(t\right)-o_{\mu }\right),$$

(49)

where

$${\gamma }_{\theta ,\mu }^{\eta }=\sum _{L=\theta }^{\mu }{\displaystyle \frac{L!B_{L,\mu }{H}_{\theta ,L}}{\mathsf{\Gamma }(L-\eta +1)}},$$

(50)

and

$$o_{\mu }={\displaystyle \frac{{(-1)}^{\mu }}{\mathsf{\Gamma }(1-\eta )}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1-\mu }{2}},-{\displaystyle \frac{\mu }{2}}\\ -\mu \end{array}|-4\right).$$

(51)

Proof. 

Formula (10) allows one to write $D_t^{\eta }{F}_{\mu +1}^{*}\left(t\right)$ as

$$D_t^{\eta }{F}_{\mu +1}^{*}\left(t\right)=\sum _{\theta =1}^{\mu }{\displaystyle \frac{\theta !B_{\theta ,\mu }}{\mathsf{\Gamma }(\theta -\eta +1)}}{t}^{\theta -\eta },$$

(52)

which can be written with the aid of Equation (21) as

$$D_t^{\eta }{F}_{\mu +1}^{*}\left(t\right)=t^{-\eta }\sum _{\theta =1}^{\mu }\sum _{L=0}^{\theta }{\displaystyle \frac{\theta !B_{\theta ,\mu }{H}_{L,\theta }}{\mathsf{\Gamma }(\theta -\eta +1)}}{F}_{L+1}^{*}\left(t\right)$$

(53)

that can be transformed into

$$D_t^{\eta }{F}_{\mu +1}^{*}\left(t\right)=t^{-\eta }\left(\sum _{\theta =0}^{\mu }{\gamma }_{\theta ,\mu }^{\eta }{F}_{\mu +1}^{*}\left(t\right)-o_{\mu }\right),$$

(54)

where

$${\gamma }_{\theta ,\mu }^{\eta }=\sum _{L=\theta }^{\mu }{\displaystyle \frac{L!B_{L,\mu }{H}_{\theta ,L}}{\mathsf{\Gamma }(L-\eta +1)}},$$

(55)

and

$$o_{\mu }={\displaystyle \frac{{(-1)}^{\mu }}{\mathsf{\Gamma }(1-\eta )}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1-\mu }{2}},-{\displaystyle \frac{\mu }{2}}\\ -\mu \end{array}|-4\right).$$

(56)

The proof is now complete.    □

The following corollary gives the operational matrix of the fractional derivative of the SFPs.

Corollary 4. 

The fractional derivative of ${\mathcal{F}}^{*}\left(t\right)$ can be written in matrix form as

$$D_t^{\eta }{\mathcal{F}}^{*}\left(t\right)=t^{-\eta }(\mathit{\gamma }{\mathcal{F}}^{*}\left(t\right)-\mathit{O}),$$

(57)

where ${\mathcal{F}}^{*}\left(t\right)={[F_1^{*}\left(t\right),F_2^{*}\left(t\right),\dots ,F_{\mathcal{N}+1}^{*}\left(t\right)]}^T$, $\mathit{O}={[o_0,o_1,\dots ,o_{\mathcal{N}}]}^T$. Additionally, $\mathit{\gamma }=\left({\gamma }_{\theta ,\mu }^{\eta }\right)$ is the operational matrix of the fractional derivative of order $(\mathcal{N}+1)\times (\mathcal{N}+1)$ whose entries are given in (51) and (50), respectively.

Proof. 

The matrix form in (57) directly results from the expression in (49).    □

4. Collocation Procedure for the TFBE

In this section, we are interested in handling the following TFBE [63]:

$${\displaystyle \frac{{\partial }^{\eta }u(\zeta ,t)}{\partial t^{\eta }}}+u(\zeta ,t){\displaystyle \frac{\partial u(\zeta ,t)}{\partial \zeta }}-\mathsf{\Psi }{\displaystyle \frac{{\partial }^{2}u(\zeta ,t)}{\partial {\zeta }^2}}=S(\zeta ,t),0<\eta <1,$$

(58)

governed by the following constraints:

$$\begin{array}{cc}\hfill u(\zeta ,0)=& u_1\left(\zeta \right),0<\zeta \le 1,\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill u(0,t)=& u_2\left(t\right),u(1,t)=u_3\left(t\right),0<t\le 1,\hfill \end{array}$$

(60)

where the source term is $S(\zeta ,t)$ and the kinematic viscosity is $\mathsf{\Psi }$.

The Proposed Collocation Algorithm

This section is confined to analyzing a numerical algorithm to solve (58) governed by the conditions (59) and (60).

Now, one may set

$${\mathcal{U}}^{\mathcal{N}}\left(\mathsf{\Omega }\right)=\mathrm{span}\{F_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right):0\le \nu ,\mu \le \mathcal{N}\},$$

where $\mathsf{\Omega }=(0,1)\times (0,1).$

Consequently, any function $u^{\mathcal{N}}(\zeta ,t)\in {\mathcal{U}}^{\mathcal{N}}\left(\mathsf{\Omega }\right)$ can be represented as

$$u^{\mathcal{N}}(\zeta ,t)=\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =0}^{\mathcal{N}}{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)=\mathcal{F}{\left(\zeta \right)}^T\mathit{C}\mathcal{F}\left(t\right),$$

(61)

where, $\mathcal{F}\left(\zeta \right)={[F_1^{*}\left(\zeta \right),F_2^{*}\left(\zeta \right),\dots ,F_{\mathcal{N}+1}^{*}\left(\zeta \right)]}^T,\mathcal{F}\left(t\right)={[F_1^{*}\left(t\right),F_2^{*}\left(t\right),\dots ,F_{\mathcal{N}+1}^{*}\left(t\right)]}^T,$ and $\mathit{C}={\left(c_{\nu \mu }\right)}_{0\le \nu ,\mu \le \mathcal{N}}$ is the unknown matrix whose order is ${(\mathcal{N}+1)}^2$.

Now, the residual $\mathbf{R}(\zeta ,t)$ of Equation (58) may be expressed as

$$\begin{array}{cc}\hfill \mathbf{R}(\zeta ,t)& ={\displaystyle \frac{{\partial }^{\eta }{u}^{\mathcal{N}}(\zeta ,t)}{\partial t^{\eta }}}+u^{\mathcal{N}}(\zeta ,t){\displaystyle \frac{\partial u^{\mathcal{N}}(\zeta ,t)}{\partial \zeta }}-\mathsf{\Psi }{\displaystyle \frac{{\partial }^2{u}^{\mathcal{N}}(\zeta ,t)}{\partial {\zeta }^2}}-S(\zeta ,t)\hfill \\ & =\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =0}^{\mathcal{N}}{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right){\displaystyle \frac{{\partial }^{\eta }{F}_{\mu +1}^{*}}{\partial t^{\eta }}}\left(t\right)+\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =0}^{\mathcal{N}}\sum _{m=0}^{\mathcal{N}}\sum _{n=0}^{\mathcal{N}}{c}_{\nu \mu }{c}_{mn}{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right){\displaystyle \frac{\partial F_{m+1}^{*}\left(\zeta \right)}{\partial \zeta }}{F}_{n+1}^{*}\left(t\right)\hfill \\ & -\mathsf{\Psi }\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =0}^{\mathcal{N}}{c}_{\nu \mu }{\displaystyle \frac{{\partial }^2{F}_{\nu +1}^{*}\left(\zeta \right)}{\partial {\zeta }^2}}{F}_{\mu +1}^{*}\left(t\right)-S(\zeta ,t).\hfill \end{array}$$

(62)

In virtue of Corollaries 3 and 4, along with the expansion (61), we may write $\mathbf{R}(\zeta ,t)$ (62) in matrix form as

$$\begin{array}{cc}\hfill \mathbf{R}(\zeta ,t)=& {\left({\mathcal{F}}^{*}\left(\zeta \right)\right)}^T\mathit{C}{t}^{-\eta }[\mathit{\gamma }{\mathcal{F}}^{*}\left(t\right)-\mathit{O}]+\left[{\left({\mathcal{F}}^{*}\left(\zeta \right)\right)}^T\mathit{C}{\mathcal{F}}^{*}\left(t\right)\right]\left[\mathit{\lambda }{\left({\mathcal{F}}^{*}\left(\zeta \right)\right)}^T\mathit{C}{\mathcal{F}}^{*}\left(t\right)\right]\hfill \\ & -\mathsf{\Psi }\left[\mathit{\beta }{\left({\mathcal{F}}^{*}\left(\zeta \right)\right)}^T\right]\mathit{C}{\mathcal{F}}^{*}\left(t\right)-S(\zeta ,t).\hfill \end{array}$$

(63)

To obtain $c_{\nu \mu }$, the collocation method is used in the following way: at specific nodes $({\zeta }_{\nu },t_{\mu })$, the residual $\mathbf{R}(\zeta ,t)$ is forced to be zero; that is

$$\mathbf{R}\left({\displaystyle \frac{\nu +1}{\mathcal{N}+2}},{\displaystyle \frac{\mu +1}{\mathcal{N}+2}}\right)=0,0\le \nu \le \mathcal{N}-2,0\le \mu \le \mathcal{N}-1.$$

(64)

Moreover, the constraints in (59) and (60) lead to

$$\begin{array}{cc}& \mathcal{F}{\left({\displaystyle \frac{\nu +1}{\mathcal{N}+2}}\right)}^T\mathcal{C}\mathcal{F}\left(0\right)=u_1\left({\displaystyle \frac{\nu +1}{\mathcal{N}+2}}\right),\nu =0,1,2,\dots ,\mathcal{N},\hfill \\ & \mathcal{F}{\left(0\right)}^T\mathcal{C}\mathcal{F}\left({\displaystyle \frac{\mu +1}{\mathcal{N}+2}}\right)=u_2\left({\displaystyle \frac{\mu +1}{\mathcal{N}+2}}\right),\mu =0,1,2,\dots ,\mathcal{N}-1,\hfill \\ & \mathcal{F}{\left(1\right)}^T\mathcal{C}\mathcal{F}\left({\displaystyle \frac{\mu +1}{\mathcal{N}+2}}\right)=u_3\left({\displaystyle \frac{\mu +1}{\mathcal{N}+2}}\right),\mu =0,1,2,\dots ,\mathcal{N}-1.\hfill \end{array}$$

(65)

Now, Equations (64) and (65) form a nonlinear system of equations with dimension ${(\mathcal{N}+1)}^2$ in the unknown expansion coefficients $c_{\nu \mu }$, which can be solved by applying Newton’s iterative approach.

Remark 3. 

The following Algorithm 1 lists the required steps to obtain the desired approximate solution to the fractional Burgers equation.

Algorithm 1 Coding algorithm for our proposed technique

Input  $\eta ,\mathcal{N},u_1\left(\zeta \right),u_2\left(t\right),u_3\left(t\right)$, and $S(\zeta ,t)$. Step 1. Assume an approximate solution $u^{\mathcal{N}}(\zeta ,t)$ as in (61). Step 2. Use Corollaries 3 and 4 along with the expansion (61) to obtain the residual $\mathbf{R}(\zeta ,t)$ in (63). Step 3. Apply the collocation method to obtain the system in (64) and (65). Step 4. Use FindRoot command with initial guess $\{c_{\nu \mu }={10}^{-\nu -\mu }:\nu ,\mu =0,1,\cdots ,\mathcal{N}\}$                 to solve the system in (64) and (65) to obtain $c_{\nu \mu }$. Output  $u^{\mathcal{N}}(\zeta ,t).$

5. The Convergence and Error Analysis

In this section, we study the convergence of the shifted Fibonacci expansion. So, some necessary lemmas will be used in this study. In addition, three theorems will be stated and proved as follows:

• The first lemma gives an upper estimate for $F_{\nu +1}^{*}\left(\zeta \right)$.
• The second lemma expresses the infinitely differentiable function $f\left(\zeta \right)$ in terms of $F_{\nu +1}^{*}\left(\zeta \right)$.
• The first theorem gives an upper estimate for the unknown expansion coefficients $a_n$.
• The second theorem gives an upper estimate for the unknown double expansion coefficients $c_{\nu \mu }$.
• The third theorem gives an upper estimate for the truncation error.

Lemma 2. 

Let $\zeta \in [0,1]$. The estimate is valid as follows:

$$F_{\nu +1}^{*}\left(\zeta \right)\le {\varphi }^{\nu +1},\forall \nu \ge 0,$$

where $\varphi ={\displaystyle \frac{1+\sqrt{5}}{2}}$ is the golden ratio [14].

Proof. 

Noting that $F_{\nu +1}^{*}\left(\zeta \right)\le F_{\nu +1}^{*}\left(1\right)=F_{\nu +1}$, where $F_i$ represents the Fibonacci numbers, and since $F_{\nu +1}$ is the closest integer to ${\displaystyle \frac{{\varphi }^{\nu +1}}{\sqrt{5}}}$, we obtain the desired result. □

Lemma 3. 

Consider the infinitely differentiable function $f\left(\zeta \right)$ at $\zeta =0$. We have

$$f\left(\zeta \right)=\sum _{n=0}^{\infty }\sum _{s=n}^{\infty }{\displaystyle \frac{f^{\left(s\right)}\left(0\right)H_{n,s}}{s!}}{F}_{n+1}^{*}\left(\zeta \right).$$

(66)

Proof. 

Assume the following expansion of $f\left(\zeta \right)$:

$$f\left(\zeta \right)=\sum _{\theta =0}^{\infty }{\displaystyle \frac{f^{\left(\theta \right)}\left(0\right)}{\theta !}}{\zeta }^{\theta }.$$

(67)

As a result of the inversion Formula (21), the last equation becomes

$$f\left(\zeta \right)=\sum _{\theta =0}^{\infty }\sum _{L=0}^{\theta }{\displaystyle \frac{f^{\left(\theta \right)}\left(0\right)H_{L,\theta }}{\theta !}}{F}_{L+1}^{*}\left(\zeta \right),$$

(68)

that can be written alternatively as

$$f\left(\zeta \right)=\sum _{n=0}^{\infty }\sum _{s=n}^{\infty }{\displaystyle \frac{f^{\left(s\right)}\left(0\right)H_{n,s}}{s!}}{F}_{n+1}^{*}\left(\zeta \right).$$

(69)

This finalizes the proof. □

Theorem 5. 

If $f\left(\zeta \right)$ is defined on $[0,1]$ and $|f^{\left(\nu \right)}\left(0\right)|\le {\lambda }^{\nu },\nu >0,$ where $0<\lambda <1$ and ${\displaystyle f\left(\zeta \right)=\sum _{n=0}^{\infty }{a}_n{F}_{n+1}^{*}\left(\zeta \right)}$, we obtain

$$|a_n|\le {\left({\displaystyle \frac{\lambda }{2}}\right)}^n.$$

(70)

Moreover, the series converges absolutely.

Proof. 

In virtue of Lemma 3, one has

$$a_n=\sum _{s=n}^{\infty }{\displaystyle \frac{f^{\left(s\right)}\left(0\right)H_{n,s}}{s!}},$$

(71)

and, therefore

$$\begin{array}{cc}\hfill |a_n|& \le \sum _{s=n}^{\infty }{\displaystyle \frac{{\lambda }^s\left|H_{n,s}\right|}{s!}}\hfill \\ \hfill & \le \sum _{s=n}^{\infty }\left\{\begin{array}{cc}{\displaystyle \frac{(n+1)2^{-s}{\lambda }^s}{\left(\frac{s-n}{2}\right)!\mathsf{\Gamma }\left(\frac{1}{2}(n+s+4)\right)}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1}{2}}(-n-s-2),{\displaystyle \frac{n-s}{2}}\\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right),\hfill & (n+s)\mathrm{even},\hfill \\ {\displaystyle \frac{(n+1)2^{-s}}{\mathsf{\Gamma }\left(\frac{1}{2}(-n+s+1)\right)\mathsf{\Gamma }\left(\frac{1}{2}(n+s+3)\right)}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1}{2}}(-n-s-1),{\displaystyle \frac{1}{2}}(n-s+1)\\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right),\hfill & (n+s)\mathrm{odd}.\hfill \end{array}\right.\hfill \end{array}$$

(72)

Now, for all $0\le n<\infty$ and $n\le s<\infty$, the following estimations can be deduced

$$\begin{array}{cc}& {\displaystyle \frac{(n+1)2^{-s}{\lambda }^s}{\left(\frac{s-n}{2}\right)!\mathsf{\Gamma }\left(\frac{1}{2}(n+s+4)\right)}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1}{2}}(-n-s-2),{\displaystyle \frac{n-s}{2}}\\ {\displaystyle \frac{1}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\le {\left({\displaystyle \frac{\lambda }{2}}\right)}^n,\hfill \\ & {\displaystyle \frac{(n+1)2^{-s}}{\mathsf{\Gamma }\left(\frac{1}{2}(-n+s+1)\right)\mathsf{\Gamma }\left(\frac{1}{2}(n+s+3)\right)}}{}_{2}F_1\left(\begin{array}{c}{\displaystyle \frac{1}{2}}(-n-s-1),{\displaystyle \frac{1}{2}}(n-s+1)\\ {\displaystyle \frac{3}{2}}\end{array}|-{\displaystyle \frac{1}{4}}\right)\le {\left({\displaystyle \frac{\lambda }{2}}\right)}^n.\hfill \end{array}$$

(73)

Therefore, we obtain the estimation

$$|a_n|\le {\left({\displaystyle \frac{\lambda }{2}}\right)}^n.$$

(74)

To prove the second part of the theorem, since

$$\sum _{n=0}^{\infty }|a_n{F}_{n+1}^{*}\left(\zeta \right)|=\sum _{n=0}^{\infty }|a_n||F_{n+1}^{*}\left(\zeta \right)|\le {\displaystyle \frac{{\lambda }^n{\varphi }^{n+1}}{2^n}}={\displaystyle \frac{2\varphi }{2-\lambda \varphi }},$$

(75)

so the series converges absolutely. □

Theorem 6. 

If a function ${\displaystyle u(\zeta ,t)=u_1\left(\zeta \right)u_2\left(t\right)=\sum _{\nu =0}^{\infty }\sum _{\mu =0}^{\infty }{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)}$, with $|u_1^{\left(\nu \right)}\left(0\right)|\le {\lambda }_1^{\nu }$ and $|u_2^{\left(\nu \right)}\left(0\right)|\le {\lambda }_2^{\nu }$, where $0<{\lambda }_1,{\lambda }_2<1$, then one has

$$|c_{\nu \mu }|\le {\displaystyle \frac{{{\lambda }_1}^{\nu }{{\lambda }_2}^{\mu }}{2^{\nu +\mu }}}.$$

(76)

Moreover, the series converges absolutely.

Proof. 

Using Lemma 3 along with the assumptions of the theorem $u(\zeta ,t)=u_1\left(\zeta \right)u_2\left(t\right)$, we can write

$$u(\zeta ,t)=\sum _{\nu =0}^{\infty }\sum _{\mu =0}^{\infty }\sum _{s=\nu }^{\infty }\sum _{r=\mu }^{\infty }{\displaystyle \frac{u_1^{\left(s\right)}\left(0\right)u_2^{\left(r\right)}\left(0\right)H_{\mu ,r}{H}_{\nu ,s}}{s!r!}}{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right).$$

(77)

Now, the expansion coefficients $c_{\nu \mu }$ can be written as

$$c_{\nu \mu }=\sum _{s=\nu }^{\infty }\sum _{r=\mu }^{\infty }{\displaystyle \frac{u_1^{\left(s\right)}\left(0\right)u_2^{\left(r\right)}\left(0\right)H_{\mu ,r}{H}_{\nu ,s}}{s!r!}}.$$

(78)

Using the assumption $|u_1^{\left(\nu \right)}\left(0\right)|\le {\lambda }_1^{\nu }$ and $|u_2^{\left(\nu \right)}\left(0\right)|\le {\lambda }_2^{\nu }$, one obtains

$$|c_{\nu \mu }|\le \sum _{s=\nu }^{\infty }{\displaystyle \frac{{\lambda }_1^s\left|H_{\nu ,s}\right|}{s!}}\times \sum _{r=\mu }^{\infty }{\displaystyle \frac{{\lambda }_2^r\left|H_{\mu ,r}\right|}{r!}}.$$

(79)

We obtain the desired result by performing similar steps as in the proof of Theorem 5. □

Theorem 7. 

Let $u(\zeta ,t)$ meet the assumptions of Theorem 6. This upper estimate on the truncation error holds as follows:

$$|u(\zeta ,t)-u^{\mathcal{N}}(\zeta ,t)|<{\displaystyle \frac{2{\varphi }^{\mathcal{N}+3}({\lambda }_1^{\mathcal{N}+1}+{\lambda }_2^{\mathcal{N}+1})}{2^{\mathcal{N}}(2-{\lambda }_1\varphi )(2-{\lambda }_2\varphi )}}.$$

Proof. 

The definitions of $u(\zeta ,t)$ and $u^{\mathcal{N}}(\zeta ,t)$ imply that

$$\begin{array}{cc}\hfill \left|u(\zeta ,t)-u^{\mathcal{N}}(\zeta ,t)\right|& =\left|\sum _{\nu =0}^{\infty }\sum _{\mu =0}^{\infty }{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)-\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =0}^{\mathcal{N}}{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)\right|\hfill \\ & \le \left|\sum _{\nu =0}^{\mathcal{N}}\sum _{\mu =\mathcal{N}+1}^{\infty }{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)\right|+\left|\sum _{\nu =\mathcal{N}+1}^{\infty }\sum _{\mu =0}^{\infty }{c}_{\nu \mu }{F}_{\nu +1}^{*}\left(\zeta \right)F_{\mu +1}^{*}\left(t\right)\right|.\hfill \end{array}$$

(80)

If we apply Theorem 6, Lemma 2, and the following inequalities

$$\begin{array}{cc}& \sum _{\nu =0}^{\mathcal{N}}{\displaystyle \frac{{\lambda }_1^{\nu }{\varphi }^{\nu +1}}{2^{\nu }}}={\displaystyle \frac{2\varphi -2^{-\mathcal{N}}{\lambda }_1^{\mathcal{N}+1}{\varphi }^{\mathcal{N}+2}}{2-{\lambda }_1\varphi }}<{\displaystyle \frac{2\varphi }{2-{\lambda }_1\varphi }},\hfill \\ & \sum _{\nu =\mathcal{N}+1}^{\infty }{\displaystyle \frac{{\lambda }_1^{\nu }{\varphi }^{\nu +1}}{2^{\nu }}}={\displaystyle \frac{{\left({\lambda }_1\varphi \right)}^{\mathcal{N}+1}\varphi }{2^{\mathcal{N}}(2-{\lambda }_1\varphi )}},\hfill \\ & \sum _{\nu =0}^{\infty }{\displaystyle \frac{{\lambda }_1^{\nu }{\varphi }^{\nu +1}}{2^{\nu }}}={\displaystyle \frac{2\varphi }{2-{\lambda }_1\varphi }},\hfill \end{array}$$

(81)

we obtain the following desired estimation:

$$\begin{array}{cc}\hfill \left|u(\zeta ,t)-u^{\mathcal{N}}(\zeta ,t)\right|& <{\displaystyle \frac{{\left({\lambda }_2\varphi \right)}^{\mathcal{N}+1}\varphi }{2^{\mathcal{N}}(2-{\lambda }_2\varphi )}}{\displaystyle \frac{2\varphi }{(2-{\lambda }_1\varphi )}}+{\displaystyle \frac{{\left({\lambda }_1\varphi \right)}^{\mathcal{N}+1}\varphi }{2^{\mathcal{N}}(2-{\lambda }_1\varphi )}}{\displaystyle \frac{2\varphi }{(2-{\lambda }_2\varphi )}}\hfill \\ & ={\displaystyle \frac{2{\varphi }^{\mathcal{N}+3}({\lambda }_1^{\mathcal{N}+1}+{\lambda }_2^{\mathcal{N}+1})}{2^{\mathcal{N}}(2-{\lambda }_1\varphi )(2-{\lambda }_2\varphi )}}.\hfill \end{array}$$

(82)

This completes the proof of this theorem. □

6. Illustrative Examples

This section is confined to presenting some test problems to show the performance and applicability of our proposed algorithm. All codes were written and debugged by Mathematica 11 on HP Z420 Workstation, Processor: Intel (R) Xeon(R) CPU E5-1620 v2—3.70 GHz, 16 GB Ram DDR3, and 512 GB storage.

Example 1 

([63,64]). Consider the TFBE of the form

$${\displaystyle \frac{{\partial }^{\eta }u(\zeta ,t)}{\partial t^{\eta }}}+u(\zeta ,t){\displaystyle \frac{\partial u(\zeta ,t)}{\partial \zeta }}-2{\displaystyle \frac{{\partial }^{2}u(\zeta ,t)}{\partial {\zeta }^2}}=S(\zeta ,t),$$

(83)

controlled by

$$\begin{array}{cc}& u(\zeta ,0)=0,0<\zeta \le 1,\hfill \\ & u(0,t)=u(1,t)=0,0<t\le 1,\hfill \end{array}$$

(84)

where $S(\zeta ,t)={\displaystyle \frac{2}{\mathsf{\Gamma }(3-\eta )}}{t}^{2-\eta }sin\left(\pi \zeta \right)+\pi t^{4}sin\left(\pi \zeta \right)cos\left(\pi \zeta \right)+2{\pi }^2{t}^{2}sin\left(\pi \zeta \right)$. In addition, $u(\zeta ,t)=t^{2}sin\left(\pi \zeta \right)$ is the exact solution to (83).

At various η values,

Table 1. Errors computed in $L^{\infty }$—sense of Example 1.

$\mathit{\eta }$	Method in [64] ($\mathit{N}={\mathbf{2}}^{\mathbf{7}}$, $\mathit{M}={\mathbf{2}}^{\mathbf{12}}$)	Method in [63] ($\mathit{N}={\mathbf{2}}^{\mathbf{7}}$, $\mathit{M}={\mathbf{2}}^{\mathbf{12}}$)	Our Method ($\mathcal{N}=\mathbf{12}$)	Our CPU Time
0.2	$5.47·{10}^{-5}$	$4.77·{10}^{-5}$	$6.09·{10}^{-9}$	70.484
0.3	$5.46·{10}^{-5}$	$4.75·{10}^{-5}$	$4.77·{10}^{-9}$	74.671
0.4	$5.45·{10}^{-5}$	$4.74·{10}^{-5}$	$3.89·{10}^{-9}$	74.78
0.5	$5.43·{10}^{-5}$	$4.72·{10}^{-5}$	$3.88·{10}^{-9}$	72.704

compares our technique to the techniques in [63,64] with respect to the $L^{\infty }$-error. In addition, the CPU time (in seconds) for our proposed method is presented in this table. Our technique and the method in [63] are compared at $\mathcal{N}=12$ and various t values when $\eta =0.2$ and $\eta =0.5$, respectively, in

Table 2. Errors computed in $L^{\infty }$—sense of Example 1.

	Method in [63]			Our Method at $\mathcal{N}=12$
$\mathit{t}$	$\mathit{\eta }=\mathbf{0}.\mathbf{2},\mathit{N}=\mathbf{500}$	$\mathit{\eta }=\mathbf{0}.\mathbf{5},\mathit{N}=\mathbf{250}$		$\mathit{\eta }=\mathbf{0}.\mathbf{2}$	$\mathit{\eta }=\mathbf{0}.\mathbf{5}$
0.2	$3.09·{10}^{-7}$	$3.24·{10}^{-7}$		$7.88·{10}^{-10}$	$1.43·{10}^{-10}$
0.4	$1.49·{10}^{-6}$	$2.20·{10}^{-6}$		$9.93·{10}^{-10}$	$9.36·{10}^{-10}$
0.6	$3.47·{10}^{-6}$	$5.60·{10}^{-6}$		$1.63·{10}^{-9}$	$1.99·{10}^{-9}$
0.8	$6.22·{10}^{-6}$	$1.04·{10}^{-6}$		$2.12·{10}^{-9}$	$3.27·{10}^{-9}$
1.0	$9.71·{10}^{-6}$	$1.67·{10}^{-6}$		$6.09·{10}^{-9}$	$4.28·{10}^{-9}$

. At $\eta =0.7$ and $\eta =0.8$,

Table 3. Comparison of AE for Example 1 at $t=1$.

	$\mathit{\eta }=0.7$			$\mathit{\eta }=0.8$
$\mathit{\zeta }$	Method in [63]	Our Method at$\mathcal{N}=\mathbf{12}$		Method in [63]	Our Method at$\mathcal{N}=\mathbf{12}$
0.1	$2.19·{10}^{-5}$	$4.06·{10}^{-11}$		$8.11·{10}^{-7}$	$5.46·{10}^{-10}$
0.2	$4.21·{10}^{-5}$	$6.36·{10}^{-10}$		$1.57·{10}^{-6}$	$1.09·{10}^{-9}$
0.3	$5.89·{10}^{-5}$	$1.24·{10}^{-9}$		$2.23·{10}^{-6}$	$1.66·{10}^{-9}$
0.4	$7.07·{10}^{-5}$	$1.88·{10}^{-9}$		$2.72·{10}^{-6}$	$2.30·{10}^{-9}$
0.5	$7.61·{10}^{-5}$	$2.57·{10}^{-9}$		$2.97·{10}^{-6}$	$3.03·{10}^{-9}$
0.6	$7.42·{10}^{-5}$	$3.30·{10}^{-9}$		$2.94·{10}^{-6}$	$3.86·{10}^{-9}$
0.7	$6.45·{10}^{-5}$	$4.07·{10}^{-9}$		$2.59·{10}^{-6}$	$4.80·{10}^{-9}$
0.8	$4.77·{10}^{-5}$	$4.86·{10}^{-9}$		$1.94·{10}^{-6}$	$5.85·{10}^{-9}$
0.9	$2.54·{10}^{-5}$	$5.43·{10}^{-9}$		$1.04·{10}^{-6}$	$6.75·{10}^{-9}$

compares our method to the method in [63] in terms of absolute error (AE). When $\mathcal{N}=12$, the AE is shown in Figure 1 for various η values.

Example 2 

([63]). Consider the TFBE of the form

$${\displaystyle \frac{{\partial }^{\eta }u(\zeta ,t)}{\partial t^{\eta }}}+u(\zeta ,t){\displaystyle \frac{\partial u(\zeta ,t)}{\partial \zeta }}-{\displaystyle \frac{{\partial }^{2}u(\zeta ,t)}{\partial {\zeta }^2}}=S(\zeta ,t),$$

(85)

controlled by

$$\begin{array}{cc}& u(\zeta ,0)=0,0<\zeta \le 1,\hfill \\ & u(0,t)=t^2,u(1,t)=-t^2,0<t\le 1,\hfill \end{array}$$

(86)

where $S(\zeta ,t)={\displaystyle \frac{2}{\mathsf{\Gamma }(3-\eta )}}{t}^{2-\eta }cos\left(\pi \zeta \right)-\pi t^{4}sin\left(\pi \zeta \right)cos\left(\pi \zeta \right)+{\pi }^2{t}^{2}cos\left(\pi \zeta \right)$. In addition, $u(\zeta ,t)=t^{2}cos\left(\pi \zeta \right)$ is the exact solution to (85).

At various η values,

Table 4. Errors computed in $L^{\infty }$—sense of Example 2.

	Method in [63]
$\mathit{\eta }$	$\mathit{N}={\mathbf{2}}^{\mathbf{7}}$,$\mathit{M}={\mathbf{2}}^{\mathbf{12}}$	$\mathit{M}={\mathbf{2}}^{\mathbf{7}}$,$\mathit{N}={\mathbf{2}}^{\mathbf{11}}$	Our Method at $\mathcal{N}=\mathbf{11}$	Our CPU Time
0.2	$1.06·{10}^{-5}$	$9.67·{10}^{-6}$	$4.70·{10}^{-8}$	54.813
0.3	$1.06·{10}^{-5}$	$9.61·{10}^{-6}$	$4.67·{10}^{-8}$	54.703
0.4	$1.06·{10}^{-5}$	$9.53·{10}^{-6}$	$4.65·{10}^{-8}$	56.171
0.5	$1.06·{10}^{-5}$	$9.42·{10}^{-6}$	$4.63·{10}^{-8}$	51.283

compares our method to the method in [63] in the sense of $L^{\infty }$-error. Moreover, the CPU time (in seconds) for our proposed method is presented in this table. The $L^{\infty }$-error at $\mathcal{N}=11$ compared to that obtained using the technique described in [63] at $\eta =0.9$ is shown in

Table 5. Errors computed in $L^{\infty }$—sense of Example 2 at $\eta =0.9$.

	Method in [63]
$\mathit{t}$	$\mathit{N}=\mathbf{500}$, $\Delta \mathit{t}=\mathbf{0.001}$	$\mathit{N}=\mathbf{210}$, $\Delta \mathit{t}=\mathbf{0.0025}$	$\mathit{N}=\mathbf{130}$, $\Delta \mathit{t}=\mathbf{0.004}$	Our Method at $\mathcal{N}=\mathbf{11}$
0.2	$2.11·{10}^{-9}$	$5.57·{10}^{-9}$	$1.01·{10}^{-8}$	$1.83·{10}^{-9}$
0.4	$9.22·{10}^{-8}$	$5.16·{10}^{-7}$	$1.35·{10}^{-6}$	$8.03·{10}^{-8}$
0.6	$2.54·{10}^{-7}$	$1.45·{10}^{-6}$	$3.78·{10}^{-6}$	$1.90·{10}^{-9}$
0.8	$4.89·{10}^{-7}$	$2.81·{10}^{-6}$	$7.30·{10}^{-6}$	$3.51·{10}^{-8}$
1.0	$7.99·{10}^{-7}$	$4.60·{10}^{-6}$	$1.20·{10}^{-5}$	$5.71·{10}^{-8}$

. When $\eta =0.7$, the AE is reported in

Table 6. The AE of Example 2 at different values of $\mathcal{N}$ at $\eta =0.7$.

$(\mathit{\zeta },\mathit{t})$	$\mathcal{N}=7$	$\mathcal{N}=9$	$\mathcal{N}=11$
(0.1, 0.1)	$1.33·{10}^{-7}$	$2.15·{10}^{-9}$	$2.41·{10}^{-11}$
(0.3, 0.3)	$4.94·{10}^{-6}$	$8.95·{10}^{-8}$	$1.07·{10}^{-9}$
(0.5, 0.5)	$2.91·{10}^{-5}$	$5.35·{10}^{-7}$	$6.44·{10}^{-9}$
(0.7, 0.7)	$9.48·{10}^{-5}$	$1.75·{10}^{-6}$	$2.11·{10}^{-8}$
(0.9, 0.9)	$1.65·{10}^{-4}$	$3.53·{10}^{-6}$	$4.58·{10}^{-8}$

for various $\mathcal{N}$ values. The AE at various $\mathcal{N}$ values when $\eta =0.9$ is illustrated in Figure 2.

Example 3 

([63]). Consider the TFBE of the form

$${\displaystyle \frac{{\partial }^{\eta }u(\zeta ,t)}{\partial t^{\eta }}}+u(\zeta ,t){\displaystyle \frac{\partial u(\zeta ,t)}{\partial \zeta }}-{\displaystyle \frac{{\partial }^{2}u(\zeta ,t)}{\partial {\zeta }^2}}=S(\zeta ,t),$$

(87)

controlled by

$$\begin{array}{cc}& u(\zeta ,0)=0,0<\zeta \le 1,\hfill \\ & u(0,t)=t^2,u(1,t)=e{t}^2,0<t\le 1,\hfill \end{array}$$

(88)

where $S(\zeta ,t)={\displaystyle \frac{2}{\mathsf{\Gamma }(3-\eta )}}{e}^{\zeta }{t}^{2-\eta }+t^4{e}^{2\zeta }-t^2{e}^{\zeta }$ and $u(\zeta ,t)=t^2{e}^{\zeta }$ is the exact solution to this problem.

Table 7. Comparison of AE for Example 3 at $\eta =0.9$ at $t=1$.

	Method in [63]
$\mathit{\zeta }$	$\mathit{N}=\mathbf{200},\Delta \mathit{t}=\mathbf{0}.\mathbf{0005}$	$\mathit{N}=\mathbf{80},\Delta \mathit{t}=\mathbf{0}.\mathbf{001}$	Our Method at  $\mathcal{N}=\mathbf{11}$
0.1	$6.18·{10}^{-9}$	$2.24·{10}^{-7}$	$1.45·{10}^{-12}$
0.2	$7.15·{10}^{-9}$	$4.00·{10}^{-7}$	$8.41·{10}^{-13}$
0.3	$3.50·{10}^{-9}$	$5.30·{10}^{-7}$	$2.77·{10}^{-12}$
0.4	$4.00·{10}^{-9}$	$6.15·{10}^{-7}$	$1.17·{10}^{-11}$
0.5	$1.42·{10}^{-8}$	$6.52·{10}^{-7}$	$3.01·{10}^{-11}$
0.6	$2.56·{10}^{-8}$	$6.41·{10}^{-7}$	$6.35·{10}^{-11}$
0.7	$3.55·{10}^{-8}$	$5.77·{10}^{-7}$	$1.21·{10}^{-10}$
0.8	$3.98·{10}^{-8}$	$4.55·{10}^{-7}$	$2.18·{10}^{-10}$
0.9	$3.18·{10}^{-8}$	$2.67·{10}^{-7}$	$3.45·{10}^{-10}$

presents a comparison of the AE between our method and the technique described in [63] for various values of ζ and $\eta =0.9$ when $t=1$.

Table 8. Errors computed in $L^{\infty }$—sense of Example 3.

	Method in [63]
$\mathit{\eta }$	$\mathit{N}={\mathbf{2}}^{\mathbf{7}}$,  $\mathit{M}={\mathbf{2}}^{\mathbf{11}}$	$\mathit{M}={\mathbf{2}}^{\mathbf{7}}$,  $\mathit{N}={\mathbf{2}}^{\mathbf{4}}$	Our Method at $\mathcal{N}=\mathbf{11}$	Our CPU Time
0.2	$5.69·{10}^{-7}$	$1.97·{10}^{-5}$	$7.23·{10}^{-12}$	53.639
0.3	$5.66·{10}^{-7}$	$2.00·{10}^{-5}$	$1.45·{10}^{-12}$	54.438
0.4	$5.61·{10}^{-7}$	$2.04·{10}^{-5}$	$1.21·{10}^{-12}$	54.376
0.5	$5.55·{10}^{-7}$	$2.09·{10}^{-5}$	$1.67·{10}^{-12}$	52.47

presents a comparison of the $L^{\infty }$-error between our method with $\mathcal{N}=11$ and the technique described in [63] at various values of η. In addition, the CPU time (in seconds) for our proposed method is presented in this table. Figure 3 displays the AE at various values of η when $\mathcal{N}=11$.

Remark 4. 

In light of Figure 1, Figure 2 and Figure 3, we can conclude the following benefit:

Excellently precise approximations can be obtained by selecting shifted Fibonacci polynomials $F_{\nu +1}^{*}\left(\zeta \right)$ as basis functions and taking specific terms of the retained modes.

7. Concluding Remarks

In this work, we presented new numerical solutions for the fractional Burgers equation. New basis functions of the shifted Fibonacci polynomials were introduced and utilized to obtain the desired numerical solutions. We established new integer and fractional derivatives formulas for the shifted Fibonacci polynomials to be capable of designing the matrix collocation algorithm. The accuracy of the shifted Fibonacci expansion is tested theoretically and numerically by presenting some examples. As far as we know, this is the first time this kind of polynomial is utilized in the scope of numerical solutions of DEs. Additionally, we think this study has unveiled a new horizon in the numerical analysis of FDEs. In the near future, we hope to utilize these polynomials to treat some other important differential equations of different types.