Chebyshev Pseudospectral Method for Fractional Differential Equations in Non-Overlapping Partitioned Domains

: Fractional differential operators are inherently non-local, so global methods, such as spectral methods, are well suited for handling these non-local operators. Long-time integration of differential models such as chaotic dynamical systems poses specific challenges and considerations that make multi-domain numerical methods advantageous when dealing with such problems. This study proposes a novel multi-domain pseudospectral method based on the first kind of Chebyshev polynomials and the Gauss–Lobatto quadrature for fractional initial value problems.The proposed technique involves partitioning the problem’s domain into non-overlapping sub-domains, calculating the fractional differential operator in each sub-domain as the sum of the ‘local’ and ‘memory’ parts and deriving the corresponding differentiation matrices to develop the numerical schemes. The linear stability analysis indicates that the numerical scheme is absolutely stable for certain values of arbitrary non-integer order and conditionally stable for others. Numerical examples, ranging from single linear equations to systems of non-linear equations, demonstrate that the multi-domain approach is more appropriate, efficient and accurate than the single-domain scheme, particularly for problems with long-term dynamics.


Introduction
Fractional differential equations (FDEs) are gaining significant attention in engineering, science, finance and other fields.Unlike classical differential equations, which are limited to integer-order derivatives, FDEs introduce a higher degree of complexity into differential equations by incorporating arbitrary non-integer-order derivatives.Fractional differential operators can be defined in various ways, with the most common definitions being those involving singular kernels, such as the Riemann-Liouville and Caputo fractional differential operators [1,2].These definitions rely on integrals with power-law kernels, which can sometimes lead to singularities and complicate numerical approximations.These differential operators with power law kernels are, however, advantageous for developing spectral-based methods using basis functions such as the Chebyshev and Legendre polynomials, as these polynomials can be expressed in power series form.Therefore, the Caputo derivative will be used in the numerical method proposed in this study.In recent years, however, there has been a significant interest in non-singular kernels, which provide an alternative framework for defining fractional operators.Examples include the Caputo-Fabrizio and Atangana-Baleanu fractional operators, which use the exponential and Mittag-Leffler functions as kernels, respectively [3,4].These non-singular kernels were introduced to broaden the scope and applicability of fractional calculus in various fields.This added complexity allows fractional differential operators to capture the memory effects of a system or describe memory-inherent phenomena.As a result, FDEs offer a more robust framework for modeling the intricate dynamics of complex systems across various scientific disciplines [5,6].Finding analytical solutions for fractional differential equations (FDEs) is particularly challenging; however, several studies have proposed analytical solutions and methods for FDEs, for instance, Ahmadian and Darvishi [7] derived solitons and traveling wave solutions using the tanh-coth, tan-cot, sech-csch and sec-csc expansion methods for the (1 + 1)-dimensional fractional Biswas-Milovic equation [8].Kumar et al. [9] proposed a modified Laplace decomposition method, which hybridizes the Laplace transform and Adomian decomposition methods to analyze the fractional Navier-Stokes equation that models the unsteady flow of a Newtonian fluid in a tube.The fractional symmetric regularized long wave equation has been resolved analytically using methods such as the G ′ /G expansion method [10], the G ′ /G 2 expansion method [11], the extended complex method [11], the direct algebraic method [12] and the Jacobi elliptic function method [13].Other analytical methods that have been used for FDEs include the Sumudu transform [14], the Kudryashov method [15] and various semi-analytical techniques such as the exp expansion method [16,17] and the sine-cosine method [18], among others.Numerical methods also play an important role in obtaining accurate solutions for FDEs, and they provide valuable insight into the dynamics of systems described by FDEs [19,20].
Many numerical methods for approximating the solutions of FDEs have been proposed.For example, Hattaf et al. [21] introduced a novel numerical approach for approximating the generalized Hattaf fractional derivative (GHF), which involves a fractional differential operator with a non-singular kernel.The numerical method developed was tested on FDEs with analytical solutions to verify its accuracy by comparing the approximate solutions to the exact solutions.Non-linear FDEs describing the dynamics of HIV infection were also solved to demonstrate the applicability of the method to real-world problems.The predictor-corrector methods (PCMs) for FDEs were proposed by Kumar et al. [22], Jhinga and Daftardar-Gejji [23] and Lee et al. [24].These studies demonstrated that PCMs solve FDEs faster and drastically reduce computational time.Kumar et al. [25] introduced two computationally efficient numerical techniques based on the mid-point and Heun method for linear and non-linear fractional initial value problems (IVPs).The study concluded that the proposed methods offer significant advantages over the traditional and improved Euler methods.
The non-local nature of fractional differential operators sometimes makes analysis and finding solutions of FDEs herculean tasks [26].Spectral methods have been a powerful numerical technique for approximating the solutions of FDEs [27][28][29].Unlike other traditional numerical techniques, spectral-based methods are global and highly accurate, making them efficient for solving differential equations [30].Numerous studies have leveraged the global nature of spectral-based methods to discretize fractional differential operators and obtain numerical solutions of FDEs [31].Zayernouri and Karniadakis [32] and Zhao et al. [33] proposed an exponentially convergent Galerkin spectral method for FDEs.Zayernouri and Karniadakis [32] compared the performance of the developed spectral-based method and the established finite difference method (FDM) and concluded that the proposed method performs better than the FDM in terms of convergence rate.Khader and Sweilam [34] used the Legendre pseudospectral method to solve the fractional advection-dispersion equation with the Caputo fractional differential operator.The study analyzed the accuracy of the new discrete fractional differential operator by examining its convergence rate and estimating the maximum possible error, and the results were consistent with existing results from the literature.Stochastic fractional differential equations were solved by Cardone et al. [35] using the spectral collocation method.The method involves obtaining a fully discrete version of the original continuous operator and solving a system of non-linear algebraic equations.Tests on several problems demonstrated that the proposed method was effective for the chosen basis functions and collocation points.Spectral methods exhibit exponential convergence for differential equations with sufficiently smooth solutions.This high level of accuracy has made spectral methods, particularly the Chebyshev pseudospectral method, a popular choice for solving a variety of fractional differential equations.For instance, Doha [36] used the Chebyshev pseudospectral method to approximate the solutions of fractional multi-term differential equations.Similarly, Sweilam [37] applied this method to numerically solve fractional integro-differential equations.Other studies that have used the Chebyshev pseudospectral method to approximate the solution of fractional differential equations include those of Dabiri [38], Oloniiju et al. [39,40,41], among others.However, this desirable property of spectral convergence can deteriorate with large time domains or non-smooth solutions [42].The multi-domain approach tackles these challenges by offering an accurate and computationally efficient way to manage potentially non-smooth solutions.It can also improve numerical stability and accuracy by breaking down a problem's domain, especially time-dependent problems, into smaller sub-domains [43].Samuel and Motsa [44] described a multi-domain spectral collocation method for solving hyperbolic partial differential equations (PDEs) with a large time domain.The study showed that dividing the computational domain into smaller sub-domains gives accurate results and reduces the computational time.Birem and Klein [45] presented a multi-domain spectral method for the Schrödinger equation.Jung and Olmos-Liceaga [46] introduced a non-conforming and non-overlapping multi-domain spectral method to simulate traveling wave solutions in the reaction-diffusion equations.The method has three key features: high-order accuracy and convergence rate, adaptive mesh and the multi-domain technique.Qin et al. [47] proposed a multi-domain Legendre-Galerkin least square method for linear differential equations with variable coefficients.The study established two key properties of the method: continuity, ensuring consistency and smoothness at the endpoints of each sub-domain, and coercivity, ensuring stability and convergence.Additionally, the study established the optimal error estimate in the H 1 norm.
This study explores the Chebyshev pseudospectral method in a multi-domain setting as a novel approach to solving FDEs.Here, the domain of the problem is divided into a finite number of non-overlapping sub-domains.The multi-domain approach for FDEs has been applied to problems in biology [48] and finance [49,50].Maleki and Kajani [48] used the multi-domain Legendre-Gauss pseudospectral method to approximate solutions for the fractional Volterra population growth model.The study demonstrated the effectiveness of the proposed method in tackling a wide range of integro-differential equations, including those with integer and fractional orders.Bambe Moutsinga et al. [49] developed a spectral-based method for fractional differential equations which models valuation within an affine jump-diffusion framework, chaotic and hyper-chaotic systems and fractional pricing models for cryptocurrencies.The study established that the developed spectral method outperforms the MATLAB numerical packages based on the Euler method.Kharazmi [50] investigated a data-driven approach to solving stochastic fractional partial differential equations (SFPDEs) using multi-domain spectral methods.Zhao et al. [51] proposed a multi-domain spectral collocation method (MDSCM) for solving non-linear FDEs.The MDSCM exhibited a significant advantage over the single-domain spectral methods in achieving high accuracy for differential equations with low-regular solutions.Klein and Stoilov [52] proposed a numerical technique for the fractional Laplacian, based on the Riesz fractional operators, defined on R. The integrands are transformed using the underlying Z q curve to ensure accurate computation within each sub-interval.The transformed integrals within each sub-interval are then solved using the Clenshaw-Curtis algorithm.Tavassoli Kajani [53] proposed a modified Müntz-Legendre collocation technique on a single domain and in a multi-domain setting for the fractional pantograph equation.Xu and Hesthaven [54] proposed a highly accurate and stable multi-domain spectral method with a penalty method for fractional partial advection and diffusion equations.The study concluded that the proposed scheme's high accuracy makes it well suited for long-term integration of FDEs.It was also noted that the convergence rate aligns with the theoretical analysis of the scheme.
The present study aims to introduce a novel numerical method, the non-overlapping multi-domain Chebyshev pseudospectral method, based on the first kind of Chebyshev polynomials and the Gauss-Lobatto quadrature for fractional differential equations.The choice of Chebyshev polynomials is due to their orthogonality property with respect to their weight function, which simplifies the projection of the approximate solution onto the polynomial space.Additionally, the computation of Chebyshev-Gauss-Lobatto quadrature weights and nodes is quite straightforward, and this reduces computational overhead compared to other orthogonal polynomials.The proposed method works by decomposing the domain of a fractional differential equation into smaller, non-overlapping sub-domains.The novelty of this study lies in its approach of partitioning the domain of a fractional initial value problem into non-overlapping sub-intervals and projecting the solution of the problem onto the space of Chebyshev orthogonal polynomials of the first kind.This technique involves separating the fractional differential operator into the sum of the 'local' and 'memory' components within each sub-domain.To demonstrate the versatility and efficiency of the proposed method, stability analysis and numerical experimentation on various types of fractional differential equations are carried out.Numerical error analysis is conducted to demonstrate the advantages of the domain decomposition approach over the single-domain Chebyshev pseudospectral method.
The rest of the article is organized as follows: Section 2 provides certain preliminary definitions which are necessary for developing the numerical schemes; Section 3 is devoted to developing and analyzing the numerical scheme for FDEs using the Chebyshev polynomials on the Gauss-Lobatto quadrature in non-overlapping partitioned domains; Section 4 demonstrates the accuracy and effectiveness of the numerical method by presenting some examples ranging from linear single equations to systems of non-linear equations that require long-term integration; and Section 5 ultimately concludes the study.

