Chebyshev fifth-kind series approximation for generalized space fractional partial differential equations

1.   Introduction

Many phenomena such as biology, physics, and fluid mechanics can be modeled by certain fractional order partial differential equations (FPDEs) ^[1,2,3]. So that, fractional calculus becomes a central branch of mathematical analysis. The importance of the numerical solution of FPDEs becomes a major, because of the difficulty of obtaining their analytical solutions ^[4,5].

Spectral methods have been developed through the past few decades by a huge number of researchers see for example ^[6,7,8], and many others. The principal feature of these methods lies in their ability to reach acceptably accurate results with substantially fewer degrees of freedom. In recent years, Chebyshev polynomials have become increasingly important again in numerical analysis, when a new two classes of polynomials appear, namely fifth and sixth kinds ^{[9,10,11,12,13,14]}. In the Ph.D. thesis of Masjed-Jamei ^[15,16], 2006 he introduces a generalized polynomial using an extended Sturm-Liouville problem. These generalized polynomials generate Chebyshev polynomials of the first, second, third, and fourth kinds, in addition to the two new classes fifth and sixth kinds obtained at special values of a given parameters.

The objective of this research paper is to present a spectral scheme according to the collocation method for the generalized space-fractional partial differential equations (GFPDEs) that we have introduced. The proposed GFPDEs are chosen to be linear and the fractional derivatives are expressed in terms of Caputo's definition. The method of solution is to apply Shifted fifth kind Chebyshev polynomials using the collocation method to discretize the proposed equation, and then generate a linear system of ordinary differential equations (SODEs), which reduces the proposed problem. Additionally, to treat the generated SODEs, the classical fourth-order Runge-Kutta method (RK4) and the finite difference method (FDM) as well, are used. The proposed equation is presented as:

on a finite domain $0 < x \leqslant L; \; \; 0 < t \leqslant T$ and the parameters $\gamma_k$ refers to the fractional orders of a special derivative with $k < \gamma_k < (k+1) \leqslant n$. The function $f(x, t)$ is the source term, the functions $Q_{k}(x)$ are well defined and known, in addition $P$ is real constant. We also assume the initial condition (IC) as:

and the boundary conditions (BCs):

The introduced GFPDEs (1.1) represent a great generalization of a significant types of many physical models. As a special cases: At $\gamma_1\neq0, \; \; \gamma_k = 0,$ Eq (1.1) reduces to space-fractional order diffusion equation, and when $\gamma_0, \; \gamma_1\neq0, \; \; \gamma_k = 0,$ then, (1.1) becomes space-fractional order advection-dispersion equation, which they well study in the application section.

Concerning the existence and uniqueness of the solution of Eq (1.1), we refer the reader to references ^[17,18,19], these studies considered the existence and uniqueness of the solution for generalized linear and non-linear models of FPDEs, and we note that (1.1) represents a special case from ^[19]. In particular, we mention the references which prove the existence and uniqueness of the main examples mentioned in this work, the fractional order diffusion equation see ^[20,21], and advection-dispersion equation in ^[22].

2.   General notations

In this section, some definitions and properties for the fractional derivative and fifth kind Chebyshev polynomials are listed ^{[10,15,16,23]}.

2.1.   The Caputo fractional derivative

The Caputo's fractional derivative operator $D^{\gamma}_t$ (insted of $^{C}_{0}D^{\gamma}_t$ for short) of order $\gamma$ is characterize in the following form:

where $\ x > 0$, $n-1 < \gamma\leq n, \ n \in \mathbb{N}_0,$ and $\mathbb {N}_0 = \mathbb {N}-\{0\}$.

${D^{\gamma}_t}\, \sum_{i = 0}^{m}\lambda_i \, \varPsi_i(x) = \sum_{i = 0}^{m}\lambda_i\, {D^{\gamma}_t}\varPsi_i(x)$, where $\lambda_i$ and $\gamma$ are constants.

The Caputo fractional differentiation of a constant is zero.

