A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel

Mohamed Obeid; Mohamed A. Abd El Salam; Mohamed S. Mohamed; Mohamed Obeid; Mohamed A. Abd El Salam; Mohamed S. Mohamed

doi:10.3934/math.2023855

AIMS Mathematics

2023, Volume 8, Issue 7: 16724-16747. doi: 10.3934/math.2023855

Previous Article Next Article

Research article Special Issues

A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel

1.
Department of Mathematics, Faculty of Science, Al-Azhar University, Nasr City 11884, Cairo, Egypt
2.
Basic Science Department, October High Institute for Engineering and Technology, Giza, Egypt
3.
Department of Mathematics, College of Science, Taif University, P.O. Box 11099, Taif, 21944, Saudi Arabia

Received: 18 February 2023 Revised: 13 April 2023 Accepted: 27 April 2023 Published: 12 May 2023
MSC : 35G05, 65N35

Polynomial based numerical techniques usually provide the best choice for approximating the solution of fractional differential equations (FDEs). The choice of the basis at which the solution is expanded might affect the results significantly. However, there is no general approach to determine which basis will perform better with a particular problem. The aim of this paper is to develop a novel generalized symmetric orthogonal basis which has not been discussed in the context of numerical analysis before to establish a general numerical treatment for the FDEs with a singular kernel. The operational matrix with four free parameters was derived for the left-sided Caputo fractional operator in order to transform the FDEs into the corresponding algebraic system with the aid of spectral Galerkin method. Several families of the existing polynomials can be obtained as a special case from the new basis beside other new families generated according to the value of the free parameters. Consequently, the operational matrix in terms of these families was derived as a special case from the generalized one up to a coefficient diagonal matrix. Furthermore, different properties relevant to the new generalized basis were derived and the error associated with function approximation by the new basis was performed based on the generalized Taylor's formula.

Keywords:

Citation: Mohamed Obeid, Mohamed A. Abd El Salam, Mohamed S. Mohamed. A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel[J]. AIMS Mathematics, 2023, 8(7): 16724-16747. doi: 10.3934/math.2023855

Related Papers:

[1]	Imran Talib, Md. Nur Alam, Dumitru Baleanu, Danish Zaidi, Ammarah Marriyam . A new integral operational matrix with applications to multi-order fractional differential equations. AIMS Mathematics, 2021, 6(8): 8742-8771. doi: 10.3934/math.2021508
[2]	Zhengang Zhao, Yunying Zheng, Xianglin Zeng . Finite element approximation of fractional hyperbolic integro-differential equation. AIMS Mathematics, 2022, 7(8): 15348-15369. doi: 10.3934/math.2022841
[3]	Manal Alqhtani, Khaled M. Saad . Numerical solutions of space-fractional diffusion equations via the exponential decay kernel. AIMS Mathematics, 2022, 7(4): 6535-6549. doi: 10.3934/math.2022364
[4]	Ahmed F. S. Aboubakr, Gamal M. Ismail, Mohamed M. Khader, Mahmoud A. E. Abdelrahman, Ahmed M. T. AbdEl-Bar, Mohamed Adel . Derivation of an approximate formula of the Rabotnov fractional-exponential kernel fractional derivative and applied for numerically solving the blood ethanol concentration system. AIMS Mathematics, 2023, 8(12): 30704-30716. doi: 10.3934/math.20231569
[5]	Hailun Wang, Fei Wu, Dongge Lei . A novel numerical approach for solving fractional order differential equations using hybrid functions. AIMS Mathematics, 2021, 6(6): 5596-5611. doi: 10.3934/math.2021331
[6]	Waleed Mohamed Abd-Elhameed, Ahad M. Al-Sady, Omar Mazen Alqubori, Ahmed Gamal Atta . Numerical treatment of the fractional Rayleigh-Stokes problem using some orthogonal combinations of Chebyshev polynomials. AIMS Mathematics, 2024, 9(9): 25457-25481. doi: 10.3934/math.20241243
[7]	Martin Stoll, Hamdullah Yücel . Symmetric interior penalty Galerkin method for fractional-in-space phase-field equations. AIMS Mathematics, 2018, 3(1): 66-95. doi: 10.3934/Math.2018.1.66
[8]	Chang Phang, Abdulnasir Isah, Yoke Teng Toh . Poly-Genocchi polynomials and its applications. AIMS Mathematics, 2021, 6(8): 8221-8238. doi: 10.3934/math.2021476
[9]	Xiaojun Zhou, Yue Dai . A spectral collocation method for the coupled system of nonlinear fractional differential equations. AIMS Mathematics, 2022, 7(4): 5670-5689. doi: 10.3934/math.2022314
[10]	Dongjie Gao, Peiguo Zhang, Longqin Wang, Zhenlong Dai, Yonglei Fang . A novel high-order symmetric and energy-preserving continuous-stage Runge-Kutta-Nyström Fourier pseudo-spectral scheme for solving the two-dimensional nonlinear wave equation. AIMS Mathematics, 2025, 10(3): 6764-6787. doi: 10.3934/math.2025310

Abstract

1. Introduction

Fractional modeling has demonstrated its effectiveness in describing various time-evolving dynamical systems ^[1,2,3]. Over the years, researchers have developed numerous numerical techniques to approximate the solutions of these fractional models ^[4,5,6,7,8]. Among these methods, the most commonly used techniques involve orthogonal functions, in which the solution is expanded as a weighted infinite sum of basis functions ^[9]. By truncating the polynomial series, the fractional model under consideration can be transformed into a system of algebraic equations. This transformation is achieved by constructing a matrix representation of the differential operator, known as an operational matrix. Spectral methods can be employed to address the residual associated with the fractional model, aiding in the derivation of the algebraic system that governs the unknowns in the finite weighted sum of the solution.

However, most of the developed numerical schemes depend on the non generalized orthogonal polynomials such as Legendre ^[10], Laguerre ^[11], Hermite ^[12], and different kinds of Chebyshev polynomials ^[13,14,15]. For example, Secer and Altun ^[16] utilized Legendre wavelets operational matrix technique to solve a system of fractional differential equations (FDEs). Baishya and Veeresha ^[17] construed a numerical scheme based on Lagurre polynomials operational matrix together with spectral collocation method for solving FDE with Mittag-Leffler kernel. Tural-Polat and Dincel ^[14] used the third-kind Chebyshev polynomials along with spectral collocation method to approximate the solution of multi-term variable order FDEs. In ^[18], Abd-Elhameed and Youssri used the modified tau method with the shifted Chebyshev polynomials of the fifth-kind to examine the solution of FDEs. Sadri and Aminikhah developed a numerical scheme based on the shifted Chebyshev polynomials of the sixth-kind operational matrix for studying a class of delay fractional order partial differential equations (FPDEs) ^[19]. Due to the mathematical difficulties, relatively few studies involve the generalized families of orthogonal polynomials, such as Jacobi and Ultraspherical polynomials (also known as Gegenbauer polynomials), compared with the researches involving the traditional polynomial ^{[20,21,22,23,24,25]}. This deficiency may affect the development of accurate numerical schemes by limiting the number of basis available for certain applications.

Recently, Masjed-Jamei ^[26] introduced a more general family of orthogonal polynomials with four free parameters known as the basic class of symmetric orthogonal polynomials (BCSPs) due to its symmetry property. This family not only includes a vast number of traditional orthogonal polynomials even the generalized one it also can be used to generate new classes of polynomials. The development of the BCSPs started with studying a generalization of the regular Sturm-Liouville problem known as:

$\begin{equation} \frac{d}{dt}\left(p\left(t\right)\frac{dy_{i}}{dt}\left(t\right)\right)+\left(\lambda_{i}q\left(t\right)-w\left(t\right)\right)y_{i}\left(t\right) = 0, \quad p\left(t\right) > 0, \ q\left(t\right) > 0, \end{equation}$

(1.1)

which, defined on the open interval $\left(a, b\right)$ and satisfies certain boundary conditions. Whereas for any two eigenfunctions $y_{i}\left(t\right), y_{j}\left(t\right)$ , the following orthogonality relation holds

$\begin{equation} \int_{\left(a, b\right)}w\left(t\right)y_{i}\left(t\right)y_{j}\left(t\right)dt = \left(\int_{\left(a, b\right)}w\left(t\right)y^{2}_{i}\left(t\right)dt\right)\delta_{i, j}, \end{equation}$

(1.2)

where $\delta_{i, j}$ is the Kronecker function. Many non generalized orthogonal polynomials that have been mentioned earlier are eigenfunction of Eq (1.1). Considering the symmetric functions $\Theta_{i}\left(t\right)$ satisfying $\Theta_{i}\left(-t\right) = \left(-1\right)^{i}\Theta_{i}\left(t\right)$ , the regular Sturm-Liouville equation can be generalized to the form:

$\begin{equation} \mathcal{A}\left(t\right)\Theta_{i}^{\prime\prime}\left(t\right)+\mathcal{B}\left(t\right)\Theta_{i}^{\prime}+\left(\lambda_{i}\mathcal{C}\left(t\right)+\mathcal{D}\left(t\right)+\frac{1-\left(-1\right)^{i}}{2}\mathcal{E}\left(t\right)\right)\Theta_{i}\left(t\right) = 0, \end{equation}$

(1.3)

where $\mathcal{A}\left(t\right), \mathcal{C}\left(t\right), \mathcal{D}\left(t\right)$ and $\mathcal{E}\left(t\right)$ are even arbitrary functions, and $\mathcal{B}\left(t\right)$ is an odd function. The functions $\Theta_{i}\left(t\right)$ obey the orthogonality relation with respect to the weight function $W\left(t\right)$ over the symmetric closed interval $\left[-\theta, \theta\right]$ :

$\begin{equation} \int_{\left[-\theta, \theta\right]}W\left(t\right)\Theta_{i}\left(t\right)\Theta_{j}\left(t\right)dt = \left(\int_{\left[-\theta, \theta\right]}W\left(t\right)\Theta_{i}^{2}dt\right)\delta_{i, j}, \end{equation}$

(1.4)

where the function $W\left(t\right)$ is defined as:

$\begin{equation} W\left(t\right) = \mathcal{C}\left(t\right)\cdot\exp\left(\int\frac{\mathcal{B}\left(t\right)-\mathcal{A}^{\prime}\left(t\right)}{\mathcal{A}\left(t\right)}dt\right). \end{equation}$

(1.5)

As a special case of (1.3), the generalized BCSPs with four free parameters, will denoted as $\mathcal{S}_{i}^{\left(r, s, p, q\right)}\left(t\right)$ through this paper, can be constructed, whereas the arbitrary function in (1.3) takes the form ^[26]:

$\begin{align*} &\mathcal{A}\left(t\right) = t^{2}\left(p t^{2}+q\right), \ \mathcal{B}\left(t\right) = t\left(r t^{2}+s\right), \\&\mathcal{C}\left(t\right) = t^{2}, \ \mathcal{D}\left(t\right) = 0, \mathcal{E}\left(t\right) = -s; \end{align*}$

where $r, s, p$ and $q$ are free parameters.