Preliminaries
This section introduces and defines the basic concepts and fractional operators required to develop the numerical scheme.Definition 1. Suppose u : [a, b] → R is a continuous function, γ > 0 and the Riemann-Liouville fractional order integral of u(t) is defined as [55] RL a I where Γ is the Euler gamma function and t γ−1 * u(t) is the convolution product of the power function, t γ−1 , and the function u(t).
Definition 2 (Fractional integral of (t − a) j ).For the power function (t − a) j , the fractional integral is defined as [56] RL a I Definition 3. The left-sided Caputo fractional differential operator is defined as [55] C a D where ⌈γ⌉ is the ceiling of γ, ⌈γ⌉ ≥ γ and t ⌈γ⌉−γ−1 * u ⌈γ⌉ (t) is the convolution product of the power function, t ⌈γ⌉−γ−1 , and the ⌈γ⌉-th integer-order derivative of the function u(t).
Definition 4 (Fractional derivative of (t − a) j ).The fractional derivative of (t − a) j is defined as [56] C a D Definition 5.The set of collocation points used in this study is that of the Chebyshev-Gauss-Lobatto points defined as t r = − cos rπ N ; r = 0(1)N [57].

Numerical Discretization and Analysis via the Pseudospectral Method
This section introduces the numerical discretization of a fractional initial value problem using Chebyshev polynomials of the first kind and the Gauss-Lobatto quadrature.It covers the approximations of the fractional derivative of an arbitrary function, both in a single domain and in partitioned, non-overlapping domains.Firstly, the polynomial approximation of an arbitrary function u(t), continuous on an interval [0, h], is presented in terms of shifted first-kind Chebyshev polynomials.These polynomials are obtained through the following recurrence formula: The function u(t) is represented as follows: where the coefficients ûj in Equation ( 6) satisfy the ℓ 2 orthogonality condition of the shifted first-kind Chebyshev polynomials, which are defined with respect to the inner product where If a truncated N-th order approximation of the function is considered and collocated at N + 1 distinct Chebyshev-Gauss-Lobatto (CGL) points, with t r ∈ [0, h]; r = 0(1)N, then where ϱ i is the Christoffel number associated with the Chebyshev-Gauss-Lobatto quadrature and is defined as follows: Using the following matrix representations Equation ( 9) can be expressed as where

