Adopted spectral tau approach for the time-fractional diffusion equation via seventh-kind Chebyshev polynomials

This study utilizes a spectral tau method to acquire an accurate numerical solution of the time-fractional diﬀusion equation. The central point of this approach is to use double basis functions in terms of certain Chebyshev polynomials, namely Chebyshev polynomials of the seventh-kind and their shifted ones. Some new formulas concerned with these polynomials are derived in this study. A rigorous error analysis of the proposed double expansion further corroborates our research. This analysis is based on establishing some inequalities regarding the selected basis functions. Several numerical examples validate the precision and eﬀectiveness of the suggested method.


Introduction
Because of their versatility in describing a broad spectrum of physical, chemical, biological, and other systems, fractional differential equations (DEs) have gained much attention in recent decades [1,2].These equations have become indispensable for modeling complicated phenomena with long-range interactions and memory effects [3,4].Recent articles have demonstrated an interest in solving fractional differential models encountered in various applied sciences.For instance, some fractional order epidemic models were examined in [5,6].In [7], the authors developed some solutions for a specific fractional financial model.An additional fractional differential equation that arises in fluid mechanics was examined in [8].
The novel results for approximating solutions are significant since analytic solutions for fractional DEs are typically not explicitly found.As a result, numerous numerical approaches for solving fractional DEs have been of great interest.Researchers used a variety of numerical algorithms to solve various fractional DEs.For example, two finite difference methods are followed in [9] for treating the time-fractional Fokker-Planck equation.The Galerkin approach is followed in [10] to solve a linear hyperbolic first-order partial DEs.The collocation method is applied in [11] to solve a specific fractional model.The homotopy analysis method treats partial fractional DEs in [12].
The time-fractional diffusion equation (TFDE), a basic modification of the classical diffusion equation that considers fractional derivatives in time, is at the core of this mathematical framework [13].The importance of the TFDE extends beyond academic fields, providing essential insights into anomalous diffusion processes found in biological, physical, and engineering systems, including viscoelastic materials, subdiffusion phenomena, and porous media flow [14][15][16].
One class of older polynomials defined by Pafnuty Chebyshev is the set of orthogonal polynomials known as Chebyshev polynomials (CPs), see [17].Most CPs utilized in various applied science domains are of the first and second kinds.These polynomials are helpful in many fields.Their significance in electrohydrodynamics, optimum control, and signal processing has been noted [18][19][20].Many authors used the first-and second-kinds of CPs.For example, the authors of [21] handled some real-life problems via first-derivative CPs of the first kind.Lane-Emden equations are solved via CPs of the first kind in [22].The authors in [23] used the second-kind CPs for higher-order DEs.The kinds of CPs are not restricted to the first-and second-kinds.Other contributions were presented, but other types of CPs were used.The third and fourth kinds of CPs are special non-symmetric Jacobi polynomials.The authors of [24] solved some variable-order fractional DEs by using the shifted CPs of the third kind; the authors in [25] treated some fractional DEs.Some singular Lane-Emden equations were treated using the third-and fourth-kinds in [26].Several applications have recently been utilized by the other four types of CPs, classified as special ones of the generalized Gegenbauer polynomials ( [27,28]).For example, the authors in [29] solved some mixed-type fractional DEs using the fifth-kind CPs.The authors of [30] employed the sixth-kind CPs to treat some telegraph type-equations.More recently, the authors of [31] introduced and utilized the seventh-kind CPs to solve some fractional delay DEs.Furthermore, the polynomials, namely, the eighth-kind CPs, solve the generalized Kawahara DEs using two different approaches in [32,33].
Due to their incredible accuracy and computing efficiency, spectral methods have become famous for solving DEs [34].Many contributions employ the different versions of spectral methods; see, for example, [35][36][37][38][39][40].Among these techniques for fractional DEs is the spectral tau approach, which converts the original problems into algebraic systems that can be solved computationally effectively [41,42].In this study, we suggest a new use of the spectral tau approach for solving the time-fractional DEs, with basis functions being seventh-kind CPs.
Since CPs are a desirable option for spectral approximations because of their rapid convergence, numerical stability, and orthogonality, we create a strong basis expansion designed explicitly for solving TFDE.Our main goal is to provide accurate numerical estimates for TFDE using the spectral tau approach with seventh-kind CPs.We will investigate the theoretical and numerical aspects of our method.The suggested method's accuracy and convergence behavior will be carefully studied.Moreover, we evaluate our methodology's effectiveness and computational power using a set of numerical examples that demonstrate its faithful representation of the complex dynamics of TFDE solutions.
This paper is organized into the following sections: Sect. 2 presents some fundamentals and useful preliminaries.A tau approach for solving the TFDE is followed in Sect.3.
A thorough error analysis of our methodology is provided in Sect. 4. Three examples are provided in Sect. 5. Lastly, Sect.6 reports a few conclusions.