However, since the time BCPs were introduced, it had never been used in any numerical techniques. The main aim of this study is to establish a novel shifted subclass of this family applied to a spectral Galerkin scheme for examining the numerical solution of the multi-term FDE of the form:

$\begin{equation} \sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right)\frac{d^{\alpha_{k}}u\left(t\right)}{dt^{\alpha_{k}}} = f\left(t, u\left(t\right)\right), \quad t\in\left[0, T\right], \ T\in\mathbb{R}, \ N\in\mathbb{N}_{0}, \end{equation}$

(1.6)

where the fractional differentiation of order, $\alpha_{k}\in\mathbb{R}^{+}$ , is defined in Caputo sense, $u\left(t\right)$ satisfies the boundary conditions $u\left(0\right) = u_{0}, u\left(T\right) = u_{1}$ and $\gamma_{N}\left(t\right)\neq0$ . Moreover, the proposed Caputo operational matrix introduced in the numerical scheme is a novel generalized matrix from which we will inferred different operational matrices with respect to other families of orthogonal polynomials.

The rest of this study is organized as follows: The next section provides the relevant definitions of fractional calculus used in this paper. In Section 3, we begin the development of the shifted family of polynomials after illustrating the main results of the BCSPs. After that, Section 4 provides the generalized operational matrix, also, the convergence analysis will be introduced in Section 5. Finally, in Section 6 the numerical scheme will be discussed, supported with different numerical examples.

2. Preliminaries

In this section, some basic definition of fractional calculus is listed which will be used further in this paper.

Definition 2.1. Let $g\left(t\right):\left[a, b\right]\rightarrow \mathbb{R}$ , $a < b$ , be Lebesgue square integrable function, and let $\alpha$ be the order of integration. The Riemann-Liouville fractional integral operator is defined as ^[27]:

$\begin{equation*} \mathcal{I}^{\alpha}_{c}g\left(t\right) = \left\{ {\begin{array}{*{20}{l}} \frac{1}{\Gamma\left(\alpha\right)}\int_{c}^{t}\frac{g\left(\tau\right)}{\left(t-\tau\right)^{1-\alpha}}d\tau, & \alpha > 0, \\ g\left(t\right), & \alpha = 0. \end{array}} \right. \end{equation*}$

Definition 2.2. The left Caputo's fractional derivative with singular kernel of order $\alpha > 0, \alpha\notin\mathbb{N}_{0}$ for the function $g\left(t\right):\left[a, b\right]\rightarrow \mathbb{R}$ is defined as:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}g\left(t\right) = \frac{1}{\Gamma\left(n-\alpha\right)}\int_{c}^{t}\frac{g^{\left(n\right)}\left(\tau\right)}{\left(t-\tau\right)^{1+\alpha-n}}d\tau, \end{equation}$

(2.1)

for $n = -\lceil-\alpha\rceil$ , ^[27].

Remark 2.1. The left Caputo's fractional derivative of the power function $t^{n}$ , for $n\in\mathbb{N}_{0}$ , has the following two cases:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}t^{n} = \left\{ {\begin{array}{*{20}{l}} 0, & for\ n < \lceil\alpha\rceil, \\ \frac{\Gamma\left(n+1\right)}{\Gamma\left(n+1-\alpha\right)}t^{n-\alpha}, & for\ n\geq\lceil\alpha\rceil. \end{array}} \right. \end{equation}$

(2.2)

Similar to normal derivatives, the left sided Riemann-Liouville and Caputo fractional derivatives obey the linearity relation:

$\begin{equation} {^{RL}}{\mathcal{D}}^{\alpha}_{a+}\sum\limits_{s = 0}^{m}b_{s}g_{s}\left(t\right) = \sum\limits_{s = 0}^{m}b_{s} {^{RL}}{\mathcal{D}}^{\alpha}_{a+}g_{s}\left(t\right), \end{equation}$

(2.3)

and

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{a+}}\sum\limits_{s = 0}^{m}b_{s}g_{s}\left(t\right) = \sum\limits_{s = 0}^{m}b_{s} {^C}{\mathcal{D}^{\alpha}_{a+}}g_{s}\left(t\right), \end{equation}$

(2.4)

where $\left\{b_{s}\right\}_{s = 0}^{m}$ are constants.

The following generalized Taylor's formula was introduced by Odibat and Momani in ^[28], which have the following statement, (see also ^[29,30]).

Definition 2.3. Suppose that $g\left(t\right)$ is differentiable, and ${^C}{\mathcal{D}^{k\alpha}_{0+}}g\left(t\right)\in \mathcal{C}\left(0, 1\right]$ for $k = 0, 1, \cdots, m$ , where $0 < \alpha\leq 1$ , the generalized Taylor's formula takes the form:

$\begin{equation} g\left(t\right) = \sum\limits_{i = 0}^{m-1}\frac{t^{i\alpha}}{\Gamma\left(i\alpha+1\right)}{^C}{\mathcal{D}^{\alpha}_{0+}}g\left(0^{+}\right)+\frac{t^{m\alpha}}{\Gamma\left(m\alpha+1\right)}{^C}{\mathcal{D}^{\alpha}_{0+}}g\left(\xi\right), \end{equation}$

(2.5)

with $0\leq \xi\leq t$ for all $t\in\left(0, 1\right]$ .

3. Basic symmetric polynomials class

This section is concerned with recalling the relations relevant to the basic class of symmetric polynomials (BCSPs) along with developing the main concepts of the novel shifted subclass.

3.1. BCSPs properties

As mentioned earlier, the BSCPs are the eigenfunctions of the symmetric generalization of the Sturm-Liouville problem in (1.3) of the form:

$\begin{array}{l} t^{2}\left(pt^{2}+q\right)\mathcal{S}_{i}^{\prime\prime^{\left(r, s, p, q\right)}}\left(t\right)+t\left(rt^{2}+s\right)\mathcal{S}_{i}^{\prime^{ \left(r, s, p, q\right)}}\left(t\right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\left[i\left(r+\left(i-1\right)p\right)t^{2}+\left(1-\left(-1\right)^{i}\right)\frac{s}{2}\right]\mathcal{S}_{i}^{\left(r, s, p, q\right)}\left(t\right) = 0. \end{array}$

(3.1)

The BCSPs can be defined over the symmetric interval $\left[-\rho, \rho\right]$ , which fulfill certain conditions that will be mentioned, through the following analytic form.

Definition 3.1. Suppose that $t\in\left[-\rho, \rho\right]$ , and $i\in\mathbb{N}_{0}$ . For the parameters $r, s, p, q\in\mathbb{R}$ such that $s$ and $q$ nor $r$ and $p$ can vanish together, the BCSPs of degree $i$ , can be defined as (see ^[26]):

$\begin{equation} \mathcal{S}_{i}^{\left(r, s, p, q\right)}\left(t\right) = \sum\limits_{k = 0}^{\lfloor \frac{i}{2}\rfloor}\binom{\lfloor \frac{i}{2}\rfloor}{k}\left(\prod\limits_{\mu = 0}^{\lfloor \frac{i}{2}\rfloor-\left(k+1\right)}\xi_{\mu, i}^{r, s, p, q}\right)t^{i-2k}, \end{equation}$

(3.2)

where $\xi_{\mu, i}^{r, s, p, q}$ is given by:

$\begin{equation} \xi_{\mu, i}^{r, s, p, q} = \frac{\left(2\mu+\left(-1\right)^{i+1}+2\lfloor \frac{i}{2}\rfloor\right)p+r}{\left(2\mu+\left(-1\right)^{i+1}+2\right)q+s}. \end{equation}$

(3.3)

In the purpose of defining the monic type of this family, the leading term is multiplied by the inverse of the coefficient $\xi_{\mu, i}^{r, s, p, q}$ as follows ^[26].

Definition 3.2. The monic type of BCSPs (MBCSPs) is defined through the following representation:

$\begin{equation} \bar{\mathcal{S}}_{i}^{\left(r, s, p, q\right)}\left(t\right) = \left(\prod\limits_{\nu = 0}^{\lfloor \frac{i}{2}\rfloor-1}\frac{1}{\xi_{\nu, i}^{r, s, p, q}}\right)\mathcal{S}_{i}^{\left(r, s, p, q\right)}\left(t\right). \end{equation}$

(3.4)

indicates the graph of $\bar{\mathcal{S}}_{5}^{\left(r, s, p, q\right)}\left(t\right)$ for different values of $s$ , with a fixed $r, q$ and $p$ .

Figure 1. The graph of MBCSP with parameter

$r = -7, p = -1 = -q$ and

$s\in\left\{2, 4, 6, 8\right\}$ .

DownLoad: Full-Size Img PowerPoint

By the means of the following three-term recurrence relation, the monic type of BSCPs can be generated as:

$\begin{equation} \bar{\mathcal{S}}_{i+1}^{\left(r, s, p, q\right)}\left(t\right) = t\bar{\mathcal{S}}_{i}^{\left(r, s, p, q\right)}\left(t\right)+\varrho_{i}^{r, s, p, q}\bar{\mathcal{S}}_{i-1}^{\left(r, s, p, q\right)}\left(t\right), \end{equation}$

(3.5)

with the leading members $\bar{\mathcal{S}}_{0}^{\left(r, s, p, q\right)}\left(t\right) = 1, \bar{\mathcal{S}}_{0}^{\left(r, s, p, q\right)}\left(t\right) = t$ , and $\varrho_{i}^{r, s, p, q}$ is defied by:

$\begin{equation} \varrho_{i}^{r, s, p, q} = \frac{\left(2p q\right)i^{2}+2\left(\left(r-2p\right)q-\left(-1\right)^{i}p s\right)i+\left(r-2p\right)s\left(1-\left(-1\right)^{i}\right)}{2\left(2p i+r-p\right)\left(2 p i+r-3 p\right)}. \end{equation}$

(3.6)

According to Favard's theorem ^[31], the above three-term recurrence relation implies the next orthogonality property ^[26]:

$\begin{array}{l} \langle\bar{\mathcal{S}}_{i}^{\left(r, s, p, q\right)}\left(t\right), \bar{\mathcal{S}}_{j}^{\left(r, s, p, q\right)}\left(t\right)\rangle = \int_{-\rho}^{\rho}\Xi^{\left(r, s, p, q\right)}\left(t\right)\bar{\mathcal{S}}_{i}^{\left(r, s, p, q\right)}\left(t\right)\bar{\mathcal{S}}_{j}^{\left(r, s, p, q\right)}\left(t\right)dt\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \left(\left(-1\right)^{i}\prod\limits_{\nu = 0}^{i}\varrho_{\nu}^{r, s, p, q}\int_{-\rho}^{\rho}\Xi^{\left(r, s, p, q\right)}\left(t\right)dt\right)\delta_{i, j}, \end{array}$

(3.7)

where, $\varrho_{\nu}^{r, s, p, q}$ as in (3.6), $\delta_{i, j}$ is the Kronecker function, and the weight function $\Xi^{\left(r, s, p, q\right)}\left(t\right)$ is defined as:

$\begin{equation} \Xi^{\left(r, s, p, q\right)}\left(t\right) = \exp\left(\int\frac{s+t^{2}\left(r-2p\right)}{t\left(q+pt^{2}\right)}\right), \end{equation}$

(3.8)

although, the expression $\left(pt^{2}+q\right)\Xi^{\left(r, s, p, q\right)}\left(t\right)$ should vanishes at the boundaries of the orthogonality interval $\left[-\rho, \rho\right]$ .

3.2. The shifted class of MBCSPs

In this part of the present article, we introduce the shifted class of the MBCSPs defined over the interval $\Omega = \left[0, \rho\right]$ by means of a certain transformation. To the best of our knowledge, this class of polynomials are being discussed here for the first time. The shifted class will be denoted by $\mathcal{G}_{i}^{\left(r, s, p, q\right)}\left(t\right)$ , for $i\in\mathbb{N}_{0}$ .

Definition 3.3. The shifted class of the MBCSPs of order $i\in\mathbb{N}_{0}$ , defined over the interval $\Omega$ , is related to the MBCSPs via the following transformation:

$\begin{equation} \mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right) = \bar{\mathcal{S}}_{i}^{\left(r, s, p, q\right)}\left(2t-\rho\right), \end{equation}$