The Discretization of the Fractional Differential Operator in a Single Domain
Now, let us approximate the arbitrary non-integer order of the function u(t).Assuming that u(t) ∈ C p ([0, h]), where p ≥ N + ⌈γ⌉ + 1, representing the space of continuously differentiable functions, t > 0 and γ > 0, then the approximation for the arbitrary noninteger derivative of u(t) is given by The fractional-order derivative of the Chebyshev polynomial C 0 D γ t Π h,j (t) is obtained by using the series form of the polynomial [56], as follows: Proceeding to approximate t k−γ as a linear combination of shifted first-kind Chebyshev polynomials at the CGL points, Equation ( 14) can now be expressed as [39]: where h D (γ) = ∆ Ω Φ ∆ T ϑ and ∆, Ω and Φ are matrices whose entries are defined as follows: where 0,h D γ = ρ Π T h ϑ h D (γ) Π h .We refer to ( 0,h D γ ) T as the fractional differentiation matrix in the domain [0, h].

Discretization of Fractional Differential Operator in Non-Overlapping Partitioned Domains
Assume that u(t) is defined in t ∈ [0, T], partitioned into q non-overlapping subintervals, such that 0 = t 0 < t 1 < . . .< t q = T and h The first term in Equation ( 18) is the 'memory', and the latter term of the equation is the local arbitrary-order derivative of u(t) with respect to t ∈ [t l , t l+1 ].The arbitrary-order derivative in the first sub-domain is local and calculated with no memory term.The two terms in Equation ( 18) are approximated separately.If the length, h l , of each sub-domain, [t l , t l+1 ], is equivalent to h (from discretization in a single domain) in Section 3.1, and we collocate at N + 1 distinct CGL points, t l r ∈ [t l , t l+1 ], with t l N ≡ t l+1 0 , then the second term in Equation ( 18) is approximated as follows: The first term in Equation ( 18), the 'memory', which is the arbitrary non-integer-order derivative, C t e−1 D γ t e u(t), is approximated as The derivative C t e−1 D γ t e Π h l ,j (t) at the CGL points t l r is approximated as C t e−1 D γ t e Π h l ,j (t l r ) ≊ ∆ Θ, where the entries of Θ are defined as The integral in Equation ( 21) can be computed using any quadrature formula, such as Newton-Cotes or Gaussian quadrature formulas.Therefore, the first term of Equation ( 18) is then calculated as where where κ ∈ C and γ ∈ (0, 1).Suppose the domain [0, T] is partitioned into q non-overlapping domains, such that 0 = t 0 < t 1 < t 2 < . . .< t q = T, and as before, we define h l = t l − t l−1 , l = 1(1)q.To obtain the solution of Equation ( 23) on an arbitrary domain [t l , t l+1 ], we rewrite the equation as Assume that the solution u(t) of Equation ( 24) can be approximated as a linear combination of functions which are in the space of shifted Chebyshev polynomials of at most degree N on the CGL quadrature.Then, using the discretizations of the differential operator given in Sections 3.1 and 3.2, the FIVP ( 24) is now equivalent to the following algebraic systems: and where f = [ f (t l 0 ), f (t l 1 ), . . ., f (t l N )] T and l = 1(1)q.

