Solution Method for Systems of Nonlinear Fractional Differential Equations Using Third Kind Chebyshev Wavelets

: Chebyshev Wavelets of the third kind are proposed in this study to solve nonlinear systems of FDEs. The main goal of the method is to convert the nonlinear FDE into a nonlinear system of algebraic equations that can be easily solved using matrix methods. In order to achieve this, we ﬁrst generate the operational matrices for the fractional integration using third kind Chebyshev Wavelets and block-pulse functions (BPF) for function approximation. Since the obtained operational matrices are sparse, the obtained numerical method is fast and computationally efﬁcient. The original nonlinear FDE is transformed into a system of algebraic equations in a vector-matrix form using the obtained operational matrices. The collocation points are then used to solve the system of algebraic equations. Numerical results for various examples and comparisons are


Introduction
Fractional differential equations (FDEs) can use real-valued orders for derivatives and integrals. FDEs applied to the evolution of physical processes reveal behaviors that classical differential equations do not allow. On the other hand, because of the additional complexity brought on by arbitrary orders of derivation and integration, it is extremely difficult to derive analytical solutions to many sorts of such FDEs. Finding precise and efficient numerical solution methods for the FDEs is therefore crucial. Recent years have seen the use of specific numerical techniques for FDEs, including Adomian decomposition method [1,2], predictor-corrector methods [3,4], finite difference method [5], Adams-Bashforth-Moulton method [6,7], F-expansion method [8], B-spline collocation method [9], reproducing kernel method [10,11], the homotopy perturbation transform method [12,13], and the residual power series method [14,15].
Similar research is being conducted on several wavelet types for various computationally demanding problems. Data can be divided into numerous time-frequency components using wavelet analysis. These functions are produced by stretching and moving a so-called mother wavelet function. The obvious advantage of the wavelet basis is that it simplifies the FDE problem's solution to the solution of a system of algebraic equations. The other benefits can be their orthogonality, compact support, and simultaneous representation of data in several resolutions. A wide variety of FDEs have been solved using numerous wavelet basis functions [16][17][18][19][20][21][22][23][24][25].
In this study, we attmpt to solve the system of FDEs in the form of: D α i * u i (x) = f i (x, u 1 , u 2 , . . . , u n ), u (r)

Third Kind Chebyshev Wavelets
Wavelets consist of a localized wave-like main function and sub-functions can be obtained from this function with the properties of zero-average and finite-energy. The main function is called the mother wavelet. Then, the shifted and stretched versions of this mother wavelet is obtained using a dilation parameter (a) and a translation parameter (b). Using the shifted and stretched versions of the wavelet function provides information on both the frequency content of the analyzed signal and where in time this frequency component occurs. The family of continuous wavelets is defined by, If parameters a and b are discrete as a = a −k 0 , b = nb 0 a −k 0 , a 0 > 1, b 0 > 0, the family of discrete wavelets are obtained as where n and k are positive integers. Discrete Chebyshev Wavelets of the third kind ψ nm (t) = ψ(k, n, m, t) can be defined on the interval [0, 1) as [25], where k is a positive integer, n = 1, 2, 3, . . . , 2 k−1 , and t denotes the normalized time. V m (t) are the Chebyshev polynomials of the third kind of the degree m, which are orthogonal with respect to the weight function ω(t) = 1+t 1−t on the interval [-1, 1] and satisfy the following recursive formula:

Block Pulse Functions
An m set of Block Pulse Functions (BPFs) is defined as where i = 1, 2, 3, . . . , m . The functions b i (t) are disjoint and orthogonal. They have the following properties for t ∈ [0, 1) Any function f (t) defined in [0, 1) with the property of squarely integrable in the interval can be expanded into an m set of BPFs as where The third kind Chebyshev wavelet matrix can also be expanded to an m set of BPFs as, The block -pulse operational matrix for fractional integration F α is defined as [27] where, with ξ k = (k + 1) α+1 − 2k α+1 + (k − 1) α+1 . Now we obtain the operational matrix of fractional integration for the third kind Chebyshev wavelet method (TKCWM): where the P α m ×m matrix of size m × m is called the TKCWM. Using (14) and (15) we obtain, Employing (14), (15), (17) and (18) results in The resulting operational matrix for TKCWM, denoted as P α m ×m yields This method yields a fractional integration operational matrix with a large number of zero entries, which speeds up the simulation. The error estimation of the third kind Chebyshev Wavelets and the convergence analysis can be found in [25].