(3.9)

where, $\rho\in\mathbb{R}$ .

The set $\left\{\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)\right\}_{i = 0}^{\infty}$ forms a complete set of orthogonal functions with respect to the weight function $\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)$ defined as:

$\begin{equation} \mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right) = \exp\left(\int\frac{s+t^{2}\left(r-2p\right)}{t\left(q+pt^{2}\right)}dt\right)\Big|_{2t-\rho}, \end{equation}$

(3.10)

which obeys the inner product of the form:

$\begin{array}{l} \langle \mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right), \mathcal{G}_{j, \rho}^{\left(r, s, p, q\right)}\left(t\right)\rangle = \int_{\Omega}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)\mathcal{G}_{j, \rho}^{\left(r, s, p, q\right)}\left(t\right)dt\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;= \left(\left(-1\right)^{i}\prod\limits_{\nu = 0}^{i}\varrho_{\nu}^{r, s, p, q}\int_{\Omega}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)dt\right)\delta_{i, j}. \end{array}$

(3.11)

In order to establish the analytic form of the polynomials $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ , we introduce the series representation in the next lemma.

Lemma 3.1. Suppose that $t\in\Omega$ and $i\in\mathbb{N}_{0}$ , the polynomials $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ can be defined through the following series representation:

$\begin{equation} \mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right) = \sum\limits_{k = 0}^{\lfloor \frac{i}{2}\rfloor}\sum\limits_{l = 0}^{i-2k}h_{k, l, i}^{r, s, p, q}t^{l}, \end{equation}$

(3.12)

where

$\begin{equation} h_{k, l, i}^{r, s, p, q} = \binom{\lfloor \frac{i}{2}\rfloor}{k}\binom{i-2k}{l}\prod\limits_{\nu = 0}^{\lfloor \frac{i}{2}\rfloor-1}\frac{1}{\xi_{\nu, i}^{r, s, p, q}}\prod\limits_{\mu = 0}^{\lfloor \frac{i}{2}\rfloor-\left(k+1\right)}\xi_{\mu, i}^{r, s, p, q}\frac{\left(-1\right)^{i-l}}{2^{-l}\rho^{2k+l-i}}. \end{equation}$

(3.13)

Proof. According to the power series of $\mathcal{S}_{i}^{\left(r, s, p, q\right)}\left(t\right)$ in (3.2) and the connection between the monic class (3.4). By substituting the transfomation in (3.3), we arrive at the relation:

$\begin{equation*} \mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right) = \sum\limits_{k = 0}^{\lfloor \frac{i}{2}\rfloor}\binom{\lfloor \frac{i}{2}\rfloor}{k}\prod\limits_{\nu = 0}^{\lfloor \frac{i}{2}\rfloor-1}\frac{1}{\xi_{\nu, i}^{r, s, p, q}}\prod\limits_{\mu = 0}^{\lfloor \frac{i}{2}\rfloor-\left(k+1\right)}\xi_{\mu, i}^{r, s, p, q}\left(2t-\rho\right)^{i-2k}, \end{equation*}$

by applying the binomial theorem, we get the desired result. □

Considering the column vector of the shifted polynomials $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ defined by:

$\begin{equation} \mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = \left[\mathcal{G}_{0, \rho}^{\left(r, s, p, q\right)}\left(t\right), \mathcal{G}_{1, \rho}^{\left(r, s, p, q\right)}\left(t\right), \cdots, \mathcal{G}_{n, \rho}^{\left(r, s, p, q\right)}\left(t\right)\right]_{1\times\left(n+1\right)}^{T}, \end{equation}$

(3.14)

this vector has a matrix representation that will be derived in the next lemma.

Lemma 3.2. Let $\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ be the shifted MBCSPs vector, where $n\in\mathbb{N}_{0}$ . For $t\in\Omega$ , the vector $\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ can be written as:

$\begin{equation} \mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = Q^{r, s, p, q}_{\rho, n}\mathfrak{T}_{n}\left(t\right), \end{equation}$

(3.15)

where, $\mathfrak{T}_{n}\left(t\right) = \left[t^{0}, t^{1}, \cdots, t^{n}\right]^{T}$ , and $Q^{r, s, p, q}_{\rho, n}$ is the $\left(n+1\times n+1\right)$ lower triangle matrix and

$\begin{equation} Q_{\rho, n}^{r, s, p, q^{T}} = \begin{pmatrix} q^{r, s, p, q}_{\rho, 00}&q^{r, s, p, q}_{\rho, 10}&q^{r, s, p, q}_{\rho, 20}&\cdots&q^{r, s, p, q}_{\rho, n0}\\ 0&q^{r, s, p, q}_{\rho, 11}&q^{r, s, p, q}_{\rho, 21}&\cdots&q^{r, s, p, q}_{\rho, n1}\\ 0&0&q^{r, s, p, q}_{\rho, 22}&\cdots&q^{r, s, p, q}_{\rho, n2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&q^{r, s, p, q}_{\rho, nn} \end{pmatrix}, \end{equation}$

(3.16)

and, the entries $q^{r, s, p, q}_{\rho, ij}$ are given by:

$\begin{equation} q^{r, s, p, q}_{\rho, ij} = \sum\limits_{k = 0}^{\lfloor \frac{j}{2}\rfloor-\varrho_{i}}h_{k, i, j}^{r, s, p, q}, \end{equation}$

(3.17)

where, the $\varrho_{i}$ is defined as:

$\begin{equation} \varrho_{i} = \left\{ {\begin{array}{*{20}{l}} \lceil\frac{i}{2}\rceil & i\ even, \\ \lceil\frac{i}{2}\rceil-1 & i\ odd, \end{array}} \right. \end{equation}$

(3.18)

and, the coefficients $h_{k, i, j}^{r, s, p, q}$ as in (3.13).

Proof. By careful examination for the coefficients of the power form given in Lemma 3.1, the lemma can be proven. □

For instance, if $r = -5, s = 8, q = 1 = -p, n = 5,$ and $\rho = 1$ , the vector $\mathfrak{G}_{1, 5}^{\left(-5, 8, -1, 1\right)}\left(t\right)$ has the form:

$\begin{equation*} \mathfrak{G}_{1, 5}^{\left(-5, 8, -1, 1\right)}\left(t\right) = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 2 & 0 & 0 & 0 & 0 \\ \frac{1}{6} & -4 & 4 & 0 & 0 & 0 \\ -\frac{1}{8} & \frac{17}{4} & -12 & 8 & 0 & 0 \\ \frac{3}{80} & -\frac{12}{5} & \frac{92}{5} & -32 & 16 & 0 \\ -\frac{1}{40} & \frac{41}{20} & -22 & 68 & -80 & 32 \\ \end{pmatrix} \begin{pmatrix} t^{0}\\t^{1}\\t^{2}\\t^{3}\\t^{4}\\t^{5} \end{pmatrix}. \end{equation*}$

The essence advantage of working with the shifted MBCSPs that it includes as a special cases most of the traditional shifted orthogonal polynomials. The following propositions gives the relation between $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ and some of other families.

Proposition 3.1. The shifted Legendre polynomials over the interval $\left[0, 1\right]$ is defined in terms of the shifted MBCSPs as:

$\begin{equation} P_{i}^{\ast}\left(t\right) = \frac{\left(2i\right)!}{\left(i!\right)^{2}2^{i}}\mathcal{G}_{i, 1}^{\left(-2, 0, -1, 1\right)}\left(t\right). \end{equation}$

(3.19)

Proof. For the even case, taking $i = 2n$ in equation (3.4) and setting the parameters $r, s, p, q, \rho$ to be $-2, 0, -1, 1, 1$ respectively, we get

$\begin{equation*} \mathcal{G}_{2n, 1}^{\left(-2, 0, -1, 1\right)}\left(t\right) = \left(\prod\limits_{\nu = 0}^{n-1}\frac{2\nu+1}{-\left(2\nu+2n-1\right)-2}\right)\mathcal{S}_{2n}^{\left(-2, 0, -1, 1\right)}\left(2t-1\right). \end{equation*}$

Thus, the analytic form of the MBCSPs reads:

$\begin{equation*} \mathcal{G}_{2n, 1}^{\left(-2, 0, -1, 1\right)}\left(t\right) = \sum\limits_{k = 0}^{n}\frac{\left(-1\right)^{k}n!\Gamma\left(\frac{1}{2}+n\right)\left(\frac{3}{2}+n\right)_{-1-k+n}}{k!\left(n-k\right)!\Gamma\left(\frac{1}{2}+n-k\right)\left(\frac{3}{2}+n\right)_{n-1}}\left(2t-1\right)^{2n-2k}, \end{equation*}$

where $\left(s\right)_{k}$ is the Pochhammer symbol (for $s\in\mathbb{C}$ ) see ^[32]. Using the identities

$\begin{equation} \left(s\right)_{m} = \frac{\Gamma\left(s+m\right)}{\Gamma\left(s\right)}, \ \mathit{\text{and}}\ \Gamma\left(\frac{1}{2}+m\right) = \frac{\left(2m\right)!\sqrt{\pi}}{2^{2m}m!}, \end{equation}$

(3.20)

we get

$\begin{equation*} \mathcal{G}_{2n, 1}^{\left(-2, 0, -1, 1\right)}\left(t\right) = \sum\limits_{k = 0}^{n}\frac{\left(-1\right)^{k}\left(2n!\right)^{2}\left(4n-2k\right)!}{k!\left(4n\right)!\left(2n-k\right)!\left(2n-2k\right)!}\left(2t-1\right)^{2n-2k}, \end{equation*}$

multiplying both sides by $\frac{\left(4n\right)!}{\left(2n!\right)^{2}2^{2n}}$ , we have

$\begin{equation} \frac{\left(4n\right)!}{\left(2n!\right)^{2}2^{2n}}\mathcal{G}_{2n, 1}^{\left(-2, 0, -1, 1\right)}\left(t\right) = \sum\limits_{k = 0}^{n}\frac{\left(-1\right)^{k}\left(4n-2k\right)!}{2^{2n}k!\left(2n-k\right)!\left(2n-2k\right)!}\left(2t-1\right)^{2n-2k}, \end{equation}$

(3.21)

where the right hand side is the analytic form of the shifted Legendre polynomials $P_{2n}^{\ast}\left(t\right)$ . Similarly, the odd case can be obtained and the proof is complete. □

Proposition 3.2. The shifted Chebyshev polynomials of the first and second-kind over the interval $\left[0, 1\right]$ are defined in terms of the shifted MBCSPs as:

$\begin{equation} T_{i}^{\ast}\left(t\right) = 2^{i-1}\mathcal{G}_{i, 1}^{\left(-1, 0, -1, 1\right)}\left(t\right), \end{equation}$

(3.22)

$\begin{equation} U_{i}^{\ast}\left(t\right) = 2^{i}\mathcal{G}_{i, 1}^{\left(-3, 0, -1, 1\right)}\left(t\right). \end{equation}$

(3.23)

Proposition 3.3. The shifted generalized Ultraspherical polynomials over the interval $\left[0, 1\right]$ is defined in terms of the shifted MBCSPs as:

$\begin{equation} C_{i}^{\lambda^{\ast}}\left(t\right) = \frac{2^{i}\left(\lambda\right)_{i}}{i!}\mathcal{G}_{i, 1}^{\left(-2\lambda-1, 0, -1, 1\right)}\left(t\right). \end{equation}$

(3.24)

Also, another significant monic types of orthogonal polynomials can be obtained directly from the shifted MBCSPs.

Proposition 3.4. The shifted monic Chebyshev polynomials of the fifth and sixth-kind over the interval $\left[0, 1\right]$ are defined in terms of the shifted MBCSPs as:

$\begin{equation} \bar{X}_{i}^{\ast}\left(t\right) = \mathcal{G}_{i, 1}^{\left(-3, 2, -1, 1\right)}\left(t\right), \end{equation}$

(3.25)

$\begin{equation} \bar{Y}_{i}^{\ast}\left(t\right) = \mathcal{G}_{i, 1}^{\left(-5, 2, -1, 1\right)}\left(t\right). \end{equation}$

(3.26)

Remark 3.1. For obtaining the matrix representation of all other shifted orthogonal polynomials related to the shifted MBCSPs in terms of $Q_{\rho, n}^{r, s, p, q}$ the coefficients correlating these polynomials with $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ must be taken into account. As an example, by using of Proposition 3.2, the coefficient matrix of the shifted Chebyshev polynomials of the second-kind, for $n = 4$ , can be written as:

$\begin{equation} diag\left(2^{j}, 4\right)Q_{1, 4}^{-3, 0, -1, 1} = \begin{pmatrix} 2^{0} & 0 & 0 & 0 & 0 \\ 0 & 2^{1} & 0 & 0 & 0 \\ 0 & 0 & 2^{2} & 0 & 0 \\ 0 & 0 & 0 & 2^{3} & 0 \\ 0 & 0 & 0 & 0 & 2^{4} \\ \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ -1 & 2 & 0 & 0 & 0 \\ \frac{3}{4} & -4 & 4 & 0 & 0 \\ -\frac{1}{2} & 5 & -12 & 8 & 0 \\ \frac{5}{16} & -5 & 21 & -32 & 16 \\ \end{pmatrix}, \end{equation}$

(3.27)

where $diag\left(\cdot, n\right)$ denotes the $\left(n+1\right)\times\left(n+1\right)$ diagonal matrix.

According to the value of $r, s, p,$ and $q$ . The four parameter weight function in (3.10) reduces to the existing weight function of the related polynomials. For example, in the first-kind Chebyshev polynomials case, we have

$\begin{align*} \label{four_parameter} \mathcal{W}^{\left(-1, 0, -1, 1\right)}_{1}\left(t\right)& = \exp\left(\int\frac{t}{1-t^2}dt\right)\Big|_{2t-1}\\& = exp\left(-\frac{1}{2} \ln \left(1-\left(2t-1\right)^2\right)\right) \\& = \left(1-\left(2t-1\right)^{2}\right)^{\frac{-1}{2}}, \end{align*}$

which is the weight function of the shifted Chebychev polynomials of the first-kind.

Lemma 3.3. The shifted monic BCSPs has the following analytic form:

$\begin{equation} \mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right) = \sum\limits_{k = 0}^{i}q^{r, s, p, q}_{\rho, ik}t^{k}, \end{equation}$

(3.28)

where $q^{r, s, p, q}_{\rho, ij}$ , is the entry in (3.17).

Proof. Having the matrix representation in Lemma 3.2, the present lemma can be proven. □

Let $u\left(t\right)\in L^{2}\left(\Omega\right)$ , the function $u\left(t\right)$ can be expanded as a linear combination terms of the orthogonal basis $\left\{\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)\right\}_{i = 0}^{\infty}$ as follows:

$\begin{equation} u\left(t\right) = \sum\limits_{i = 0}^{\infty}c_{i}\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right), \end{equation}$

(3.29)

where $c_{i}$ are the coefficients of the series (3.29). By truncating the series at $i = n\in\mathbb{N}_{0}$ , the first $n+1$ terms are taken as an appropriate approximation of $u\left(t\right)$ as:

$\begin{equation} u\left(t\right)\simeq u_{n}\left(t\right) = \sum\limits_{i = 0}^{n}c_{i}\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right) = \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) \end{equation}$

(3.30)

where $\Lambda^{T} = \left[c_{0}, c_{1}, \cdots, c_{n}\right]^{T}$ is the coefficient vector having the entries:

$\begin{equation} c_{i} = \frac{1}{\sigma_{i, \rho}^{r, s, p, q}}\int_{\Omega}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)u\left(t\right)\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)dt, \end{equation}$

(3.31)

and the coefficients $\sigma_{i}^{r, s, p, q}$ are defined as:

$\begin{equation} \sigma_{i, \rho}^{r, s, p, q} = \left(-1\right)^{i}\prod\limits_{\nu = 0}^{i}\varrho_{\nu}^{r, s, p, q}\int_{\Omega}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)dt. \end{equation}$

(3.32)

Through the upcoming section we introduce the shifted MBCSPs operational matrix of the left-sided Caputo's fractional operator in the purpose of investigating the numerical solution of the FDEs in (1.6) with the help of spectral Galerkin technique.

4. Shifted MBCSPs operational matrix

The shifted monic operational matrix of the left-sided Caputo fractional operator will be introduced in this section. Consequently, the operational matrix of this operator with respect to other orthogonal polynomials can be constructed as a special case of the generalized one.

Theorem 4.1. Suppose that ${^C}{\mathcal{D}^{\alpha}_{c+}}$ is the Caputo's fractional differentiation of order $\alpha\in\mathbb{R}$ . The derivative of the shifted MBSCPs vector ${^C}{\mathcal{D}^{\alpha}_{c+}}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = \mathfrak{O}_{\rho, n, \alpha}^{\left(r, s, p, q\right)}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right), \end{equation}$

(4.1)

where $\mathfrak{O}_{\rho, n, \alpha}^{\left(r, s, p, q\right)}$ is the shifted monic BCSPs operational matrix of the Caputo fractional derivative of order $\alpha$ , defined as:

$\begin{equation*} \mathfrak{O}_{\rho, n, \alpha}^{\left(r, s, p, q\right)} = Q^{r, s, p, q}_{\rho, n}\mathcal{H}_{n, \alpha}\left(t\right)\left(Q_{\rho, n}^{r, s, p, q}\right)^{-1}, \end{equation*}$

and $\mathcal{H}_{n, \alpha}\left(t\right)$ is the diagonal matrix of the form:

$\begin{equation} \mathcal{H}_{n, \alpha}\left(t\right) = \begin{pmatrix} 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & \cdots & \frac{\Gamma (\lceil \alpha \rceil +1)}{\Gamma (\lceil \alpha \rceil-\alpha +1)}t^{-\alpha} & \cdots & 0 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & \frac{\Gamma (n+1)}{\Gamma (n-\alpha + 1)}t^{-\alpha} \\ \end{pmatrix}. \end{equation}$

(4.2)

Proof. According to the matrix equation in Lemma 3.2, the derivative of $\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ can be written as:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = Q^{r, s, p, q}_{\rho, n}{^C}{\mathcal{D}^{\alpha}_{c+}}\mathfrak{T}_{n}\left(t\right), \end{equation}$

(4.3)

whereas the derivative of $\mathfrak{T}_{n}\left(t\right)$ , by applying (2.2), is given by:

$\begin{equation} \begin{pmatrix} 0 & \cdots & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & \cdots & \frac{\Gamma (\lceil \alpha \rceil +1)}{\Gamma (\lceil \alpha \rceil-\alpha +1)}t^{\lceil \alpha \rceil-\alpha} & \cdots & 0 \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & \frac{\Gamma (n+1)}{\Gamma (n-\alpha + 1)}t^{n-\alpha} \\ \end{pmatrix}, \end{equation}$

(4.4)

which reads $\mathcal{H}_{n, \alpha}\left(t\right)\mathfrak{T}_{n}\left(t\right)$ . Since the vector $\mathfrak{T}_{n}\left(t\right)$ can be written terms of $\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ as:

$\begin{equation*} \mathfrak{T}_{n}\left(t\right) = \left(Q_{\rho, n}^{r, s, p, q}\right)^{-1}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right). \end{equation*}$

Substituting the last equation in (4.3), we get

$\begin{equation*} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = Q^{r, s, p, q}_{\rho, n}\mathcal{H}_{n, \alpha}\left(t\right)\left(Q_{\rho, n}^{r, s, p, q}\right)^{-1}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right). \end{equation*}$

The previous equation completes the proof. □

The operational matrix with respect to vast families of orthogonal polynomials related to the shifted monic BCSPs can be expressed in terms of the generalized operational matrix given in Theorem 4.1. The following corollaries provide the methodology taken to derive such operational matrices.

Corollary 4.1. Let $\mathcal{P}_{n}^{\ast}\left(t\right) = \left[P_{0}^{\ast}\left(t\right), P_{1}^{\ast}\left(t\right), \cdots, P_{n}^{\ast}\left(t\right)\right]$ be the shifted Legendre polynomials vector defined over the interval $\left[0, 1\right]$ . The shifted Legendre operational matrix of the Caputo fractional operator of order $\alpha$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathcal{P}_{n}^{\ast}\left(t\right) = \left(\mathfrak{D}_{p, n}Q^{-2, 0, -1, 1}_{1, n}\right)\mathcal{H}_{n, \alpha}\left(t\right)\left(\mathfrak{D}_{p, n}Q_{1, n}^{-2, 0, -1, 1}\right)^{-1}\mathcal{P}_{n}^{\ast}\left(t\right), \end{equation}$

(4.5)

where $\mathfrak{D}_{p, n}$ is the $\left(n+1\right)\times\left(n+1\right)$ diagonal matrix given by:

$\begin{equation*} \begin{pmatrix} \frac{(2\cdot0)!}{(0!)^2 2^0} & 0 & \cdots & 0 \\ 0 & \frac{(2\cdot1)!}{(1!)^2 2^1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{(2\cdot n)!}{(n!)^2 2^n} \\ \end{pmatrix}. \end{equation*}$

Proof. Since the shifted Legendre polynomials is related to $\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)$ as in (3.19), by applying Remark 3.1, the matrix representation of $\mathcal{P}_{n}^{\ast}\left(t\right)$ would be:

$\begin{equation} \mathcal{P}_{n}^{\ast}\left(t\right) = diag\left(\frac{\left(2i\right)!}{\left(i!\right)2^{i}}, n\right)Q^{-2, 0, -1, 1}_{1, n}\mathfrak{T}_{n}\left(t\right). \end{equation}$

(4.6)

Acting on the vector $\mathfrak{T}_{n}\left(t\right)$ as in (4.4), we get the desired result. □

Corollary 4.2. Let $\mathcal{U}_{n}^{\lambda^\ast}\left(t\right) = \left[C_{0}^{\lambda^\ast}\left(t\right), C_{1}^{\lambda^\ast}\left(t\right), \cdots, C_{n}^{\lambda^\ast}\left(t\right)\right]$ be the shifted Ultraspherical polynomials vector defined over the interval $\left[0, 1\right]$ . The shifted Ultraspherical operational matrix of the Caputo fractional operator of order $\alpha$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathcal{U}_{n}^{\lambda^\ast}\left(t\right) = \left(\mathfrak{D}_{u, n}Q^{-2\lambda-1, 0, -1, 1}_{1, n}\right)\mathcal{H}_{n, \alpha}\left(t\right)\left(\mathfrak{D}_{u, n}Q_{1, n}^{-2\lambda-1, 0, -1, 1}\right)^{-1}\mathcal{U}_{n}^{\lambda^\ast}\left(t\right), \end{equation}$

(4.7)

where $\mathfrak{D}_{u, n}$ is the $\left(n+1\right)\times\left(n+1\right)$ diagonal matrix given by:

$\begin{equation*} \begin{pmatrix} \frac{2^{0}\left(\lambda\right)_{0}}{(0!)} & 0 & \cdots & 0 \\ 0 & \frac{2^{1}\left(\lambda\right)_{1}}{(1!)} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \frac{2^{n}\left(\lambda\right)_{n}}{(n!)} \\ \end{pmatrix}. \end{equation*}$

Corollary 4.3. Let $\mathcal{T}_{n}^{\ast}\left(t\right) = \left[T_{0}^{\ast}\left(t\right), T_{1}^{\ast}\left(t\right), \cdots, T_{n}^{\ast}\left(t\right)\right]$ be the first-kind shifted Chebyshev polynomials vector defined over the interval $\left[0, 1\right]$ . The shifted Chebyshev of the first-kind operational matrix of the Caputo fractional operator of order $\alpha$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\mathcal{T}_{n}^{\ast}\left(t\right) = \left(\mathfrak{D}_{t, n}Q^{-1, 0, -1, 1}_{1, n}\right)\mathcal{H}_{n, \alpha}\left(t\right)\left(\mathfrak{D}_{t, n}Q_{1, n}^{-1, 0, -1, 1}\right)^{-1}\mathcal{T}_{n}^{\ast}\left(t\right), \end{equation}$

(4.8)

where $\mathfrak{D}_{t, n}$ is the $\left(n+1\right)\times\left(n+1\right)$ diagonal matrix given by:

$\begin{equation*} \begin{pmatrix} 2^{0-1} & 0 & \cdots & 0 \\ 0 & 2^{1-1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 2^{n-1} \\ \end{pmatrix}. \end{equation*}$

Corollary 4.4. Let $\bar{\mathcal{X}}_{n}^{\ast}\left(t\right) = \left[\bar{X}_{0}^{\ast}\left(t\right), \bar{X}_{1}^{\ast}\left(t\right), \cdots, \bar{X}_{n}^{\ast}\left(t\right)\right]$ be the shifted moinc Chebyshev polynomials vector of the fifth-kind defined over the interval $\left[0, 1\right]$ . The shifted moinc Chebyshev of the fifth-kind operational matrix of the Caputo fractional operator of order $\alpha$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\bar{\mathcal{X}}_{n}^{\ast}\left(t\right) = Q^{-3, 2, -1, 1}_{1, n}\mathcal{H}_{n, \alpha}\left(t\right)\left(Q_{1, n}^{-3, 2, -1, 1}\right)^{-1}\bar{\mathcal{X}}_{n}^{\ast}\left(t\right). \end{equation}$

(4.9)

Corollary 4.5. Let $\bar{\mathcal{Y}}_{n}^{\ast}\left(t\right) = \left[\bar{Y}_{0}^{\ast}\left(t\right), \bar{Y}_{1}^{\ast}\left(t\right), \cdots, \bar{Y}_{n}^{\ast}\left(t\right)\right]$ be the shifted moinc Chebyshev polynomials vector of the sixth-kind defined over the interval $\left[0, 1\right]$ . The shifted moinc Chebyshev of the sixth-kind operational matrix of the Caputo fractional operator of order $\alpha$ is given by:

$\begin{equation} {^C}{\mathcal{D}^{\alpha}_{c+}}\bar{\mathcal{Y}}_{n}^{\ast}\left(t\right) = Q^{-5, 2, -1, 1}_{1, n}\mathcal{H}_{n, \alpha}\left(t\right)\left(Q_{1, n}^{-5, 2, -1, 1}\right)^{-1}\bar{\mathcal{Y}}_{n}^{\ast}\left(t\right). \end{equation}$

(4.10)

5. Discussion of error estimate

Here, the error estimate will be examined by the help of the generalized Taylor's formula.

Theorem 5.1. Assume that ${^C}{}{\mathcal{D}}^{i\alpha}_{0+}g\left(t\right)\in \mathcal{C}\left(\Omega\right)$ for $i = 0, 1, \cdots, n$ , and let $u\left(t\right)\in \mathcal{C}^{n}\left[0, \rho\right]$ . Consider $\mathfrak{J}_{n} =$ Span $\left\{\mathcal{G}_{0, \rho}^{\left(r, s, p, q\right)}\left(t\right), \mathcal{G}_{1, \rho}^{\left(r, s, p, q\right)}\left(t\right), \cdots, \mathcal{G}_{n, \rho}^{\left(r, s, p, q\right)}\left(t\right)\right\}$ . If $u_{n}\left(t\right) = \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)$ is the best approximation of the function $u\left(t\right)$ out of $\mathfrak{J}_{n}$ , then the error bound is estimated by:

$\begin{equation} \left\lVert u\left(t\right)-u_{n}\left(t\right)\right\rVert_{{L^{2}_{\mathcal{W}}\left(\Omega\right)}}\leq\frac{M}{\Gamma\left(n\alpha+1\right)}\sqrt{I^{n, \rho, \alpha}_{r, s, p, q}}, \end{equation}$

(5.1)

where $I^{n, \rho, \alpha}_{r, s, p, q}$ is defined as:

$\begin{equation} I^{n, \rho, \alpha}_{r, s, p, q} = \int_{\Omega}t^{2n\alpha}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)dt, \end{equation}$

(5.2)

which converges for all $t\in\Omega$ , and $M = \underset{t\in\Omega}{\sup}\left|{^C}{}{\mathcal{D}}^{n\alpha}_{0+}u\left(t\right)\right|$ .

Proof. By expanding $u\left(t\right)\in L^{2}\left(\Omega\right)$ in terms of the generalized Taylor's formula as given in Definition 2.3

$\begin{equation*} u\left(t\right) = \sum\limits_{i = 0}^{n-1}\frac{t^{i\alpha}}{\Gamma\left(i\alpha+1\right)}{^C}{}{\mathcal{D}}^{\alpha}_{0+}u\left(0^{+}\right)+\frac{t^{n\alpha}}{\Gamma\left(n\alpha+1\right)}{^C}{}{\mathcal{D}}^{\alpha}_{0+}u\left(\xi\right), \end{equation*}$

for $\xi\in\left[0, t\right]$ . Consider the function $\varphi\left(t\right)\in \mathfrak{J}_{n}$ defined as:

$\begin{equation*} \varphi\left(t\right) = \sum\limits_{i = 0}^{n-1}\frac{t^{i\alpha}}{\Gamma\left(i\alpha+1\right)}{^C}{\mathcal{D}^{\alpha}_{0+}}u\left(0^{+}\right), \end{equation*}$

then

$\begin{equation*} \left|u\left(t\right)-\varphi\left(t\right)\right|\leq\frac{t^{n\alpha}}{\Gamma\left(n\alpha+1\right)}\left|{^C}{}{\mathcal{D}}^{\alpha}_{0+}u\left(\xi\right)\right| = \frac{Mt^{n\alpha}}{\Gamma\left(n\alpha+1\right)}. \end{equation*}$

Having that $u_{n}\left(t\right)$ is the best approximation of $u\left(t\right)$ in $\mathfrak{J}_{n}$ , then

$\begin{align*} \left\lVert u\left(t\right)-u_{n}\left(t\right)\right\rVert^{2}_{{L^{2}_{\mathcal{W}}\left(\Omega\right)}}&\leq \left\lVert u\left(t\right)-\varphi\left(t\right)\right\rVert^{2}_{{L^{2}_{\mathcal{W}}\left(\Omega\right)}} = \int_{\Omega}\left|u\left(t\right)-\varphi\left(t\right)\right|^{2}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)dt\\&\leq\frac{M^{2}}{\Gamma\left(n\alpha+1\right)^{2}}\int_{\Omega}t^{2n\alpha}\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)dt, \end{align*}$

which proves the theorem. □

The next theorem gives the residual error associated with the FDEs considered in Eq (1.6).

Theorem 5.2. Suppose that $f:\left[0, T\right]\times\Omega\rightarrow \mathbb{R}$ satisfies a Lipschitz condition in $u$ with Lipschitz constant $L$ , that is

$\begin{equation} \left\|f\left(t, u\right)-f\left(t, \tilde{u}\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\leq L\left\|u\left(t\right)-\tilde{u}\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}. \end{equation}$

(5.3)

Then, the residual error is bounded by:

$\begin{equation} \left\|R\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\leq\left(\sum\limits_{k = 0}^{N}\vartheta_{k}\xi_{k}+L\right)\left\|u_{n}\left(t\right)-u\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}, \end{equation}$

(5.4)

where $\left\|u_{n}\left(t\right)-u\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}$ is the truncation error in Theorem 5.1.

Proof. Consider the approximate spectral solution $u_{n}\left(t\right)$ of (1.6). Since the residual error is defined as:

$\begin{equation*} R\left(t\right) = \sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}u_{n}\left(t\right)-f\left(t, u_{n}\left(t\right)\right), \end{equation*}$

then