The Numerical Stability of the Scheme
The stability of the proposed numerical method is demonstrated by considering the following linear fractional initial value problem (27) whose exact solution is u(t) = u 0 E γ,1 (κt γ ), where E γ,1 (•) is the two-parameter Mittag-Leffler function.If we consider a partition of the domain [0, h] into q = h equal sub-domains, the discretization of the model Equation (27) through the proposed numerical method in an arbitrary sub-domain [t l , t l+1 ] ≡ [t l , t l + h l ] is obtained as Here, ( t e−1 ,t e D γ e ) T is the memory differentiation matrix as in Equation ( 22), ( 0,h l D γ l ) T is the fractional differentiation matrix in a single domain [0, h l ] and To study the stability of the proposed pseudospectral numerical scheme, we define the stability matrix where w = κh γ l and P, Q and R are the matrices defined, respectively, as The stability of the numerical scheme is determined by the eigenvalues of S(w).The stability region of the numerical scheme ( 28) is typically defined by the set of points in the C-plane for which the spectral radius or the modulus of dominant eigenvalue is less than or equal to 1.If there is more than one non-zero eigenvalue with the potential to cross the stability threshold of 1, the stability region is then determined by considering these eigenvalues and the regions where these eigenvalues have a magnitude less than or equal to 1.The overall stability region of the numerical scheme would then be the intersection of these regions.If the stability region includes Re(w) ≤ 0 (the entire left half of the C-plane), then the numerical scheme is said to be absolutely stable (A-stable) [58].If, in addition to A-stability, the limit of the dominant eigenvalue tends to zero as w tends to −∞, then the numerical scheme is said to be Lubich-stable (L-stable) [59].The L-stability property means that any components that may cause instability in the numerical scheme are damped as the length of the sub-domain, h l , increases, which is particularly important when handling longterm integration of stiff differential equations.To illustrate, Table 1 shows the dominant eigenvalue, stability region and stability property for a linear discretization (N = 1) of Equation ( 28) with two sub-domains (q = 2) for some selected values of fractional order γ.As shown in the table, the stability region of the numerical scheme is the exterior of some closed disks in C; the boundaries of these disks are illustrated in Figure 1.As indicated by the figure, the stability regions for γ = 0.01, 0.10, 0.21, 0.32, 0.43, 0.54, 0.65 satisfy the conditions for absolute and Lubich stability, and so the numerical scheme for these values of γ is absolutely and Lubich-stable.However, the numerical schemes for γ = 0.76, 0.87, 0.98 are not absolutely stable and, by extension, not L-stable.
In Figure 2, the stability regions for a quadratic approximation with q = 3 and q = 5 for γ = 0.21, 0.54, 0.98 are presented.The stability region in the figure is the area outside the elliptical disks, which shows that the stability regions for q = 3 and q = 5 sub-domains are the same.However, in both cases, a varying number of eigenvalues with the potential to cross the stability threshold are obtained.For q = 3, there is one non-zero dominant eigenvalue for γ = 0.21 and 0.54, and there are four non-zero eigenvalues for γ = 0.98 with the potential to cross the stability threshold of 1.For q = 5, γ = 0.21, 0.54, 0.98 all have eight non-zero eigenvalues that can cross the stability threshold.However, when the intersection of the resulting stability regions from these eigenvalues is taken, the same stability region with q = 3 is obtained for each illustrated γ value, as shown in Figure 2. The figure shows that the numerical schemes for γ = 0.21 and γ = 0.54 are absolutely stable, and the scheme for γ = 0.98 is not absolutely stable (conditional stability).Although a general conclusion about the stability of the numerical scheme cannot be made, it can be inferred from this analysis that whether the numerical scheme is absolutely or conditionally stable depends on the value of the arbitrary non-integer order of the fractional differential equation.

