Effective Method for Solving Higher-Order Fractional Differential Equations Using Spline Functions

: In this paper, appropriate spline functions in polynomial form are outlined and applied to solve higher-order linear fractional differential equations (H-OLFDEs). Several techniques have been proposed to solve this type of equations. A description of the projected methodology is first introduced. The method involves transforming the H-OLFDE of order α , with given initial conditions, into a system of fractional differential equations with 0 < α i < 2 , denoting the order of the Caputo fractional derivative for each equation i in the transformed system, where the number of equations in the resulting system is equal to that of initial conditions. The numerical results obtained using this strategy demonstrate substantial agreement with the precise analytical solutions. The results in this paper affirm the robustness and ease of application of the proposed technique.


Introduction
Differential equations with fractional order are generalizations of typical differential equations with integer order.These types of fractional differential equations (FDEs) are increasingly used to represent problems in viscoelasticity, fluid flow, mechanics, biology, physics, engineering, and many other applications [1,2].Finding analytical and numerical methods that can be relied on, and efficient to obtain accurate solutions to FDEs has been an interest of scientific researchers.Since exact solutions for the fractional differential equations cannot be found simply, and sometimes not possible, then numerical solutions are the available options using some analytical and reliable methods.
Modeling complex processes mathematically is a great challenge for scientists.On the other hand, for basic classical systems, the theory of differential equations with integer order is sufficient to describe the dynamics of these systems.Fractional derivatives are important since they provide a crucial tool for the outline of memory and innate properties within diverse complex systems and materials.It is noticed that, in the last few decades, the amount of scientific and engineering problems represented by fractional derivatives equations has surged and maybe the fractional differentiation and integration could be the calculus of the current century [3][4][5][6][7][8].
Investigations of the analytic aspects concerning the existence, and the uniqueness of the solutions to the FDEs have been done by several researchers, as exemplified in the references [9,10].Several analytical and numerical strategies to solve FDEs have been proposed throughout the last few decades.The foremost commonly used techniques are fractional difference method [11,12], decomposition method [13], Adomian decomposition method (ADM) [1,14], variational iteration method (VIM) [15,16], and the Adams-Bashforth-Moulton method (ABMM) [17][18][19], Chelyshkov matrix method [20].Lately, the homotopy perturbation method (HPM) [21], and the Lagrange multiplier method have been applied to numerically solve multi-order fractional differential equations [22].A higher-order numerical technique was used to solve a nonlocal neutron diffusion differential equation incorporating delayed neutrons that represents neutron transport in a nuclear reactor, the used procedure was based on utilizing an approximation of the temporal derivative L1-2 methodology, together with a space discretization through the compact difference method [23].In [24], the authors addressed the numerical solution for the 3-dimensional non-local evolution equation featuring a weakly singular kernel, their method was based on constructing the backward Euler alternating direction implicit numerical algorithms.Orthogonal spline collocation technique was used to solve the 4th order fractional sub-diffusion equation [25,26].New operational matrices for fractional derivatives of a first class of classical orthogonal polynomials were used to solve fractional differential equations [27].
The method of spline functions (SFs) in polynomial form has emerged as a powerful approach to solving FDEs, providing a flexible and efficient numerical framework for approximating solutions.SFs offer a natural way to represent functions over intervals.One of the key motivations behind using SFs lies in their ability to efficiently capture the behavior of complex systems governed by fractional calculus.By expressing the solution of an FDE as a spline function, researchers can discretize the problem and apply standard numerical techniques to obtain approximate solutions.
By considering the successful applications of the SFs in polynomial form method in various areas such as system analysis [28], delay differential equations [29], linear and nonlinear FDEs [30], systems of linear and nonlinear FDEs [31], we have a tendency to expect that the spline functions in polynomial form method can be applicable to address linear higher-order fractional differential equations.
In the current work, our intention is to expand the application of the SFs in polynomial form in order to address the higher-order fractional differential equations of the general form: with the corresponding initial conditions where and D γ * y(x) represents the fractional derivative of order γ > 0 in Caputo sense.
We assume that α > Numerous types of visco-elastic damping problems have been modeled using higher-order fractional differential equations [32].Notably, the proposed model equations have predominantly taken on a linear form.So, our numerical tests in this work will concentrate on equations of linear nature subject to the initial conditions in equation (2).
The proposed method is simple and straightforward; the procedure is broken down into three steps, which are discussed at the end of Section 3., and it may be implemented in any programming language.We shall show that the numerical results are substantially consistent with the actual solutions.
The paper is structured as follows: in Section 2. , we present some well-known definitions and fundamental theorems used in the current work.In Section 3. , we present some methods used to solve the H-OLFDEs, and outline our proposed technique.In Section 4. , we consider three examples to show the effectiveness of the discussed method, accompanied by data and figures.Section 5 is the conclusion.