$\begin{align*} \left\|R\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}& = \left\|\sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}u_{n}\left(t\right)-f\left(t, u_{n}\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\\& = \Bigg\|\sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}u_{n}\left(t\right)-f\left(t, u_{n}\right)\\& \qquad\qquad\qquad\qquad\qquad-\left(\sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}u\left(t\right)-f\left(t, u\right)\right)\Bigg\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\\& = \left\|\sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}\left(u_{n}\left(t\right)-u\left(t\right)\right)-\left(f\left(t, u_{n}\right)-f\left(t, u\right)\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\\&\leq\left\|\sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}\left(u_{n}\left(t\right)-u\left(t\right)\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}+\left\|f\left(t, u_{n}\right)-f\left(t, u\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\\&\leq\sum\limits_{k = 0}^{N}\vartheta_{k}\left\|{^C}{\mathcal{D}^{\alpha_{k}}_{0+}}\left(u_{n}\left(t\right)-u\left(t\right)\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}+ L\left\|u_{n}\left(t\right)-u\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}, \end{align*}$

where

$\begin{equation} \vartheta_{k} = \underset{t\in\Omega}{\sup}\left\|\gamma_{k}\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}. \end{equation}$

(5.5)

Now, from the boundedness of ${^C}{\mathcal{D}^{\alpha_{k}}_{0+}}$ (see ^[33]), for each $k$ in $\left\{0, 1, \cdots, N\right\}$ , we have

$\begin{equation} \left\|{^C}{\mathcal{D}^{\alpha_{k}}_{0+}}\left(u_{n}\left(t\right)-u\left(t\right)\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\leq\xi_{k}\left\|u_{n}\left(t\right)-u\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}. \end{equation}$

(5.6)

Hence,

$\begin{equation*} \left\|R\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}\leq\left(\sum\limits_{k = 0}^{N}\vartheta_{k}\xi_{k}+L\right)\left\|u_{n}\left(t\right)-u\left(t\right)\right\|_{L^{2}_{\mathcal{W}}\left(\Omega\right)}, \end{equation*}$

which tends to zero as shown in Theorem 5.1 when $n$ is sufficiently large. □

6. Numerical scheme

In this section, we introduce the method of solution of the multi-term FDE in (1.6). Unlike other traditional numerical schemes, the proposed scheme based on the shifted MBCSPs can be used to study the numerical solution with respect to numerous orthogonal polynomials without performing further calculations.

6.1. Method of solution

Consider the multi-term FDE in (1.6), by approximating the function $u\left(t\right)\in L^{2}_{\mathcal{W}}\left(\Omega\right)$ as in (3.30), the FDE reads

$\begin{equation*} \sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right){^C}{\mathcal{D}^{\alpha_{k}}_{0+}}\Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = f\left(t, \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)\right), \end{equation*}$

now, applying Theorem 4.1 implies

$\begin{equation*} \sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right)\Lambda^{T}\mathfrak{O}_{\rho, n, \alpha_{k}}^{\left(r, s, p, q\right)}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right) = f\left(t, \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)\right), \end{equation*}$

hence, the residual of the problem would be

$\begin{equation} R\left(t\right) = \sum\limits_{k = 0}^{N}\gamma_{k}\left(t\right)\Lambda^{T}\mathfrak{O}_{\rho, n, \alpha_{k}}^{\left(r, s, p, q\right)}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)-f\left(t, \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(t\right)\right). \end{equation}$

(6.1)

The application of spectral Galerkin method yields

$\begin{equation} \langle\frac{\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)}{\mathcal{W}^{\left(r, s, p, q\right)}_{\rho}\left(t\right)}, R\left(t\right)\rangle = \int_{\Omega}\mathcal{G}_{i, \rho}^{\left(r, s, p, q\right)}\left(t\right)R\left(t\right)dt = 0, \end{equation}$

(6.2)

for $i = 0, 1, \cdots, n-2$ . the above equation generates $n-1$ dimensional system of equation in the unknowns $c_{i}$ . Making use of the boundary conditions, we get the following additional equation

$\begin{equation} \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(0\right) = u_{0}, \quad and \quad \Lambda^{T}\mathfrak{G}_{\rho, n}^{\left(r, s, p, q\right)}\left(T\right) = u_{1}. \end{equation}$

(6.3)

By solving the obtained system using a suitable solver, the approximate solution can be expressed in terms of the shifted MBCSPs.

6.2. Numerical experiments

In the purpose of illustrating the accuracy and applicability of the proposed generalized approach, a numerical treatment for a different FDEs supported with comparison with other techniques will be present in this section.

The error estimate between the exact solution $u\left(t\right)$ and the numerical solution $\tilde{u}\left(t\right)$ , will be measured by calculating the absolute value $e_{i} = \left|u\left(t_{i}\right)-\tilde{u}\left(t_{i}\right)\right|$ for $i = 1, 2, \cdots, n$ .

Example 6.1. Consider the following Bagley-Torvik initial-value problem (IVP) ^[18,34]:

$\begin{equation} {^C}{}{}\mathcal{D}^{2}_{0+}u\left(t\right)+{^C}{}{}\mathcal{D}^{\alpha}_{0+}u\left(t\right)+u\left(t\right) = f\left(t\right), \quad t\in\left[0, \rho = 1\right], \end{equation}$

(6.4)

subject to $u\left(0\right) = 1$ which has the exact solution $u\left(t\right) = 1+t$ . For $\alpha = \frac{1}{2}$ , and $n = 2$ , writing $u\left(t\right)$ as in (3.30)

$\begin{equation} u_{2}\left(t\right) = c_{0}\mathcal{G}_{0}^{\left(r, s, p, q\right)}\left(t\right)+c_{1}\mathcal{G}_{1}^{\left(r, s, p, q\right)}\left(t\right)+c_{2}\mathcal{G}_{2}^{\left(r, s, p, q\right)}\left(t\right), \end{equation}$

(6.5)

by applying (6.1), the residual of given equation reads

$\begin{equation} R\left(t\right) = \Lambda^{T}\left(\mathfrak{O}_{1, 2, 2}^{\left(r, s, p, q\right)}\mathfrak{G}_{1, 2}^{\left(r, s, p, q\right)}\left(t\right)+\mathfrak{O}_{1, 2, \frac{1}{2}}^{\left(r, s, p, q\right)}\mathfrak{G}_{1, 2}^{\left(r, s, p, q\right)}\left(t\right)+\mathfrak{G}_{1, 2}^{\left(r, s, p, q\right)}\left(t\right)\right)-f\left(t\right), \end{equation}$

(6.6)

where $f\left(t\right) = t+\frac{2 \sqrt{t}}{\sqrt{\pi }}+1$ . Choosing $r = -5, s = 4, p = -q = -1$ , the proposed scheme with the weights $\mathcal{G}_{k}^{\left(-5, 4, -1, 1\right)}\left(t\right)$ , $k = 0, 1, 2$ , provide the coefficients

$\begin{equation*} u_{2}\left(t\right) = \left(\frac{3}{2}\right)\mathcal{G}_{0}^{\left(-5, 4, -1, 1\right)}\left(t\right)+\left(\frac{1}{2}\right)\mathcal{G}_{1}^{\left(-5, 4, -1, 1\right)}\left(t\right)+\left(0\right)\mathcal{G}_{2}^{\left(-5, 4, -1, 1\right)}\left(t\right) = 1+t, \end{equation*}$

which is the exact solution. For small values of $\alpha = 0.1, 0.05,$ and $0.01$ with $n = 2$ the following results are obtained in Table 1.

Table 1. The numerical results of Example 6.1 for different values of

$\alpha$ at

$n = 2$ .

$t$	Exact	Absolute error	Absolute error	Absolute error
	solution	$\left(\alpha=0.1\right)$	$\left(\alpha=0.05\right)$	$\left(\alpha=0.01\right)$
0.0	1.00000	$2.66454\times10^{-15}$	$9.10383\times10^{-15}$	$1.11022\times10^{-15}$
0.2	1.20000	$3.10862\times10^{-15}$	$9.99201\times10^{-15}$	$1.11022\times10^{-15}$
0.4	1.40000	$3.33067\times10^{-15}$	$1.04361\times10^{-14}$	$1.33227\times10^{-15}$
0.6	1.60000	$3.33067\times10^{-15}$	$1.04361\times10^{-14}$	$1.11022\times10^{-15}$
0.8	1.80000	$3.10862\times10^{-15}$	$9.99201\times10^{-15}$	$1.11022\times10^{-15}$
1.0	2.00000	$3.10862\times10^{-15}$	$9.10383\times10^{-15}$	$8.88178\times10^{-16}$

| Show Table

DownLoad: CSV

Example 6.2. Consider the linear fractional IVP of the form ^[35]:

$\begin{equation} {^C}{}{}\mathcal{D}^{2}_{0+}u\left(t\right)-2{^C}{}{}\mathcal{D}^{1}_{0+}u\left(t\right)+{^C}{}{}\mathcal{D}^{\frac{1}{2}}_{0+}u\left(t\right) = f\left(t, u\left(t\right)\right), \quad t\in\left[0, \rho = 1\right], \end{equation}$

(6.7)

where $f\left(t, u\left(t\right)\right)$ is defined as:

$\begin{equation*} f\left(t, u\left(t\right)\right) = t^{7}-14t^{6}+42t^{5}-t^{2}+\frac{8\left(256t^{5}-143\right)t^{\frac{3}{2}}}{429\sqrt{\pi}}+4\left(t-\frac{1}{2}\right)-u\left(t\right), \end{equation*}$

having the initial conditions $u\left(0\right) = u^{\prime}\left(0\right) = 0$ with the exact solution $u\left(t\right) = t^{2}\left(t^{5}-1\right)$ . Similarly as in Example 6.1, the residual of the problem can be derived. Assuming that

$\begin{equation*} u\left(t\right)\simeq u_{7}\left(t\right) = \sum\limits_{i = 0}^{7}c_{i}\mathcal{G}_{i}^{\left(-4, 5, -1, 1\right)}\left(t\right), \end{equation*}$

by constructing the algebraic system with the aid of Galerkin technique, the approximate solution $u_{7}\left(t\right)$ is found to be

$\begin{equation*} u_{7}\left(t\right) = t^{2}\left(t^{5}-1\right), \end{equation*}$

which is the same as the exact solution. Table 2 shows the comparison between the maximal absolute error of the present method and other polynomial based techniques such as Tau-Shifted Chebyshev (TC), Legendre-Gauss quadrature (LGQ), triangular functions (TF), and block pulse functions (BPF).

Table 2. The maximal absolute error of Example 6.2.

Present method $(n=8)$	TSC $(n=8)$	LGQ $(n=8)$	TF $(h=0.01)$	BPF $(h=0.01)$
0.00000	$3.07566\times10^{-15}$	$1.22696\times10^{-15}$	$3.96626\times10^{-4}$	$2.61473\times10^{-2}$

| Show Table

DownLoad: CSV

Example 6.3. Consider the Bagley-Torvik boundary value problem (BVP) ^[34]:

$\begin{equation} {^C}{}{}\mathcal{D}^{\alpha}_{c+}u\left(t\right) = f\left(t, u\left(t\right)\right), \quad t\in\left[0, \rho = 1\right], \end{equation}$

(6.8)

for $\alpha = \frac{3}{2}$ , the function $f\left(t, u\left(t\right)\right)$ takes the form:

$\begin{equation*} f\left(t, u\left(t\right)\right) = t^{5}-t^{4}+\frac{128t^{\frac{7}{2}}}{7\sqrt{\pi}}-\frac{64t^{\frac{5}{2}}}{5\sqrt{\pi}}-u\left(t\right), \end{equation*}$

