A Highly Accurate Computational Approach to Solving the Diffusion Equation of a Fractional Order

: This study aims to present and apply an effective algorithm for solving the TFDE (Time-Fractional Diffusion Equation). The Chebyshev cardinal polynomials and the operational matrix for fractional derivatives based on these bases are relied on as crucial tools to achieve this objective. By employing the pseudospectral method, the equation is transformed into an algebraic linear system. Consequently, solving this system using the GMRES method (Generalized Minimal Residual) results in obtaining the solution to the TFDE. The results obtained are very accurate, and in certain instances, the exact solution is achieved. By solving some numerical examples, the proposed method is shown to be effective and yield superior outcomes compared to existing methods for addressing this problem.


Introduction
Fractional calculus, a branch of mathematical analysis dealing with derivatives and integrals of non-integer orders, has diverse applications across various fields.These applications showcase the versatility and significance of fractional calculus in diverse fields, ranging from control systems and biology to economics [1,2].By utilizing fractional calculus, researchers can create more precise models of intricate systems and detect dynamic behaviors that could be missed by conventional integer-order calculus.A wealth of information about engineering and physical processes, along with an extensive use of fractional-order derivatives, can be discovered in [3][4][5][6].
The diffusion equation serves as a foundational partial differential equation, outlining the process by which quantities like heat, mass, or momentum disperse throughout a medium as time progresses.The diffusion equation is a widely utilized tool in physics for turbulence [22], heat conduction [23], dissipation [24], magnetic plasma [25], and electron transportation [26].Anomalous diffusion is notable for its exceptional traits such as long-distance interaction and history dependency, which differ from typical diffusion phenomena.The conventional model based on integer-order differential equations struggles to accurately capture these anomalous behaviors.Instead, the fractional derivative has emerged as a viable alternative modeling technique for representing these anomalous diffusion phenomena [4].In recent years, interest has been increasing in investigating anomalous diffusion equations using fractional derivatives [27][28][29].
This paper focuses on solving the time-fractional diffusion equation in which c D µ t with 0 < µ ≤ 1 is the Caputo fractional derivative (CFD) operator with respect to variable t, q(x) ∈ C 1 [0, 1], r(x) ∈ C[0, 1], and q(x) > 0, r(x) ≥ 0 (∀x ∈ [0, 1]).Moreover, w is considered to be a smooth function, and g ∈ C([0, 1] × [0, T]).For this equation, the Dirichlet boundary and initial conditions are as follows: where the functions Analytical solutions for fractional-derivative diffusion equations are generally scarce, except for cases involving straightforward initial and boundary conditions [30].Therefore, the numerical solution method is crucial for solving the fractional derivative diffusion equation in practical scenarios.In [31], the authors obtained the approximate solution to a fractional advection-dispersion flow equation using finite difference approximation.A numerical scheme based on the random walk method is proposed in [23] for solving the considered equation.Chen et al. [32] proposed the Kansa method for solving the considered Equation (1).A paper focused on the kernel-based scheme to solve (1) is introduced in [33] by Fardi.The subsequent sections of this paper are structured in the following manner: Chebyshev cardinal polynomials and their properties are reviewed and introduced in Section 2. The pseudospectral method is applied to solve the TFDE (1) in Section 3. Section 4 is devoted to demonstrating the practicality and precision of the method.Section 5 of this paper provides a concise summary of the findings.