Preliminaries and definitions
Here, some well-known definitions and theorems [7,11,30,33,34], that are utilized in this work; Riemann-Lioville for fractional integral, Caputo and Riemann-Lioville for fractional derivatives, also, the beta function and relation with gamma function.
Definition 1 Let α be a positive real number.The operator J α a , defined on L 1 [a, z] by is called the Riemann-Lioville (RL) fractional integral operator of order α.
We set J 0 a : = I, in the case when α = 0, which is the identity operator.Definition 2 The Caputo fractional derivative c D α a of a given function f , of order α > 0, with a nonnegative real number a, is defined by where m is a natural number, x ≥ a, and f (x) ∈ C m −1 .The applicability of RL fractional operators can be expanded to accommodate for larger values of α.To express this, we set α = ⌊α⌋ + β , where ⌊α⌋ stands for the integer part of α, and β = α − ⌊α⌋, thus we give the next definition.
Definition 3 If α > 0, and α is not a natural number, then define In the case when α is negative, the definition becomes

Definition 4
The Beta function, defined by means of integral, is defined as The relationship between the beta and the gamma functions is given by the following formula: , where Re u > 0, and Re v > 0.
The integral in Definition ( 4) is defined on the interval [0, 1], for a general interval [a, b] we obtain the following formula, Notice that, in the special case when a = 0, and z = 1, we obtain the formula in Definition (4).
Definition 5 The incomplete Beta function B(x; u, v) is defined by where Re u > 0, Re v > 0, and 0 < x < 1.When x = 1, the incomplete Beta function reduces to the usual Beta function, i.e., B(1; u, v) = B(u, v).Definition 6 The regularized incomplete Beta function denoted by I(x; u, v) is defined in terms of the complete and incomplete Beta function as:

Spline function method of polynomial form for higher-order differential equations
Considerable research efforts have been devoted to approximating solutions to systems of ordinary and delay differential equations through the application of spline functions as indicated in studies [28,29].Specifically, Micula et al. [28] considered a system of two equations with their initial solutions: where the functions In their work, they required that each of the functions g (α) i , i = 1, 2, and for α = 0, 1, . . ., r must satisfy the condition of Lipschitz, i.e., |g (α) These required conditions guarantee the existence of the unique solution of the system (4).In their research, Micula et al. explored expanding the form of SFs to approximate the solution of the system of initial value problems (IVPs) outlined in (4), aiming for unique solutions w = w(t), and v = v(t).

Volume 5 Issue 3|2024| 3349 Contemporary Mathematics
Al-Rabtah et al. [30] applied a spline function method in polynomial form for both linear and nonlinear FDEs.The formulation of the fractional spline function S(x) in polynomial form employed to approximate the actual solution was defined as where , and S −1 (x 0 ) = y 0 , for x in the subinterval [x k , x k+1 ], which was used to solve the IVP with Ω representing the uniform partitioning of the interval [a, b], defined by the nodes Ω: The formulation of the discussed method was grounded in the fractional generalization of the Taylor's theorem [35], which is defined in the next theorem: Theorem 1 Let n be a positive integer, α > 0, and the function is the remainder.
In [31], Al-Rabtah extended the utilization of the SFs method in polynomial form to address systems of linear and nonlinear FDEs, also, the applicability of the method was analyzed in the case 0 where m represents the number of equations in the system, namely, where a ≤ x ≤ b, and D α i is the derivative of y i of order α i (0 < α i ≤ 1) in the notion of Caputo, and f i 's are arbitrary linear or nonlinear functions, subject to the initial conditions In this paper, we expand the application of the SFs in polynomial form to address higher-order fractional differential equations of the form (1), considering to the initial conditions (2).The work in [31] was limited to a specific value for each α i in the interval (0, 1], however, in this research, we choose values in the interval (0, 2).
To solve this higher-order fractional differential equation, we first transform the equation into a system of fractional differential equations with 0 < α i < 2, where α i represents the order of the Caputo fractional derivative for each equation i in the system with the corresponding initial conditions Let Ω represents the uniform partitioning to the given interval [a, b], defined by the nodes Ω: Secondly, define the form of the fractional spline function S(x) in polynomial form to approximate the actual solution y by: Thirdly, evaluate the previous system at each value of k in the uniform partitioning.All calculations are implemented with MATLAB to obtain the numerical results.