governed by the boundary conditions $u\left(0\right) = u\left(1\right) = 0$ , with the exact solution $u\left(t\right) = t^{4}\left(t-1\right)$ . By applying the proposed method for $r = -5, s = 4$ and $p = -q = -1$ , the exact solution can be found with $n = 5$ , i.e.,

$\begin{equation*} u_{5}\left(t\right) = \sum\limits_{i = 0}^{5}c_{i}\mathcal{G}_{i}^{\left(-5, 4, -1, 1\right)}\left(t\right) = t^{4}\left(t-1\right). \end{equation*}$

The graph of the exact and numerical solutions is given in Figure 2(a), in addition, the comparison between the absolute error of the present method and the Gegenbauer wavelet method ^[34] is given graphically in Figure 2(b). Therefore, the numerical results is compared with the Gegenbauer wavelet method ^[34] in . Also, provides the absolute error for small values of $\alpha$ .

Figure 2. (a) The graph of exact and numerical solutions, (b) comparison between the absolute error of the present method and the method in ^[34].

DownLoad: Full-Size Img PowerPoint

Table 3. The numerical results of Example 6.3 for

$n = 5$ .

$t$	Exact	Present	Absolute error	Absolute error
	solution	method	$(n=5)$	in ^[34]
0.0	-0.00000	-0.00000	0.00000	$0.00000\times10^{-6}$
0.1	-0.00009	-0.00009	0.00000	$2.00000\times10^{-6}$
0.2	-0.00128	-0.00128	0.00000	$3.00000\times10^{-6}$
0.3	-0.00567	-0.00567	0.00000	$1.00000\times10^{-6}$
0.4	-0.01536	-0.01536	0.00000	$0.00000\times10^{-6}$
0.5	-0.03125	-0.03125	0.00000	$2.00000\times10^{-6}$
0.6	-0.05184	-0.05184	0.00000	$1.00000\times10^{-6}$
0.7	-0.07203	-0.07203	0.00000	$3.00000\times10^{-6}$
0.8	-0.08192	-0.08192	0.00000	$2.00000\times10^{-6}$
0.9	-0.06561	-0.06561	0.00000	$0.00000\times10^{-6}$
1.0	-0.00000	-0.00000	0.00000	$0.00000\times10^{-6}$

| Show Table

DownLoad: CSV

Table 4. The numerical Absolute error of Example 6.3 for different values of

$\alpha$ at

$n = 2$ .

$t$	Exact	Absolute error	Absolute error	Absolute error
	solution	$\left(\alpha=0.25\right)$	$\left(\alpha=0.1\right)$	$\left(\alpha=0.05\right)$
0.0	-0.00000	$1.04083\times10^{-17}$	$3.46945\times10^{-18}$	$1.73472\times10^{-18}$
0.2	-0.00218	$7.80626\times10^{-18}$	$1.30104\times10^{-18}$	$2.58040\times10^{-17}$
0.4	-0.01536	$8.67362\times10^{-18}$	$2.77556\times10^{-17}$	$3.12250\times10^{-17}$
0.6	-0.05184	$2.08167\times10^{-17}$	$8.32667\times10^{-17}$	$9.71445\times10^{-17}$
0.8	-0.08192	$5.55112\times10^{-17}$	$6.93889\times10^{-17}$	$8.32667\times10^{-17}$
1.0	-0.00000	$1.56125\times10^{-16}$	$2.25514\times10^{-16}$	$2.23779\times10^{-16}$

| Show Table

DownLoad: CSV

One of the main advantages of the proposed scheme appears from the ability of providing the approximate solution in terms of different families of polynomials. So, compares the absolute error between the shifted Legendre (SL), first-kind shifted Chebyshev (FSC), second-kind shifted Chebyshev (SSC) and the shown configuration of the MBCSPs for Example 6.3 with $\alpha = \frac{1}{4}$ . provides the numerical results obtained with different $r, s, p$ and $q$ configurations representing other families of polynomials for $n = 5$ .

Figure 3. Comparison of the absolute error between, SL, FSC, SSC, and MBCSP.

DownLoad: Full-Size Img PowerPoint

Table 5. The numerical results of Example 6.3 with

$\alpha = 0.25$ and

$n = 5$ for different families of polynomials.

$t$	Exact	SL	FSC	SSC
	solution	solution	solution	solution
0.0	-0.00000	$1.73472 ^{-17}$	$3.46945 ^{-18}$	$1.73472 ^{-17}$
0.2	-0.00128	$9.97466 ^{-18}$	$7.58942 ^{-18}$	$1.51788 ^{-18}$
0.4	-0.01536	$1.04083 ^{-17}$	$2.08167 ^{-17}$	$3.46945 ^{-17}$
0.6	-0.05184	$4.16334 ^{-17}$	$6.93889 ^{-17}$	$1.38778 ^{-16}$
0.8	-0.08192	$5.55112 ^{-17}$	$8.32667 ^{-17}$	$2.63678 ^{-16}$
1.0	-0.00000	$7.97973 ^{-17}$	$1.90820 ^{-16}$	$5.44703 ^{-16}$

| Show Table

DownLoad: CSV

Example 6.4. Consider the linear boundary value problem ^[36]:

$\begin{equation} \left(10+t\right)^{2}{^C}{}{}\mathcal{D}^{\frac{5}{2}}_{0+}u\left(t\right)+\frac{5}{2}\left(10+t\right){^C}{}{}\mathcal{D}^{\frac{3}{2}}_{0+}u\left(t\right)+\frac{1}{2}{^C}{}{}\mathcal{D}^{\frac{1}{2}}_{0+}u\left(t\right) = \frac{t^{\frac{3}{2}}}{100\sqrt{\pi}}, \quad t\in\left[0, \rho = 1\right], \end{equation}$

(6.9)

governed by the boundary conditions $u\left(0\right) = \ln\left(10\right), u\left(1\right) = \ln\left(11\right)$ , with the exact solution $u\left(t\right) = \ln\left(10+t\right)$ . By applying the proposed method for $r = -7, s = 4$ and $p = -q = -1$ . , obtained The graph of exact and numerical solutions at $n = 12$ , whereas, gives the absolute error for different values of $n$ , which reflects the convergence of the proposed method by increasing the number of basis. Also, gives The numerical results of the present example at $n = 8, 10$ and $12$ .

Figure 4. (a) The graph of exact and numerical solutions, (b) the absolute error of the present method for

$n = 8, 10$ and

$12$ .

DownLoad: Full-Size Img PowerPoint

Table 6. The numerical results of Example 6.4 for

$n = 8, 10$ and

$12$ .

$t$	Exact	Present	Absolute error	Absolute error	Absolute error
	solution	method	$(n=8)$	$(n=10)$	$(n=12)$
0.0	2.30259	2.30259	$8.88178\times10^{-16}$	$8.88178\times10^{-16}$	$4.44089\times10^{-16}$
0.1	2.31254	2.31254	$2.22045\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.2	2.32239	2.32239	$1.33227\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.3	2.33214	2.33214	$3.10862\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.4	2.34181	2.34181	$2.66454\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.5	2.35138	2.35138	$2.66454\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.6	2.36085	2.36085	$5.32907\times10^{-15}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.7	2.37024	2.37024	$8.88178\times10^{-16}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
0.8	2.37955	2.37955	$3.10862\times10^{-15}$	$8.88178\times10^{-16}$	$4.44089\times10^{-16}$
0.9	2.38876	2.38876	$4.44089\times10^{-16}$	$4.44089\times10^{-16}$	$0.00000\times10^{-16}$
1.0	2.39789	2.39789	$0.00000\times10^{-16}$	$1.33227\times10^{-15}$	$0.00000\times10^{-16}$

| Show Table

DownLoad: CSV

In case that we want to provide the approximate solution in terms of different families of polynomials. shows the absolute error of the approximate solution in terms of the SL, SFC, SSC, and the monic shifted sixth-kind Chebyshev polynomials (MSSC) obtained as a special case of $\mathcal{G}^{\left(r, s, p, q\right)}_{i, \rho}\left(t\right)$ . It is clear that the Chebyshev family has similar performance in this problem. Finally, gives the absolute error of different polynomials (SL solution, SFC solution, SSC solution, and MSSC solution) at $n = 8$ .

Figure 5. Graph of the absolute error in Example 6.4 of different families of polynomials for

$n = 8$ , where, (a) SL solution, (b) SFC solution, (c) SSC solution, (d) MSSC solution.

DownLoad: Full-Size Img PowerPoint

Table 7. The absolute error in Example 6.4 of different polynomials for

$n = 8$ .

$t$	SL Absolute	SFC Absolute	SSC Absolute	MSSC Absolute
	Error	Error	Error	Error
0.0	$0.00000 ^{-16}$	$8.88178 ^{-16}$	$4.44089 ^{-16}$	$4.44089 ^{-16}$
0.2	$4.44089 ^{-16}$	$1.33227 ^{-16}$	$0.00000 ^{-16}$	$0.00000 ^{-16}$
0.4	$3.55271 ^{-15}$	$2.66454 ^{-15}$	$3.99680 ^{-15}$	$3.99680 ^{-15}$
0.6	$4.44089 ^{-15}$	$5.32907 ^{-15}$	$3.99680 ^{-15}$	$4.44089 ^{-15}$
0.8	$3.99680 ^{-15}$	$3.10862 ^{-15}$	$3.99680 ^{-15}$	$3.99680 ^{-15}$
1.0	$4.44089 ^{-16}$	$8.88178 ^{-16}$	$4.44089 ^{-16}$	$0.00000 ^{-16}$

| Show Table

DownLoad: CSV

7. Conclusions

Through the discussion in this paper, the authors have introduced a novel generalized shifted symmetric orthogonal basis which was used later for obtaining a general operational matrix for the left-sided Caputo fractional derivative. This class of polynomials presumably used for the first time as trial functions with spectral methods in a numerical scheme concerning the solution of fractional differential equations (FDEs). The obtained operational matrix with four parameters can be used for expressing the operational matrix with respect to other families without needing additional calculations. For the sake of demonstration, the proposed basis were used for solving multiple FDEs showing its accuracy and applicability along with comparing it with other different techniques. The advantage of these orthogonal polynomials is that it include numerous choices for the basis of solution with different domains of definitions as well as generating a new orthogonal basis. This makes the task of choosing the most appropriate basis for the problem under study much easier by comparing the error for different four parameters configurations.

Acknowledgments

The authors are thankful to the Taif University, Taif, Saudi Arabia. The research was funded by Taif University TURSP-2020/160.

Conflict of interest

The authors declare that they have no competing interests.

References