Numerical Examples
We present three numerical examples using TKCWM in this section. Matlab R2021a is used for the simulations.

Example 1
First, let us look at the system of FDE given below: The initial values are u(0) = 0, v(0) = 1 and the exact solution for α = β = 1 is given as As a result of using the proposed approach for fractional derivatives, we now obtain, where, R T m = [r 1 , r 2 , . . . , r m ] and S T m = [s 1 , s 2 , . . . , s m ] are the unknown coefficients of size 1 × m . Using the Caputo definition for fractional derivatives [28] and using (14), (17), (22), and the initial conditions we obtain, where are also the vectors of size 1 × m . Finally, using (22) and (23) in (21) the following system of algebraic equations, which yields: Here, the R T m and S T m vectors include the total of 2m unknowns. Equation (24) is used to create a system of algebraic equations using collocation points and R T m and S T m are calculated by solving this system. Thus, the numerical solutions for u(t) and v(t) are also obtained as shown in (23).
The absolute errors obtained for TKCWM in u(t) and v(t) for several m values (α = β = 1) are given in Table 1 where the absolute errors in u(t) and v(t) are represented by E u and E v , respectively. The numerical solution results for fractional orders α = 0.75, 0.85, 0.95, β = 0.75, 0.85, 0.95 are given in Table 2. Figure 1 shows the TKCWM result graphs u(t) and v(t) for α = β = 1 with the exact solution, whereas Figure 2 displays the TKCWM result graphs u(t) and v(t) for the fractional orders α = 0.75, 0.85, 0.95, Axioms 2023, 12, 546 5 of 12 β = 0.75, 0.85, 0.95 with the integer orders α = β = 1. The proposed method substantially approximates the exact solution for the integer-orders of α = β = 1, as shown in Table 1 and Figure 1. As m increases, the absolute errors decrease. The calculated absolute errors are about 10 −4 for m = 32, 10 −4 -10 −5 for m = 64, and 10 −5 -10 −6 for m = 128. The fractional orders α, β can take any arbitrary value on the interval [0, 1], they are chosen to be equal in Table 2 and Figure 2. As Table 2 and Figure 2 demonstrate, the solution approaches to the exact solution when α, β approach 1. We provide a comparison using the Differential Transformation Method (DTM) [29] and Homotopy Perturbation Method (HPM) [30] for the fractional orders α = 0.7, β = 0.9 in Table 3. Table demonstrates that the approximate solutions for these fractional orders are close to one another for all the methods. Tables 1-3 and Figures 1 and 2 thus show that the suggested approach is a good approximation to the relevant system of FDE.

Example 2
Now consider the system of FDE given as with the initial values ( ) u 0 =0 , ( ) v 0 =0 ,and the exact solution for α=β=1 is given Employing the TKCWM method for this example, we obtain,