Numerical examples
To illustrate the effectiveness of this proposed method, we examine three specific examples.These examples are chosen since their analytical solutions are available.This allows for a direct comparison between the results obtained using this method with the actual solution.The calculations were performed on a laptop running Matlab R2021a with an Intel(R) Core(TM) i5-1135G7 CPU at 2.40 GHz, 8.00 GB Ram, 4 cores, and 8 logical processors.
In addition to comparisons with the actual solutions, the numerical results are compared to those obtained using the fractional Euler's method.The general formula for the fractional Euler's method is given by, see [36], , Example 1 [37] Consider the higher-order linear fractional differential equation subject to the initial conditions (ICs), The actual solution of this higher-order initial value problem is E α (−x α ), where x k /Γ(αk + 1), which is known as the Mittag leffler function.In order to solve this present problem, we transform the higher-order equation into a system of two FDEs, the first equation is of order α 1 , and the second one is of order α 2 , i.e.
with the corresponding ICs,

Contemporary Mathematics
Table 1 shows the numerical approximate solution compared to the actual solution in the case where α = 1.8, α 1 = 0.95 for the first equation, α 2 = 0.85 for the second one, and step size h = 0.1, the maximum relative error in this case is 2.1 × 10 −4 .The computing time for the solution of the resulting system is 3.1250 × 10 −2 seconds.Figure 1 shows the errors.
Table 1 also shows the numerical approximate solution obtained by the fractional Euler method compared to the actual solution, the maximum relative error in this case is 7.3 × 10 −2 .Additionally, comparisons between the solutions produced by the fractional Euler technique and the SFs method at the same parameters are displayed in the table.It is clear that the results obtained using the SFs method are much better than those obtained using the fractional Euler method.
Table 1.The solution of Example (1) where α = 1.8, α 1 = 0.95, α 2 = 0.85, and h = 0.1 is the step size.In the case of α = 2, equation ( 5) is an ordinary differential equation, and the actual solution is y(x) = cos(x).We also solve this equation as described earlier.Table 2 shows the numerical solution compared to the actual solution in the case where α = 2, α 1 = 0.95 for the first equation, α 2 = 1.05 for the second one, and step size h = 0.1, the maximum relative error in this case is 1.4528 × 10 −3 .The computing time for the solution of the resulting system is 1.5625 × 10 −2 seconds.
Figure 2 shows the numerical approximate solution curves generated by the SFs method and the fractional Euler method compared to the actual solution curve in the case where α = 1.9, α 1 = 0.9 for the first equation, α 2 = 1 for the second one, and step size h = 0.02, the maximum relative error in this case using the SFs method is 3.5536 × 10 −3 , and 2.4075 × 10 −1 when using the fractional Euler method.The computing time for the solution of the resulting system using the SFs method is 3.1250 × 10 −2 seconds.It is evident that the outcomes of the SFs approach are significantly superior to the fractional Euler method.The errors for the approximate solution produced using the SFs approach are displayed in Figure 3.  2) where α = 1.9, α 1 = 0.9, α 2 = 1 , and step size h = 0.02 Example 3 Consider the higher-order linear fractional differential equation Volume 5 Issue 3|2024| 3355 Contemporary Mathematics subject to the ICs, The actual solution of this IVP is y(x) = x 2 , see [17].To solve this problem, we transform the higher-order equation into a system of two FDEs, the first equation is of order α 1 , and the second one is of order α 2 , i.e.
subject to the corresponding ICs, Table 3 shows the numerical solution compared to the exact solution in the case where α = 1.75, α 1 = 0.85 for the first equation, α 2 = 0.9 for the second one, and step size h = 0.1, the maximum relative error in this case is 3.2044 × 10 −3 .The computing time for the solution of the resulting system is 3.1250 × 10 −2 seconds.
Table 3 additionally presents the numerical approximate solution obtained by the fractional Euler method compared to the actual solution, the maximum relative error in this case is 1.0024×10 −1 .The table also shows comparisons between the solutions obtained by the SFs method and the fractional Euler method at the same parameters.It is obvious that the results obtained using the SFs method are much better than those obtained using the fractional Euler method.Figure 4 shows the numerical solution curve compared to the exact solution curve in the case where α = 1.9, α 1 = 0.95 for the first equation, α 2 = 0.95 for the second one, and step size h = 0.05, the maximum relative error in this case is 1.7964 × 10 −3 .The computing time for the solution of the resulting system is 3.1250 × 10 −2 seconds.Figure 5 shows the Log-errors to demonstrate the stability of the solution.

Conclusions
In this work, we expanded and generalized the strategy of the polynomial form of spline functions, to solve higherorder linear fractional differential equations by transforming the H-OLFDE of order α, with given initial conditions, into a system of fractional differential equations, where the order of each equation is 0 < α i < 2, with α i represents the order of the fractional derivative in the Caputo notion for every equation i in the system, and the number of equations in the resulting system is equal to the number of initial conditions.The obtained approximate numerical results demonstrate substantial agreement with the actual solutions.
and the step size h = (b − a)/n.
and the step size h = (b − a)/n.
and the step size h = (b − a)/n, which is used to solve the IVP,