${D^{\gamma}_t}\, x^{k} = \left\{ \begin{array}{ll} 0, & \text{ for $k\in \mathbb{N}_{0}$ and $k < \lceil \gamma \rceil$ } \\ \frac{\Gamma(k+1)\, x^{k-\gamma}}{\Gamma(k+1-\gamma)}, & \text{ for $k\in \mathbb{N}_{0}$ and $k\geq\lceil \gamma \rceil$ } \\ \end{array} \right.$, where $\lceil\gamma\rceil$ denote to the smallest integer greater than or equal to $\gamma$.

Remark 1. In this work we write the fractional Caputo's operator symbol $D^{ \gamma}$ instead of $^{C}_{0}D^{ \gamma}_{x}$ for short.

2.2.   Fifth-kind Chebyshev polynomials

The Chebyshev polynomials of the fifth-kind $\bar{X}_{n}(x)$ defined as: an orthonormonal polynomials in $x$ of degree $n$ defined on the closed interval $[-1, 1],$ The polynomials $\bar{X}_{n}(x)$ are orthogonal and the orthogonality relation is:

In ^[16] the authors normalize the monic Chebyshev polynomials of the fifth kind, and define $X_{n}(x)$ as:

and

and (2.2) may rewrite using $X_n(x)$ as:

By the usual transformation the Shifted Chebyshev polynomials of the fifth-kind $C_{n}(x)$ defined as:

The Shifted Chebyshev fifth-kind $C_{n}(x)$ are orthogonal polynomials on the closed interval $[0, 1]$, and they can be generated by using the following recurrence relation

with $\hbar_i$ defined in (2.4) and $\beta_{i+1}$ is

from (2.5), it is not difficult to note that $C_n (x), n \geq 0,$ are orthonormal on $[0, 1]$, and they have an orthogonality relation as:

Proposition 1. The Shifted polynomials $C_n(x)$ are defined through the Shifted first kind $T^*_n(x)$ by the following formula

where

and

The proof of Proposition 1 is given in ^[16]. Nevertheless, the Shifted Chebyshev polynomials of the first kind $T^{*}_n(x)$ are defined on $[0, 1]$, and they can be generated by using the following recurrence relation:

where

Therefore, the analytic form of $T^{*}_n(x)$ of degree $n$ is given by:

According to Proposition 1 the following Corollary is easy to prove.

Corollary 2.1. Shifted Chebyshev polynomials of the fifth-kind $C_{n}(x)$ be explicitly expressed in terms of $T^*_{n} (x)$ in the following form:

and

where $\delta_k$ is defined before in Proposition 1.

Corollary 2.2. Shifted Chebyshev polynomials of the fifth-kind $C_{n}(x)$ be explicitly expressed in terms of $x^{n}$, or the analytic form in the following form:

where

where $\delta_k$ is defined before, and $\lfloor.\rfloor$ is the floor function.

According to relations (2.9), (2.11) and (2.12) the first four terms of $C_n(x)$ are:

where $C_0(x)$ and $C_1(x)$ are used as initials the recurance relation (2.7).

Proposition 1 gives the connection formulae of the fifth-kind Chebyshev polynomials and the Shifted first kind Chebyshev polynomials, therefore, the fifth-kind inherits from them its ability, boundness and convergence.

Lemma 2.1. The Shifted fifth-kind $Y^*_{n}(x)$ are bounded according to the following form:

The full proof is in ^[16,24], and it may directly given from the connection relation $(2.9)$.

3.   Procedures the approximate solution

In the spectral method, in contrast, the function $\phi(x)$ may be expanded by Chebyshev polynomials of the fifth-kind series, which $\phi(x)$ is a square-integrable in $[0, 1]$, ^[12,25]:

Lemma 3.1. The infinite series $(3.1)$ is convergent uniformly to $\phi (x)$, and the following relation holds:

therefore, $L$ is some positive constant provided from:

Lemma 3.2. The global error $e_N(x)$ for the function $\phi(x)$ defined in $(3.1)$, such that: $e_N(x) = \sum_{n = N+1}^{\infty}a_{n}C{n}(x)$, is bounded and the following relation is valid:

The proofs of Lemma 3.1 and Lemma 3.2 are found in ^[16,24], and it refers to that the error almost tends to zero in the case of a large $N$. Subsequently, by truncate series (3.1) to $N < \infty$, then the approximate $\phi(x)$ by a finite sum of $(n+1)-$terms expressed in the following form:

The coefficients $a_{n}$ in relation (3.5) are given by the following relation:

Lemma 3.3. The local error $\tilde{e}_N(x)$ for the function $\phi(x)$ defined in $(3.1)$, such that: $\tilde{e}_N(x) = |\phi_{N+1}(x)-\phi_N(x)|$, is of order $N^3$.

The proof of Lemma 3.3 is completed in ^[16], and Gangi et al. ^[24] are increase in proof the form of the supremum for the local error.

Theorem 1. The fractional derivative of order $\gamma$ for the polynomials $C_n (x)$ according to Caputo's operator is given by:

and

where, $\varrho_{k, n}$ defined in Corollary 2.2.

Proof. According to (2.1) (the Caputo's operator) and the relation given in Corollary 2.2 it is easy to obtain the result, for more details see ^[16,24].

Theorem 2. Assume that, $\phi_N(x)$ be approximated function of $\phi(x)$ in terms of Shifted Chebyshev polynomials of the fifth kind as $(3.5)$, then the Caputo fractional derivative of order $\gamma$ when operating $\phi_N(x)$ is given by:

where, $\varrho_{k, n}$ defined in Corollary 2.2.

Proof. According to Theorem.1 and relation (3.5) one optians:

then the result (3.9) easily obtained, for more details see ^[16,24].

4.   Numerical scheme

Consider the generalized space fractional partial differential equations of the type given in Eq (1.1) with their given conditions. In order to use the Chebyshev collocation method, let us approximate $u(x, t)$ as follows ^[26,27,28]:

substituting (4.1) in (1.1) we obtain:

with the help of Theorem 1 then:

Now, we turn to collocate equation (4.3) at $(N+1)$ points, the collocation points are defined in the following form:

By substituting the collocation points (4.4) in (4.3) we get:

Also, two additional equations may generate from the boundary conditions using relation (4.1) in (1.3) as:

The collocated equation (4.5), together with the generated equations of the boundary conditions (4.6), give us an ordinary system of differential equations with $\phi_k(t)$ as the unknowns, which can be solved by a suitable technique. Using the initial conditions (1.2) and by the help of relation (4.1) and the orthogonality (2.8) we can generate initial conditions for the proposed system of differential equations, the IC may take the form:

consequently, by expanding $h(x)$ in terms of $C_k(x)$ and comparing the coefficents of Eq (4.7), then, we can find the constants $\phi_k(0)$. The produced system of ordinary differential equations according to (4.5) is linear and generally has the following matrix form:

where

and $\bar{Q}$ is a square constant matrix represent the coefficients of the unknowns $\phi_k(t)$, which is featured by the first column is null. Additionally, (4.8) may written as:

now, the system (4.9) is ready to solve with a suitable solver technique with the subjected initial conditions (4.7).

5.   Numerical applications

In this section, several numerical applications (physical models) have been given to illustrate the accuracy and effectiveness of the method.

Example 1:

Consider the following space fractional order PDE:

The IC is:

and the BCs:

Equation (5.1) is obtained when $\gamma_0\neq0, \; \; \gamma_k = 0,$ in Eq (1.1), and $0 < \gamma_0 < 1$ where the exact solution of Eq (5.1) under conditions (5.2) and (5.3) is $u(x, t) = x^2(t+1)$, with $Q_{0}(x) = P = 1$ and the function $f(x, t)$ is $(1.91116+1.91116 t) x^{1.1}+x^2$ at $\gamma_0 = 0.9$. At $N = 3$ according to (4.1) we have:

By the same proccess as: Eq (4.2) to (4.9) we have:

and by expanding $h(x) = x^2$ in terms of $C_k(x)$ according to (2.8) and comparing the coefficients then, we get the initial conditions of the differential system as:

Two additional equations may generate from the boundary conditions (5.3) using relation (4.1) in (5.3), then:

System (4.9) with matrices (5.5) and the initial conditions (5.6) is a system of differential equations, (Eq (5.7) may replace with the last two equations in (4.9)) the Runge-Kutta method of the fourth-order (RK4) is used here with $h$ step size equal to $0.01$ with $100$ iterations means that $0\leq t\leq 1$, (the regular algorithm for RK4 is coded by the authors using Mathematica.10. package) the numerical results obtained as:

According to (5.4) one obtains the approximate solution $u_3(x, 1)$ (at $t = 1$) using the last row in (5.8) as:

As references ^[27,28,29], their numerical results were obtained using finite difference method (FDM) for the differential system, we turn to solve the system (4.9) with matrices (5.5) using FDM. Then,

Therefore, the system in Eq (4.9) with matrices (5.5), is discretized in the time and have the following form:

or

where

Hense, a sample of the numerical results for FDM obtained as:

In Table 1 the comparison of the absolute errors for the present method with both RK4 and FDM at $N = 3, \; \; \Delta t = 0.01$ where $\gamma_0 = 0.9$, also, shows the numerical values of the approximate solution using the proposed method (using both RK4 and FDM) with the exact solution. Also, Table 2 shows the $L_2$ error norm^[26] at $N = 3$ at different values of $T$. In Figures 1 and 2 the comparison of the exact and the approximate solutions with both RK4 and FDM methods for example 1 with $N = 3$ and $T = 1, \; 2$.

Example 2:

Consider the following generalized space fractional order diffusion equation of the following type:

if $1 < \gamma_1 < 2$, at $Q_{1}(x) = -\Gamma(1.2)x^{1.8}, \; P = 1$, $f(x, t) = -3x^2(-1+2x)e^{-t},$ then, Eq (5.13) has the exact solution of of the form $u(x, t) = x^2(1-x)e^{-t}$ at $\gamma_1 = 1.8$, which mentiented in ^[27,28,29]. The IC is:

and the BCs:

At $N = 3$, according to (4.1), (using same prosses (4.2)–(4.9)), we have:

and $C$ not changed for $N = 3$ as example 1. In addition, by expanding $h(x) = x^2(1-x)$ in terms of $C_k(x)$ according to (2.8) and comparing the coefficents then, we get the initial conditions of the differential system as:

The generated equations from the homogenuous boundary conditions (5.15) using relation (4.1) are:

System (4.9) with matrcies (5.16) and the initial conditions (5.17) is a system of differential equations, by repleceing equtions (5.18) with the last two equations in (4.9) the RK4 method used as example 1 with $0\leq t\leq 2$. The RK4 method's numeric results at $t = 1, \; \; t = 2, \; N = 3$ obtained as:

As example 1 we turn to solve the system (4.9) with matrices (5.16) using FDM. Then, using same process as (5.10), (5.11) the results are obtained. In Table 3 the comparison of the absolute errors for the present two schemes at $N = 3, \; \; T = 2$ with the methods mentioned in ^[27,28,29]. Also, the numeric absolute errors are represent in Table 3 for the collocation method with Chebyshev first ^[29] second ^[27] and third ^[28] kinds. These values show that the fifth kind gives a more accurate approximate solution using the proposed method with RK4, but less accuracy is given when using regular FDM with the present method. Table 4 shows the $L_2$ error norm at $N = 3$ at two values of $T$. In Figures 3 and 4 the comparison of the exact and the approximate solutions with both RK4 and FD methods for example 2 with $N = 3$ and $T = 1, \; 2$.

Example 3:

Consider the following space fractional-order advection-dispersion equation of the following type:

if $1 < \gamma_1 < 2$ and $0 < \gamma_0 < 1$ at $Q_{1}(x) = -1$, $Q_{0}(x) = 1$, $P = 1$ and $f(x, t) = e^{-2 t} \left(-2 \left(x^{\gamma_1 }-x^{\gamma_0 }\right)-(\Gamma(\gamma_1 +1)+\Gamma(\gamma_0 +1))+\frac{\Gamma(\gamma_1 +1)}{\Gamma(1-\gamma_0 +\gamma_1) }x^{\gamma_1 -\gamma_0 }\right)$, then, Eq (5.20) has the exact solution of the form $u(x, t) = (x^{\gamma_1 }-x^{\gamma_0 })e^{-2 t}$, this case mentiented in ^[30,31,32] with $\gamma_1 = 2, \gamma_0 = 1$, where, the IC is:

and the BCs are homogenuous as:

At $N = 3, \gamma_1 = 2, \gamma_0 = 1$, according to (4.1), (using same prosses (4.2)–(4.9)), we have:

and $C$ not changed for $N = 3$ as examples 1 and 2. In addition, by expanding $h(x) = x^{\gamma_1 }-x^{\gamma_0 }$ in terms of $C_k(x)$ according to (2.8) and comparing the coefficients then, we get the initial conditions of the differential system at $N = 3, \gamma_1 = 2, \gamma_0 = 1$ as:

The generated equations from the homogenuous boundary conditions (5.22) are same as (5.18) using relation (4.1) in example 2. The system (4.9) with matrices (5.23) and the initial conditions (5.24) is a system of differential equations, by replacing the generated equations from the homogenous boundary conditions with the last two equations in (4.9), the RK4 method may be used as examples 1 and 2 with $0\leq t\leq 2$. As references ^[31,32] the numerical results obtained using FDM except ^[30] used the non-standard FDM for the differential system. As examples 1 and 2 we turn to solve the system (4.9) with matrices (5.23) using FDM. Then we use same elements as example 2, as system (5.10), (5.11) but using matrices (5.23). In Table 5 the comparison of the absolute errors for the present method (using the two proposed schemes) at $N = 3,$ where $\gamma_1 = 2, \; \gamma_0 = 1, \; T = 2$ with the methods mentioned in ^[30,31,32]. Also, it shows the numerical values of the proposed method gives best approximate solution except ^[30] which uses a modified technique (the non-standard FDM with Vieta-Lucas polynomials), where ^[31] uses Legendre polynomials FDM and ^[32] uses fourth kind Chebyshev polynomials with FDM. Table 6 gives the $L_2$ error norm along the interval $[0, 1]$ at $N = 3$ with two values of $T$. In Figures 5 and 6 the comparison of the exact and the approximate solutions with both RK4 and FD methods for example 3 with $N = 3$ and $T = 1, \; 2$.

Example 4:

Consider the following space fractional-order advection-dispersion equation, semilar to example 3, but $\gamma_1 = 1.5, \gamma_0 = 1$ which found at ^[30,32,33]:

with $Q_{1}(x) = -1$, $Q_{0}(x) = 2$, $P = 1$ and $f(x, t) = -\frac{4 (-1+t) t \sqrt{x}}{\sqrt{\pi }}+(-1+2 t) (-1+x) x+2 (-1+t) t (-1+2 x)$, then, Eq (5.25) has the exact solution of of the form $u(x, t) = \left(x^2-x\right) \left(t^2-t\right)$, where, the IC is homogenuous as:

also, the BCs are homogenuous as:

Equation (5.25) according to (4.1), by using the same prosses (4.2)–(4.9) where $C$ not changed for $N = 3$ as the pervuose examples, we have:

Additionally, by the homogenety of the IC, then, we get zero initial conditions of the differential system as:

The generated equations from the homogenous boundary conditions (5.27) are the same as given in examples 2 and 3. The system (4.9) with matrices (5.28) has zero ICs, by replacing the generated equations from the homogenous boundary conditions with the last two equations in (4.9), the RK4 used as examples 2, 3 with $0\leq t\leq 2$. As ref ^[30] the numerical results were obtained using the non-standard FDM for the differential system with the aid of Vieta-Lucas polynomials. Therefore, as example 3 we turn to solve the system (4.9) with matrices (5.23) using FDM. Then we use same elements as examples 2, 3, for systems (5.10), (5.11) but using matrices (5.28). The numerical comparisons will hold only with ^[30] because the results in ^[32,33] (collocation method with fourth and second Chebyshev kinds) are less than $10^{-5}$, it is much less accurate than indicated in our results. In Table 7 the comparison of the absolute errors for the present two schemes (PM with RK4 and FDM) at $N = 3,$ where $\gamma_1 = 1.5, \; \gamma_0 = 1, \; T = 0.5$ with ^[30], while same comparison given in Table 8 but $T = 0.5$. Also, it shows the numerical values of the proposed method gives a highly accurate approximate solution with RK4, and ^[30] which uses a modified technique gives accuracy near PM with FDM. Table 9 gives the $L_2$ error norm along the interval $[0, 1]$ at $N = 3$ with three values of $T$. In Figures 7–9 the comparison of the exact and the approximate solutions with both RK4 and FD methods for example 1 with $N = 3$ and $T = 0.3, \; 0.5, \; 0.9$. In the end, we conclude that the Chebyshev fifth-kind series approximation gives a great accuracy when using high appropriate accurate methods, and the Runge-Kutta method remains one of the best methods in dealing with linear systems, as was shown in the last two examples.

6.   Conclusions

A numerical study for a generalized form of linear space-fractional partial differential equations is introduced using the Chebyshev fifth kind series. The suggested general form represents many fractional-order mathematical physics models, as advection-dispersion equation and diffusion equation. Additionally, the proposed scheme uses the Shifted Chebyshev polynomials of the fifth-kind, where the fractional derivatives are expressed in terms of Caputo's definition. Therefore, the collocation method is used to reduce the GFPDE to a system of ordinary differential equations which can be solved numerically. In addition, the classical fourth-order Runge-Kutta method is used to treat the differential system as well as the finite difference method which obtains a great accuracy. We have presented many numerical examples, where represent mathematical physical models, that greatly illustrate the accuracy of the presented study to the proposed GFPDE, and also show how that the fifth-kind polynomials are very competitive than others.

Acknowledgments

The authors are thankful to the Taif University (supporting project number TURSP-2020/160), Taif, Saudi Arabia.

Conflict of interest

The authors declare that they have no competing interests.