Some key relations and formulas
This section is dedicated to reviewing the fundamentals of fractional calculus.Moreover, some relations and properties of the seventh-kind CPs and their corresponding shifted ones are exhibited.

A account on seventh-kind CPs
Consider the seventh-kind CPs { Z (ξ ), ≥ 0} defined as [31] Z where (a) n is the Pochhammer function defined as (a) n = (a+n) (a) , and P (a,b) (ξ ) is the classical Jacobi polynomial of degree .
Z (ξ ) ≥0 is an orthogonal set on [-1, 1] with the following orthogonality relation: , if even, with Lemma 1 [31] The following trigonometric expressions apply: Lemma 2 [31] Let ∈ N 0 .Zj (ξ ) may be expressed as: Theorem 1 [31] Let ∈ N 0 .ξ 2 and ξ 2 +1 can be expressed in terms of Zj (ξ ) as follows: Theorem 2 [31] It is possible to express Zj (ξ ) as a combination of the first-kind CPs T i (ξ ) as Theorem 3 The second derivative of Z (ξ ) can be expressed as Proof First, we prove (17).The series expression of Zi (t) in (10) leads to the following expression for d 2 Z2 (ξ ) which, by utilizing the inversion formula (1), can be expressed as Upon expanding and rearranging the terms, the last relation can be transformed into Now, setting By directly applying Zeilberger's algorithm [43], we can derive the subsequent recurrence relation: with the initial values: The recursive formula (23) has the following exact solution Therefore, relation ( 17) may be obtained.
To prove formula (18), we use formula (11) to get Using ( 13), we get that may be turned into Now, set Apply Zeilberger's algorithm to demonstrate that the function Gk, fulfills the subsequent recurrence relation: with: The solution of ( 28) is Therefore, relation ( 18) may be obtained.

A account on the shifted CPs of seventh kind
The shifted seventh-kind CPs: According to (5), the polynomials C n (t), n ≥ 0, are orthogonal on [0, τ ], with τ 0 Theorem 4 [31] The analytic form of C i (t) is given as follows: where Theorem 5 [31] t m can be expanded in terms of C i (t) as where

Corollary 1 It is possible to express C j (t) as a combination of the shifted first-kind CPs T
Proof The proof is easily obtained by replacing ξ with ( 2 t τ -1) in Theorem 2.
Theorem 6 For 0 < α < 1, we have where Proof From Eq. (31), we have using (3), the last equation can be rewritten as In virtue of the inversion formula (33), Eq. ( 40) becomes The arrangement of Eq. ( 41) leads to the following formula which can be written in simple form as This completes the proof of the theorem.
Consider the space and assume that u(ξ , t) ∈ L 2 ω(ξ ,t) ( ) is expanded as where ω(ξ , t) . It is easy to see that b ij are given by Therefore, u(ξ , t) can be approximated in where and B = (b ij ) 0≤i,j≤M is a (M + 1) × (M + 1) matrix.Now, the following residual R(ξ , t) of Eq. ( 43) can be obtained after replacing u(ξ , t) with u M (ξ , t) Let us define where ω(ξ ) = ξ 4 1-ξ 2 and ω(t) = The application of tau method [30] leads to which leads to the following matrix equation where

Theorem 7
The elements of matrices A, E, G and Y are given by Proof By using the orthogonality relations ( 5), (30) and Theorem 3, the entries of A, E, and G may be easily found.By virtue of Theorem 6, the elements y j,s can be written as With the aid of formula (31) The integration in the right-hand side of Eq. ( 62) can be evaluated to give the following result where β(ξ , y) is the beta function.
Step 2. Apply tau method to obtain the system in (55) and (64).
Step 3. Use Theorem 7 to get the elements of matrices A, E, G and Y Step 4. Use NDsolve command to solve the system in ( 55) and (64) to get b ij .
Moreover, we get the following initial and boundary conditions Merging the M 2 -M equations in (55) with the (3 M + 1) constraints in (64), a system of order (M + 1) 2 is obtained, and thus the approximate solution can be derived using the Gaussian elimination technique.
Remark 1 The Algorithm 1 specifies the necessary steps to obtain the desired approximate solution to the time-fractional diffusion equation.

Error analysis
This section is confined to investigating the double expansion used in approximation.The following two lemmas are useful for what comes next.

Lemma 4
The following inequality is satisfied in the interval [0, τ ] Proof The proof of Lemma 4 is similar to the proof of Lemma 3 after using |T * (t)| ≤ 1.
Lemma 5 Let ≥ 0. The following inequality holds Proof The proof of this lemma is similar to the proof of Lemma 3 after using Eqs.( 14) and ( 15) along with |T m (t)| ≤ 1.
Theorem 8 Assume that a function u(ξ , t) = h(ξ ) h(t) ∈ L 2 ω(ξ ,t) ( ), with h(ξ ) and h(t) have bounded fourth derivative satisfy the expansion in (47).Then the series in (49) converges uniformly to u(ξ , t).Moreover, the following bound on the expansion coefficients in (49) satisfy the inequality where n n means that there ia a constant c with n ≤ c n.
Proof From Eq. ( 48), we have The previous equation will be transformed after using the following substitutions and the assumption u(ξ , t) = h(ξ ) h(t) along with the definition of which can be written as where and Now, we have four cases: (i) If i, j even.Firstly, we will find I 1 .
Using Lemma 1, i = 9 π 8 (i+1) (i+3) in Eq. ( 6) along with the following property and simplifying the result yields If we integrate Eq. ( 74) by parts three times, we get the following relation where (cos((i + 3)ψ) -cos((i + 5)ψ)) (cos(iψ) -cos((i + 2)ψ)) Now, integrating Eq. ( 76) again by parts and taking into consideration that h(ξ ) and h(t) have bounded fourth derivative, then the following inequality is obtained In addition, after applying the trigonometric representations to I 2 , the integration by parts three times results in where And hence, we get the following inequality after integrating Eq. ( 79) by parts and using the hypothesis of the theorem.Finally, we have from Eqs. ( 71), ( 77) and ( 80) To prove the following three cases (ii) i, j odd.
(iii) i even, j odd.
(iiii) i odd, j even.Imitating similar steps to those given in case (i), we get the following inequality This completes the proof of this theorem.

Remark 2
The following inequalities will be obtained if we follow similar steps to those given in Theorem 8 and Lemma 6 [45] Consider a positive, decreasing, continuous function f Theorem 9 If u(ξ , t) satisfies the hypothesis of Theorem 8, then the following estimate of truncation error is fulfilled: Proof We can write Inserting Eqs. ( 82) and (83) into Eq.( 85) and using the inequality Using Lemma 6 together with the inequality where f is decreasing function, we get This completes the proof of Theorem 9.

Theorem 10
The truncation error estimate can be estimated as Proof First, we can write With the aid of Lemmas 4, 5, Theorem 8 and Eqs. ( 82), (83), the last equation becomes (89) Now, noting the inequality: we get Inserting Eqs.(90) into Eq.(89), we get Theorem 11 (Stability) Under the assumptions of Theorem 8, we have Now, with the aid of Theorem 8, Eqs.( 82) and (83) and using similar steps to those followed in the proof of Theorem 9, we get the desired result.

Illustrative examples
The following section presents three examples to ensure the validity and accuracy of our presented algorithm.
Example 1 Consider the following equation governed by with the exact solution: u(ξ , t) = t 4 ξ 4 .
Figure 1 shows the absolute errors (AE) at different values of α at M = 4. Table 1 shows the AE at α = 0.5 and M = 4.In addition, Table 2 shows the L ∞ and L 2 errors for α = 0.6 at M = 4 and M = 5.The approximation of the exact solution is seen to be close.governed by

Example 2 Consider the following equation
where E p,q (t) is the Mittag-Leffler function and the exact solution of this problem is u(ξ , t) = ξ 3 e α t .In Fig. 2, we plot the AE at different values of α when M = 18.These figures show the fast reduction of AE at small values of M to achieve very high accuracy.Table 3 shows the L ∞ and L 2 errors for α = 0.5 at M = 18 and M = 19.These results show the accuracy of our method.
Example 3 Consider the following equation  governed by where E p,q (t) is the Mittag-Leffler function, and the exact solution of this problem is u(ξ , t) = ξ 2 sin(π t).
In Fig. 3, we sketch the AE (left) and Approximate solution (right) for α = 0.5 at M = 18.
In Tables 4 and 5, we list the L ∞ and L 2 errors for different values of α at M = 18.Table 6 presents the maximum AE at α = 0.3 and different values of M.These findings demonstrate how closely the approximate solution approaches the exact solution.
where E p,q (t) is the Mittag-Leffler function and the exact solution of this problem is u(ξ , t) = sin(π ξ) sin(π t).
Table 7 shows a comparison of the best L 2 errors between our at M = 18 method and method in [44] at M = 11, N = 15 when α = 0.9.Moreover, Fig. 4, illustrates the AE (left) and approximate solution (right) for α = 0.1 at M = 18.Table 8 presents the maximum AE for α = 0.3 and different values of M.These results prove that the approximate solution is very close to the exact solution.
Remark 3 In light of Figs.1-4, we can conclude the following benefit: Excellently precise approximations can be obtained by selecting modified sets of seventh-kind CPs as basis functions and taking specific terms of the retained modes.Our method at M = 18 10 -13

Concluding remarks
In conclusion, our study introduces a novel application of the spectral tau methodology tailored to tackle the TFDE by leveraging seventh-kind CPs as basis functions.The choice of CPs stems from their rapid convergence, numerical stability, and orthogonality properties, allowing for constructing a robust basis expansion explicitly designed for solving TFDE.Through a comprehensive error analysis and numerical validation, we have demonstrated our proposed method's efficacy and computational prowess in accurately capturing the complex dynamics of TFDE solutions.Moving forward, our work opens avenues for further exploration and refinement of spectral techniques for solving TFDE, offering a promising approach with diverse applications across scientific and engineering domains.Continued efforts in advancing numerical methods and innovative algorithmic developments can significantly enhance our understanding and computational capabilities in modeling complex dynamical systems governed by fractional differential equations.

Figure 1
Figure 1 The L ∞ error of Example 1

Figure 2
Figure 2 AE of Example 2

Figure 4
Figure 4 The AE (left) and Approximate solution (right) of Example 4

Table 1
AE of Example 1

Table 2 L
∞ and L 2 errors of Example 1

Table 3 L
∞ and L 2 errors of 2

Table 4 L
∞ of Example

Table 8
The maximum AE of Example 4