Chebyshev Cardinal Polynomials
Given N ≥ 0, let R := {r j : T N+1 (r j ) = 0, j ∈ N } be the set of the roots of the TChebyshev polynomial T N+1 in which N := {1, 2, . . ., N + 1}.Recall that the TChebyshev polynomials are defined on [−1, 1] by T N+1 (cos(θ)) = cos((N + 1)θ), N = 0, 1, . . .and their roots are specified by Shifted TChebyshev polynomials for generic intervals [a, b] are related to the TChebyshev polynomials by and the roots of T * N+1 in its turn are obtained by The Chebyshev cardinal function (CCF) is one of the orthogonal polynomials' most notable cardinal functions [34][35][36].Considering T * N+1,t (t j ) as the derivative of function T * N+1 (t) with respect to the variable t, Chebyshev cardinal functions can be denoted by The most striking feature of these polynomials is their cardinality, i.e., in which δ ji indicates the Kronecker delta.This property is mostly important as it enables us to approximate any function w ∈ H α ([a, b]) (the Sobolev space H α ([a, b]) will be briefly introduced) easily and without integration in finding the coefficients, viz, In what follows, since we need the definition of Sobolev spaces and their norm, we provide a brief definition of it.For α ∈ N, we denote by H α ([a, b]) the sobolev space of functions w(t) which have continuous derivatives up to order α such that and the semi-norm Lemma 1.Given N ≥ 0, if R * denotes the shifted Chebyshev nodes {t j } j∈N , then the error of approximation (8) can be bounded where the constant C is independent of N. Furthermore, it can be verified that Operational matrix of derivative Let Ψ(t) be a vector function with entries {ψ j } j∈N .We specify the operational matrix of derivatives for CCFs as To evaluate the elements of D, they can be obtained via the following process using the approximation (8).It follows from (8) that It is worth noting that there is another presentation of CCFs [37] where ϱ = 2 2N+1 /((b − a) N+1 T * N+1,t (t j )).When the operator D acts on both sides of ( 15), coming back to (14), we obtain by ( 8) It can be shown by ( 14) and ( 16) that

Operational Matrix of Fractional Integration
Considering the interval [0, 1], the fractional integral is defined as where Γ(µ) denotes the Gamma function.
Note that there is a square matrix I µ such that the acting of the fractional integral operator on Ψ(x) can be represented by it, viz, It is straightforward to show that the elements of this matrix can be obtained by After performing some simple calculations, it can be inferred from [38] that in which Motivated by (15), the CCFs can be determined by Using this definition of CCFs, (19) leads to So, it can be concluded that

Pseudospectral Method and Its Implementation
The present chapter is focused on solving the fractional diffusion equation (FDE) with the Caputo operator using an efficient and accurate scheme based on the pseudospectral method.As mentioned above, we consider the inhomogeneous TFDE with conditions ], and r(x) ≥ 0, q(x) > 0 (∀x ∈ [0, 1]).Also, we consider the unknown solution w to be an analytic function.Furthermore, the functions g(x, t), w 0 (x), f 1 (t), and f 2 (t) belong to the spaces C([0, 1] × [0, T]), C[0, 1], and C[0, T], respectively.To simplify the analysis, hereafter, we assume that T is equal to 1.
To obtain the pseudospectral discretization of Equation ( 27), the solution w(x, t) is approximated using CCFs as follows: where the (N + 1)-dimensional square matrix W consists of the unknowns {w i,j } N+1 i,j=1 .Substituting w N instead of w in (27), we have Now, we estimate each term in (31) as follows: • Taking into account the representation of CFD based on CCFs as the operational matrix D µ and (30), one can obtain • The first step in approximating the second term is to estimate ∂w N ∂x as Similarly, we can approximate ∂ 2 w N ∂x 2 as follows: Motivated by the second term, we have Substituting (33) and (35) into Equation ( 31) leads to the introduction of the residual function Selecting the roots of T * N+1 (Shifted Chebyshev polynomial) as the collocation points which are outlined in the introductory Section 2, and we denote them as (t i , t j ), i, j ∈ N , Equation (35) reduces to a linear system: R i,j = r(t i , t j ) = 0, i, j = 1, . . ., N.
In order to apply the boundary and initial conditions ( 28) and ( 29), we replace some equations of (37) by This leads to a new linear system Υ W = Ḡ, (38) in which Ḡ and W are the matrix-to-vector conversion by rows of the matrices G and W whose elements are By implementing the GMRES (generalized minimal residual) method [39], we obtain the matrix W.

Numerical Results
By providing some numerical simulations, the effectiveness of the present method is can showcased.These examples demonstrate how the method can provide practical solutions to various problems.To provide an overview of the method's efficiency, tables and figures report absolute error All examples were run on the combined use of Maple and Matlab software (version R2022a) with an Intel(R) Core(TM) i7-7700k CPU 4.20 GHz (RAM 32 GB).To obtain a higher precision, we increased the precision beyond 50 digits.
Example 1.As the first example, we utilized the presented scheme for IFDE (27) For this equation, the initial and Dirichlet boundary conditions were considered.According to [33], w(x, t) = sin(πx)(t 3 − t) is the exact solution.
Table 1 shows the CPU time and L 2 and L ∞ errors considering different values of µ and N. The results obtained showcase the capability and effectiveness of the method.The convergence of the proposed scheme is also verified via the presented results.To provide more evidence for the capability of the presented method, the results of this work are compared with those of a kernel-based method [33] in Table 2.The approximate solution and related absolute errors for µ = 0.2 are plotted in Figure 1.The values of the L ∞ -error and L 2 -error for µ = 0.2 versus N are illustrated in Figure 2.
For this equation, the Dirichlet boundary and initial conditions are considered.According to [33], the exact solution is w(x, t) = t(e −x − 1)(1 − x).
Table 3 shows the CPU time, L 2 -error, and L ∞ -error for different values of N. It is worth realizing that the error tends to zero as N → ∞.Moreover, Figure 3 is depicted to showcase the method's convergence.To compare the obtained results via the present scheme and the kernel-based method [33], Table 4 is tabulated.The L 2 -error at different times versus the number of bases N is reported in Table 5.   Example 3. To illustrate the accuracy of the method, the following equation from [33] is considered: For this equation, the initial and Dirichlet boundary conditions are considered.It follows from [33] that w(x, t) = t 2 x(2 − x) is the exact solution.
For this equation, the exact solution is obtained using the proposed scheme.To show this, the L 2 -error with N = 3 at different times is demonstrated in Table 6.
For this equation, the initial and Dirichlet boundary conditions are considered.The exact solution is w(x, t) = t 3/2 x(1 − x).
Table 7 shows the CPU time, L 2 -error, and L ∞ -error for different values of N and µ = 0.5.This table shows the convergence of the method.The approximate solution and related absolute error for N = 14 and µ = 0.5 are plotted in Figure 4.
For this equation, the Dirichlet boundary and initial conditions w(0, t) = w(1, t) = 0, w(x, 0) = 0, x, t ∈ [0, 1], are considered.Motivated by [40], the exact solution is w(x, t) = sin(πx)t 2 .Table 8 demonstrates the CPU time, L 2 -error, and L ∞ -error for different values of N and µ.To compare the presented method and the finite difference method, Table 9 is tabulated.It can be seen that the proposed method provides a better result with a lower computational cost.The approximate solution and related absolute error for µ = 0.7 are plotted in Figure 5.

Conclusions
The pseudospectral method is widely recognized as a highly effective and efficient technique for solving various equations.On the other hand, according to their inherent properties, Chebyshev cardinal polynomials are highly effective and powerful bases in numerical techniques.Hence, in this study, the spectral method based on Chebyshev bases is considered to address the time-fractional diffusion equation.Using the operational matrix of fractional derivatives in the Caputo sense, the desired equation is reduced to an algebraic linear system.The GMRES method is applied to solve this system.The results obtained are very accurate, and in certain instances, the exact solution is achieved.By solving some numerical examples, it is demonstrated that the proposed method is effective and yields superior outcomes compared to existing methods for addressing this problem.
In the future, this numerical approach to solving generalized fractional models will be expanded, including the space-fractional advection-diffusion equation [41] and timefractional diffusion equations with a time-invariant-type variable order [42], etc.

Figure 2 .
Figure 2. Plot of L 2 -error and L ∞ -error for different values of N (convergence confirmation) (Example 1).

Figure 3 .
Figure 3. Plot of L 2 -error and L ∞ -error versus the number of bases N (convergence confirmation) (Example 2).

Table 1 .
CPU time, L 2 -error, and L ∞ -error for different values of µ and N (Example 1).

Table 3 .
CPU time, L 2 -error, and L ∞ -error for different values of N (Example 2).
Example 4. Consider the following time-fractional diffusion equation: c D µ t (w)(x, t) =

Table 7 .
CPU time, L 2 -error, and L ∞ -error for different values of N (Example 4).
Example 5.The following IFDE is solved: c