Numerical Experimentation and Results
This section presents a series of numerical examples, ranging from single linear equations to systems of non-linear equations, along with the results.Detailed analyses of the results for each example are provided to illustrate key observations.For a comprehensive understanding of the numerical methods used in these examples, refer to the description provided in Section 3.
Table 2 presents the absolute error between the approximate and exact solutions using the non-overlapping multi-domain pseudospectral method with the CGL quadrature and the Chebyshev polynomials as the basis functions.The table also includes numerical results obtained without decomposing the domain.This comparative analysis highlights the accuracy and effectiveness of the multi-domain approach in reducing errors and improving the accuracy of the numerical solution.From the table, it can be seen that, for the fractional order of γ = 0.5, although the maximum absolute error in the multi-domain setup does not decrease in terms of significant zeros, it reduces as q increases.This indicates that, although minimal, there is an advantage in using more sub-domains.Comparing the error to the single-domain result, there is an additional significant zero in the error of the multi-domain approach, which indicates that the multi-domain approach is more effective for this problem.From Table 2, the effectiveness of the proposed numerical technique in obtaining an accurate approximate solution to a linear FDE is evident, so in the examples that follow, we consider non-linear FDEs to further validate the numerical method.

Example 2. Consider the non-linear fractional differential equation
where the initial condition is u(0) = −2 and the exact solution when γ = 1 is given by u .
Table 3 shows the approximate solution obtained using two different sub-domains in the proposed method.The solutions in these sub-domains match the numerical values reported in Odibat and Momani [60].It should be noted that closed-form solutions for γ = 0.2 and γ = 0.4 are not available.However, for γ = 1, the accuracy of the proposed method is unquestionable.To further give credence to the accuracy of the proposed method, Table 4 demonstrates that the residual error is consistently close to zero for various values of γ.This suggests that the proposed method effectively resolves the problem in Example 2.
Table 3. Approximate solution of Example 2 using two different sub-domains, q, and the numerical values reported in Odibat and Momani [60].
with the initial condition u(0) = 1 and whose exact solution is u(t) = t 2 − t + 1.
Table 5 presents the maximum absolute error of the difference between the exact and approximate solutions obtained using the multi-domain and single-domain methods.It is observed from the table that there is conformity with the trend observed in Table 2.As the domain of the FDE becomes larger, the error in the single domain scheme deteriorates; however, in comparison, the magnitude of the errors is better for the multi-domain scheme.We note that the magnitude of the error increases as γ increases, similar to the trend observed in Table 2.
If we replace the right-hand side of Equation ( 34) with 2 + 4 t π + t 2 and set t ∈ (0, 1], subject to the starting conditions u(0) = 0.0 and u ′ (0) = 0, the closed-form solution is given as u(t) = t 2 when γ = 1.5.Table 6 presents the errors obtained by taking the absolute value of the difference between the numerical solutions obtained through the multi-domain technique and the exact solutions, and these are compared with the errors obtained through the evolutionary artificial neural networks method in the study of Raja et al. [62].Table 6 compares the exact and approximate solutions obtained using the multidomain method at various t values and the associated errors at each t.The errors are compared with those obtained in Raja et al. [62].It is observed from the table that the approximate solutions match the exact solutions better than the results in Raja et al. [62], and this validates the efficiency of the proposed method as it outperforms a neural network approach known to be efficient and reliable.
To further validate the accuracy of the proposed multi-domain method, the infinity norm of the residual for different numbers of sub-domains is presented.Table 7 gives the infinity norm of the residual errors for Example 4. The residual error is obtained through the discrete equivalent of Equation (34), defined by Table 7 shows that the method is efficient in solving fractional-order differential equations of varying values of arbitrary non-integer order, as the residual error remains very small for various fractional orders, demonstrating the effectiveness of the multidomain technique.To further validate the results obtained, we reproduce the solution profiles obtained by Zafar et al. [63] and observe excellent agreement in the result, as shown in Figure 3.In the remaining examples in this section, we consider systems of coupled fractional initial value problems that describe chaotic systems.These equations often exhibit longterm dynamics and, invariably, require long-term integration to fully understand the underlying dynamics involved.Therefore, a numerical technique, such as the multi-domain Chebyshev pseudospectral method, that can handle long-term integration is a perfect fit for solving such systems, as corroborated by the results in the next four examples.
Example 5. Consider the fractional chaotic system of differential equations which describe a memory-inherent electrical system [64] C subject to the following initial state, (0.1, 0.0, 0.1).Here, c = 0.70, L = 3.00, β = 1.10 and α = 0.60.The voltage across the capacitor, c, is represented by the variable u 1 (t), and the current through the inductor, L, is represented by the variable u 2 (t), and the internal state of the memristive system is denoted by the variable u 3 (t).
The average maximum residual error for the system (36) is given in Table 8.We remark here that the maximum residual errors presented in Table 8 and subsequent tables in this section are obtained as the average of the maximum residuals of the equations in the system.The residuals are calculated similarly to Equation (35).It can be seen from Table 8 that for various fractional orders and seemingly large domains, [0, T], a high level of accuracy is obtained when the multi-domain Chebyshev pseudospectral method is used.Due to the inherently chaotic nature of the problem, the dynamics of the system typically require a large physical domain to fully develop, and to effectively capture these complex long-term behaviors, the physical domain is partitioned into a large number of smaller sub-domains.The table shows that the average of the maximum residual errors decreases marginally as the number of sub-domains increases.This inverse relationship between the number of sub-domains and the residual error highlights the applicability of the multi-domain technique when dealing with this kind of system that requires long-term integration.This is particularly evident when the error is compared to the error associated with the singledomain method.As already observed in the previous examples and now corroborated by Table 8, the error associated with the single-domain method is worse than the error obtained when the multi-domain approach is used.To lend further credibility to the numerical results of this example, we reproduce the dynamics, in the form of phase portraits, of the memristive chaotic electrical system as reported in the work of Sene [64].As shown in Figure 4, there is an excellent agreement between these phase portraits and those in the work of Sene [64].
under the initial state (0.2, −0.2, 0.2), with γ = 1.00,Ω 1 = 0.2, Ω 2 = 2.0, Ω 3 = 1.0 and Ω 4 = 0.1.Similar to the previous example, this system typically requires a large physical domain to fully stabilize and develop, and when solved on a single domain, the numerical solution is unstable, as evidenced by the average residual errors in a single domain presented in Table 9.We present the average maximum residual errors for various fractional orders and relatively large domains, [0, T], for different relatively large numbers of sub-domains in Table 9.There is an agreement between the level of accuracy observed for the results in Table 9 and those obtained in the previous examples.The error obtained for the various sub-domain numbers used is shown to be small, indicating that the proposed multidomain Chebyshev pseudospectral method efficiently solves this system of FDEs.As q increases, it is observed that the average residual error improves, resulting in higher accuracy.The superior effectiveness of the multi-domain approach over the single-domain numerical scheme becomes marginally clearer in this example.In this case, the singledomain numerical scheme is observed to produce quite large numerical errors, whereas the multi-domain technique can achieve a higher accuracy level.This disparity in the level of performance highlights the limitation in the scope of the single-domain scheme, while the multi-domain approach demonstrates robustness and versatility, especially for complex systems like the one considered in this example.
Figure 5 gives the behavior of the solution profiles as t evolves.The results are as expected and in complete agreement with the result of Azam et al. [65].It is also observed in Figure 6 that the chaotic system flows through the various states by turning and twisting about the quadrants.This dynamics was also captured by Azam et al. [65] and validates the numerical results being put forward in this study.(c) Figure 5.The time evolution of u 1 (t), u 2 (t), u 3 (t) of the chaotic system (37) with γ = 1.00,Ω 1 = 0.2, T = 100, N = 13 and q = 350.
Table 10 shows the average maximum residual errors for Equation (38) for selected values of arbitrary non-integer order, γ, and for two domains, (0, T].The results are presented for different numbers of domain decomposition.From Table 10, we see that the residual errors are small, indicating that the approximate solutions are accurate.As has been validated several times in this study, an increase in the number of sub-domains used results in more accurate solutions.It is also shown that using the single-domain approach results in less accuracy in the approximate solutions obtained.Figure 7 displays the phase portraits for the chaotic system, and the result agrees with the results obtained in Li et al. [66], thereby validating the present numerical solutions. with the initial state given as (0.2, 0.2, 0.2), Ψ 1 = Ψ 2 = Ψ 3 = 1.0, m 1 = 0.0 and γ = 0.9.
In Table 11, the average maximum residual errors obtained for Equation (39) are presented.The table shows that an increase in the number of partitioned domains results in higher accuracy of approximate solutions obtained for the selected fractional orders and domains.The average maximum residual errors are observed to be small, confirming that the proposed multi-domain method is highly efficient for this fractional differential model, unlike the single-domain approach, which is not as accurate.We present the phase portraits of the chaotic system in Figure 8.The figure validates the dynamics reported by Sene [67] and is in excellent agreement with the previously published result.The agreement between the phase portraits presented here and those documented in the study by Sene [67] serves to further corroborate the accuracy and reliability of the current numerical method.Figure 9 shows the phase portraits of the chaotic system in Equation (39) using two different γ values.The figure justifies that the dynamics of the chaotic system depend on the arbitrary real order.36)- (39).The figures show that the accuracy of the solutions obtained through the multi-domain Chebyshev pseudospectral method is high and sometimes reaches errors smaller than those reported so far in the tables.As can be clearly observed, as the number of sub-domains in the computational domain increases, the residual error decreases correspondingly.In some instances, the residual error is reduced to the order of 10 −16 .However, the residual errors with a small number of sub-domains already show good accuracy, thereby showing that using the multi-domain approach proposed in this study is very efficient in solving fractional differential equations that typically require long-term integration.38) and (d) Equation ( 39) using γ = 0.8, T = 50, N = 10 and q = 100.

Conclusions
This study introduced a novel approach, the multi-domain Chebyshev pseudospectral method for solving fractional differential equations.The novelty of the method stemmed from the use of non-overlapping multi-domains, Chebyshev polynomials of the first kind and Gauss-Lobatto collocation points in developing the scheme, which was specifically adapted to handle fractional differential equations.The method leverages the global properties of spectral-based methods for discretizing non-local operators such as fractional differential operators, the exponential convergence of the Chebyshev polynomial-based pseudospectral method and the multi-domain approach's ability to maintain the accuracy of solutions to problems with long-term dynamics.The stability region of the multi-domain scheme was established through linear stability analysis.The stability analysis showed that obtaining an absolute and Lubich-stable numerical scheme for some values of arbitrary noninteger order is possible.Several numerical experiments involving linear and non-linear single FDEs and systems of non-linear FDEs were conducted.The examples were carefully chosen to span varying complexity levels, including fractional initial value problems with long-term dynamics, such as chaotic systems, which were shown to significantly benefit from a multi-domain numerical scheme.These examples have established that the proposed multi-domain numerical method gives highly accurate results.To validate the method, comparisons of the numerical results obtained in this study were made with analytical solutions, where available, and approximate solutions from the literature.It was found that the proposed method gave more accurate solutions than the approximate solutions from the literature.Where an exact solution was unavailable, residual error analysis was conducted to establish the accuracy of the proposed method.In all the examples considered, the method was found to produce comparative, highly accurate solutions.In future studies, we will explore the applicability of the Chebyshev pseudospectral method in non-overlapping or overlapping partitioned domains to fractional boundary value problems with non-linear dynamics.

Table 5 . 5 N = 10 , T = 6
The maximum absolute error of the difference between the exact and approximate solutions of Example 3 for different numbers of sub-domains, q.N = 10, γ = 0.

Figure 3 .
The dynamics of the Bagley-Torvik Equation(34) for various values of γ with T = 70, N = 10 and q = 400: (a) the phase portrait for the u − u ′ plane and (b) the u(t) profile.

Table 8 .
The residual errors for Example 5 for two values of T using N = 10 and for different numbers of sub-domains q, with c = 0.70, L = 3.00, β = 1.10 and α = 0.60.

Figure 10
Figure 10 displays the accuracy of the numerical results (the maximum absolute residual in each sub-domain) of each equation in the chaotic systems presented in Equations (36)-(39).The figures show that the accuracy of the solutions obtained through the multi-domain Chebyshev pseudospectral method is high and sometimes reaches errors smaller than those reported so far in the tables.As can be clearly observed, as the number of sub-domains in the computational domain increases, the residual error decreases correspondingly.In some instances, the residual error is reduced to the order of 10 −16 .However, the residual errors with a small number of sub-domains already show good accuracy, thereby showing that using the multi-domain approach proposed in this study is very efficient in solving fractional differential equations that typically require long-term integration.

Figure 10 .
The maximum residual error in each domain for (a) Equation (36) (b) Equation (37) (c) Equation ( T , t e i are the CGL points mapped onto [t e−1 , t e ] and t e−1 ,t e D γ e = ρ Π T

Table 1 .
The dominant eigenvalue, stability region and stability property for some selected values of arbitrary order, γ, for a linear approximation (N = 1) with q = 2.

Table 4 .
The residual errors of Example 2 at selected points using two different sub-domains, q, for different values of fractional order γ.

Table 6 .
[62]absolute errors of the difference between the exact and numerical solutions of Example 4 at different t values using N = 12 and q = 400 in comparison with the errors reported in the work of Raja et al.[62].

Table 7 .
The maximum residual errors for Example 4 for different T values, fractional orders, and different numbers of multi-domains, q.

Table 9 .
The average maximum residual errors for Example 6 for two values of T and different numbers of sub-domains q using N = 10.The following parameter values are used: Ω 1 = 0.2, Ω 2 = 2.0, Ω 3 = 1.0 and Ω 4 = 0.1.

Table 10 .
The residual errors for Example 7 at two final times T using the number of grid points N = 10 for different numbers of sub-domains q and Ψ 1 = 0.40.