Example 2
Now consider the system of FDE given as with the initial values u(0)= 0, v(0)= 0, and the exact solution for α = β = 1 is given as Employing the TKCWM method for this example, we obtain, where, R T m = [r 1 , r 2 , . . . , r m ] and S T m = [s 1 , s 2 , . . . , s m ] are again the coefficients to be determined in the method. With the help of Caputo definition for fractional derivatives [28] and using (14), (17), and (26), together with the initial conditions, we obtain Finally, substituting (26)- (28) in (25) results the following system of algebraic equations where R T m and S T m coefficients are the 2m unknowns to be determined: As before, the approximate solution of the proposed method is obtained by solving (29) for R T m , S T m and substituting those coefficients in (27). The absolute errors obtained for TKCWM for u(t) and v(t) for several m values (α = β = 1) are given in Table 4 where the absolute errors in u(t) and v(t) are represented by E u and E v , respectively. As the table demonstrates, to increase the accuracy, m should be increased. Table 5 demonstrates the results of the numerical method for fractional orders for several resolutions of m . As can be seen from the table, the results approach the exact solution when fractional orders approach the integer values, which validates the results. Figure 3 shows the TKCWM result graphs u(t) and v(t) for α = β = 1 with the exact solution, whereas Figure 4 displays the TKCWM result graphs u(t) and v(t) for the fractional orders α = 0.45, 0.65, 0.85, β = 0.55, 0.75, 0.95 with the integer orders. As Figures 3 and 4 demonstrate, the proposed TKCWM method substantially approximates the exact solution for the integer-orders and the approximate solution approaches to the exact solution of the α, β approach 1 for fractional orders.

Example 3
For our last example, we consider the following system of FDEs Here, there are a total of 3m coefficients to be determined, namely T T T

Example 3
For our last example, we consider the following system of FDEs with the initial conditions u(0) = 0, v(0) = 1, w(0) = 0 and the exact solution for α = β = γ = 1 is given as u ex (t) = e t , v ex (t) = e 2t , and w ex (t) = e 3t − 1, t ∈ [0, 1]. Let us now employ TKCWM for this example and obtain Here, there are a total of 3m coefficients to be determined, namely R T Employing (31)-(32) in (30) results in the following system of algebraic equations: As before, the approximate solution is obtained by solving (33) for R T m , S T m and T T m and substituting the coefficients in (32).
The absolute errors obtained for TKCWM in u(t), v(t), and w(t) for m = 64 and m = 128 (α = β = γ = 1) are summarized in Table 6 where E u , E v , and E w denote the absolute errors in u(t), v(t), and w(t), respectively. We obtain the absolute errors around 10 −3 -10 −4 for m = 64, and around 10 −4 -10 −5 for m = 128. Again, for higher accuracy, m should be increased. The solutions for fractional orders of ff , β, γ for several resolutions is given in Table 7. When the orders α, β, γ approach 1, the TKCWM solution approaches the exact solution for integer orders, which is expected. Figure 5 shows the TKCWM result graphs of u(t), v(t), and w(t) for α = β = γ = 1, whereas Figure 6 displays the TKCWM result graphs of u(t), v(t), and w(t) for the fractional orders 0.5, 0.7, 0.9 with α = β = γ = 1. As is the case with the previous examples, the solution approaches the exact solution when α, β, γ approach 1.

Conclusions
We have developed a numerical solution approach for the systems of FDEs in this study. The operational matrices for the fractional integration are obtained using the suggested method, which makes use of discrete third kind Chebyshev Wavelets. We employ the block pulse functions and Chebyshev polynomials of the third kind in the method, which makes the operational matrices for fractional integration very sparse, which in turn is crucial for fast evaluation and reduced computational cost for the proposed method. Those sparse operational matrices are used to map fractional terms in FDEs into the discrete terms in the proposed method. The Newton-Raphson method is used to solve the algebraic equation system created from the system of FDEs in order to determine the unknown coefficients. As predicted, the accuracy of the suggested technique improves as the resolution is increased by larger collocation points. Tables 1, 4 and 6 and Figures 1, 3, and 5 verify this claim. The maximum errors for the collocation points up to 128 are typically in the range of 10 −4 to 10 −6 .
The main focus of this study is fractional orders, and the numerical solutions derived for these orders are consistent with those found for integer orders. Tables 2, 3, 5, and 7 and Figures 2, 4, and 6 show how the method accurately solves fractional orders because when the fractional orders become close to integer values, the answer likewise become close to the exact solution found for integer orders.
Systems of variable-order FDEs, fractional partial differential equations, fractional integral equations, and fractional delay differential equation systems can all benefit from the advantages of the proposed method and are thus taken into consideration as potential future applications of this research.