[1]	N. Sene, Analytical solutions and numerical schemes of certain generalized fractional diffusion models, Eur. Phys. J. Plus, 134 (2019), 199. https://doi.org/10.1140/epjp/i2019-12531-4 doi: 10.1140/epjp/i2019-12531-4
[2]	A. Atangana, S. Qureshi, Modeling attractors of chaotic dynamical systems with fractal-fractional operators, Chaos, Solitons Fract., 123 (2019), 320–337. https://doi.org/10.1016/j.chaos.2019.04.020 doi: 10.1016/j.chaos.2019.04.020
[3]	R. P. Yadav, R. Verma, A numerical simulation of fractional order mathematical modeling of COVID-19 disease in case of Wuhan China, Chaos, Solitons Fract., 140 (2020), 110124. https://doi.org/10.1016/j.chaos.2020.110124 doi: 10.1016/j.chaos.2020.110124
[4]	C. M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing sub-diffusion, J. Comput. Phys., 227 (2007), 886–897. https://doi.org/10.1016/j.jcp.2007.05.012 doi: 10.1016/j.jcp.2007.05.012
[5]	C. M. Chen, F. Liu, K. Burrage, Finite difference methods and a Fourier analysis for the fractional reaction-subdiffusion equation, Appl. Math. Comput., 198 (2008), 754–769. https://doi.org/10.1016/j.amc.2007.09.020 doi: 10.1016/j.amc.2007.09.020
[6]	H. Gul, S. Ali, K. Shah, S. Muhammad, T. Sitthiwirattham, S. Chasreechai, Application of asymptotic homotopy perturbation method to fractional order partial differential equation, Symmetry, 13 (2021), 2215. https://doi.org/10.3390/sym13112215 doi: 10.3390/sym13112215
[7]	M. S. Arshad, D. Baleanu, M. B. Riaz, M. Abbas, A novel 2-stage fractional Runge–kutta method for a time-fractional logistic growth model, Discrete Dyn. Nat. Soc., 2020 (2020), 1–8. https://doi.org/10.1155/2020/1020472 doi: 10.1155/2020/1020472
[8]	O. González-Gaxiola, A. Biswas, W-shaped optical solitons of Chen–Lee–Liu equation by Laplace–Adomian decomposition method, Opt. Quant. Electron., 50 (2018), 1–11. https://doi.org/10.1007/s11082-018-1583-0 doi: 10.1007/s11082-018-1583-0
[9]	M. Obeid, M. A. Abd El Salam, J. A. Younis, Operational matrix-based technique treating mixed type fractional differential equations via shifted fifth-kind Chebyshev polynomials, Appl. Math. Sci. Eng., 31 (2023), 2187388. https://doi.org/10.1080/27690911.2023.2187388 doi: 10.1080/27690911.2023.2187388
[10]	S. Ibrahim, A. Isah, Solving system of fractional order differential equations using Legendre operational matrix of derivatives, Eurasian J. Sci. Eng., 7 (2021), 25–37. https://doi.org/10.23918/eajse.v7i1p25 doi: 10.23918/eajse.v7i1p25
[11]	A. Zamiri, A. Borhanifar, A. Ghannadiasl, Laguerre collocation method for solving Lane–Emden type equations, Comput. Methods Differ. Equ., 9 (2021), 1176–1197. https://dx.doi.org/10.22034/cmde.2020.35895.1621 doi: 10.22034/cmde.2020.35895.1621
[12]	A. Yari, Numerical solution for fractional optimal control problems by Hermite polynomials, J. Vib. Control, 27 (2021), 698–716. https://doi.org/10.1177/1077546320933129 doi: 10.1177/1077546320933129
[13]	M. H. Derakhshan, Numerical solution of a coupled system of fractional order integro differential equations by an efficient numerical method based on the second kind Chebyshev polynomials, Math. Anal. Contemp. Appl., 3 (2021), 25–40. https://doi.org/10.30495/maca.2021.1938222.1025 doi: 10.30495/maca.2021.1938222.1025
[14]	S. N. Tural-Polat, A. T. Dincel, Numerical solution method for multi-term variable order fractional differential equations by shifted Chebyshev polynomials of the third kind, Alex. Eng. J., 61 (2022), 5145–5153. https://doi.org/10.1016/j.aej.2021.10.036 doi: 10.1016/j.aej.2021.10.036
[15]	C. Cesarano, S. Pinelas, P. E. Ricci, The third and fourth kind pseudo-Chebyshev polynomials of half-integer degree, Symmetry, 11 (2019), 274. https://doi.org/10.3390/sym11020274 doi: 10.3390/sym11020274
[16]	A. Secer, S. Altun, A new operational matrix of fractional derivatives to solve systems of fractional differential equations via legendre wavelets, Mathematics, 6 (2018), 238. https://doi.org/10.3390/math6110238 doi: 10.3390/math6110238
[17]	C. Baishya, P. Veeresha, Laguerre polynomial-based operational matrix of integration for solving fractional differential equations with non-singular kernel, Proc. R. Soc. A, 477 (2021), 20210438. https://doi.org/10.1098/rspa.2021.0438 doi: 10.1098/rspa.2021.0438
[18]	W. M. Abd-Elhameed, Y. H. Youssri, Fifth-kind orthonormal Chebyshev polynomial solutions for fractional differential equations, Comput. Appl. Math., 37 (2018), 2897–2921. https://doi.org/10.1007/s40314-017-0488-z doi: 10.1007/s40314-017-0488-z
[19]	K. Sadri, H. Aminikhah, Chebyshev polynomials of sixth kind for solving nonlinear fractional PDEs with proportional delay and its convergence analysis, J. Funct. Spaces, 2022 (2022), 1–20. https://doi.org/10.1155/2022/9512048 doi: 10.1155/2022/9512048
[20]	M. Abdelhakem, D. Mahmoud, D. Baleanu, M. El-kady, Shifted ultraspherical pseudo-Galerkin method for approximating the solutions of some types of ordinary fractional problems, Adv. Differ. Equ., 2021 (2021), 1–18. https://doi.org/10.1186/s13662-021-03247-6 doi: 10.1186/s13662-021-03247-6
[21]	F. A. Shah, R. Abass, Solution of fractional oscillator equations using ultraspherical wavelets, Int. J. Geom. Methods Mod. Phys., 16 (2019), 1950075. https://doi.org/10.1142/S0219887819500750 doi: 10.1142/S0219887819500750
[22]	D. Fortunato, N. Hale, A. Townsend, The ultraspherical spectral element method, J. Comput. Phys., 436 (2021), 110087. https://doi.org/10.1016/j.jcp.2020.110087 doi: 10.1016/j.jcp.2020.110087
[23]	A. A. El-Sayed, D. Baleanu, P. Agarwal, A novel Jacobi operational matrix for numerical solution of multi-term variable-order fractional differential equations, J. Taibah Univ. Sci., 14 (2020), 963–974. https://doi.org/10.1080/16583655.2020.1792681 doi: 10.1080/16583655.2020.1792681
[24]	K. Maleknejad, J. Rashidinia, T. Eftekhari, A new and efficient numerical method based on shifted fractional‐order Jacobi operational matrices for solving some classes of two‐dimensional nonlinear fractional integral equations, Numer. Methods Partial Differ. Equ., 37 (2021), 2687–2713. https://doi.org/10.1002/num.22762 doi: 10.1002/num.22762
[25]	J. Rashidinia, T. Eftekhari, K. Maleknejad, Numerical solutions of two-dimensional nonlinear fractional Volterra and Fredholm integral equations using shifted Jacobi operational matrices via collocation method, J. King Saud Univ.-Sci., 33 (2021), 101244. https://doi.org/10.1016/j.jksus.2020.101244 doi: 10.1016/j.jksus.2020.101244
[26]	M. Masjed-Jamei, A generalization of classical symmetric orthogonal functions using a symmetric generalization of Sturm–Liouville problems, Integr. Transf. Spec. Funct., 18 (2007), 871–883. https://doi.org/10.1080/10652460701510949 doi: 10.1080/10652460701510949
[27]	S. G. Samko, A. A. Kilbas, O. I. Marichev, Fractional integrals and derivatives, Vol. 1, Yverdon-les-Bains, Switzerland: Gordon and breach science publishers, Yverdon, 1993.
[28]	Z. M. Odibat, S. Momani, An algorithm for the numerical solution of differential equations of fractional order, J. Appl. Math. Inform., 26 (2008), 15–27.
[29]	Z. M. Odibat, N. T. Shawagfeh, Generalized Taylor's formula, Appl. Math. Comput., 186 (2007), 286–293. https://doi.org/10.1016/j.amc.2006.07.102 doi: 10.1016/j.amc.2006.07.102
[30]	A. El-Ajou, O. A. Arqub, M. Al-Smadi, A general form of the generalized Taylor's formula with some applications, Appl. Math. Comput., 256 (2015), 851–859. https://doi.org/10.1016/j.amc.2015.01.034 doi: 10.1016/j.amc.2015.01.034
[31]	A. E. Choque-Rivero, I. Area, A Favard type theorem for Hurwitz polynomials, Discrete Cont. Dyn. Syst.-Ser. B, 25 (2020), 529–544. https://doi/10.3934/dcdsb.2019252 doi: 10.3934/dcdsb.2019252
[32]	H. M. Srivastava, J. Choi, Zeta and q-Zeta functions and associated series and integrals, El-sevier, 2011.
[33]	A. W. Naylor, G. R. Sell, Linear operator theory in engineering and science, Springer Science and Business Media, 1982.
[34]	H. M. Srivastava, F. A. Shah, R. Abass, An application of the Gegenbauer wavelet method for the numerical solution of the fractional Bagley–Torvik equation, Russ. J. Math. Phys., 26 (2019), 77–93. https://doi.org/10.1134/S1061920819010096 doi: 10.1134/S1061920819010096
[35]	S. K. Damarla, M. Kundu, Numerical solution of multi-order fractional differential equations using generalized triangular function operational matrices, Appl. Math. Comput., 263 (2015), 189–203. https://doi.org/10.1016/j.amc.2015.04.051 doi: 10.1016/j.amc.2015.04.051
[36]	W. M. Abd-Elhameed, Y. H. Youssri, Sixth-kind Chebyshev spectral approach for solving fractional differential equations, Int. J. Nonlinear Sci. Numer. Simul., 20 (2019), 191–203. https://doi.org/10.1515/ijnsns-2018-0118 doi: 10.1515/ijnsns-2018-0118

Reader Comments

Your name:*

Email:*
© 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

AIMS Mathematics

1.8 3.4

Metrics

Article views(1405) PDF downloads(56) Cited by(0)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(5) / Tables(7)

AIMS Mathematics

A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel

Related Papers:

Abstract

1. Introduction

2. Preliminaries

3. Basic symmetric polynomials class

3.1. BCSPs properties

3.2. The shifted class of MBCSPs

4. Shifted MBCSPs operational matrix

5. Discussion of error estimate

6. Numerical scheme

6.1. Method of solution

6.2. Numerical experiments

7. Conclusions

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

AIMS Mathematics

A novel generalized symmetric spectral Galerkin numerical approach for solving fractional differential equations with singular kernel

Related Papers:

Abstract

1. Introduction

2. Preliminaries

3. Basic symmetric polynomials class

3.1. BCSPs properties

3.2. The shifted class of MBCSPs

4. Shifted MBCSPs operational matrix

5. Discussion of error estimate

6. Numerical scheme

6.1. Method of solution

6.2. Numerical experiments

7. Conclusions

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog