Numerical Solution of Variable-Order Fractional Differential Equations Using Bernoulli Polynomials

We introduce a new numerical method, based on Bernoulli polynomials, for solving multiterm variable-order fractional differential equations. The variable-order fractional derivative was considered in the Caputo sense, while the Riemann--Liouville integral operator was used to give approximations for the unknown function and its variable-order derivatives. An operational matrix of variable-order fractional integration was introduced for the Bernoulli functions. By assuming that the solution of the problem is sufficiently smooth, we approximated a given order of its derivative using Bernoulli polynomials. Then, we used the introduced operational matrix to find some approximations for the unknown function and its derivatives. Using these approximations and some collocation points, the problem was reduced to the solution of a system of nonlinear algebraic equations. An error estimate is given for the approximate solution obtained by the proposed method. Finally, five illustrative examples were considered to demonstrate the applicability and high accuracy of the proposed technique, comparing our results with the ones obtained by existing methods in the literature and making clear the novelty of the work. The numerical results showed that the new method is efficient, giving high-accuracy approximate solutions even with a small number of basis functions and when the solution to the problem is not infinitely differentiable, providing better results and a smaller number of basis functions when compared to state-of-the-art methods.


Introduction
In the last few decades, fractional calculus has attracted the attention of many scientists in different fields such as mathematics, physics, chemistry, and engineering. Due to the fact that fractional operators consider the evolution of the system, by taking the global correlation, and not only local characteristics, some physical phenomena are better described by fractional derivatives [1]. There are many definitions of fractional differentiation and integration in the literature (for details, see, e.g., [2][3][4][5]). However, the most commonly used are the definitions of the Caputo derivative and Riemann-Liouville fractional integral operators. For new fractional derivatives with nonlocal and nonsingular kernels, with applications in rheological models, we refer the reader to [6].
A recent generalization of the theory of fractional calculus is to allow the fractional order of the derivatives to be dependent on time, i.e., to be nonconstant or of variable order. In [7], arXiv:2111.07413v1 [math.NA] 14 Nov 2021 the authors investigated operators when the order of the fractional derivative is variable with time [8,9].
The nonlocal properties of systems are more visible with variable-order fractional calculus, and many real-world phenomena in physics, mechanics, control, and signal processing have been described by this approach [10][11][12][13]. In particular, several applications of variableorder fractional calculus are found in engineering mechanics; see [14] for an application of variable-order fractional operators to model the microscopic structure of a material, [15] for an application of the Riesz-Caputo fractional derivative of space-dependent order in continuum elasticity, [16,17] for the nonlinear viscoelastic behavior of fractional systems with variable time-dependent fractional order, and [18] for the use of variable-order fractional calculus in the static response of nonlocal beams having either a porous or a functionally graded core.
Obtaining analytic solutions for such fractional differential equations (FDEs) is, however, very difficult. Therefore, in most cases, the exact solution is not known, and one needs to seek a numerical approximation. Therefore, many researchers have introduced and developed numerical methods in order to obtain approximated solutions for this class of equations. For example, in [19,20], Legendre polynomials were used to construct a numerical solution for a class of multiterm variable-order FDEs, while [21] used Legendre wavelets and operational matrices; References [22,23] used Chebyshev polynomials, respectively of the second and fourth kinds; Reference [24] used Bernstein polynomials. The book [25] showed the usefulness of numerical methods for approximating variable-order fractional operators in the framework of the calculus of variations; the paper [26] adopted Coimbra's variable-order time-fractional operator, discussing the stability, convergence and solvability of a numerical scheme based on Fourier analysis. Reference [27] implemented a numerical method for solving a circulant Halvorsen system described by Caputo fractional variable-order derivatives. For multivariable-order differential equations with nonlocal and nonsingular kernels, we refer to [28], where a collocation method was developed based on Chebyshev polynomials of the fifth kind.
Here, we considered the following general form of a multiterm variable-order FDE: with initial conditions where n is the smallest positive integer number such that for all t ∈ [0, 1], one has 0 < α(t) ≤ n, are the (left) fractional derivatives of variable-order defined in the Caputo sense. Problems of the form (1) and (2) have a practical impact. Indeed, specific applications are found in noise reduction and signal processing [29,30], the processing of geographical data [31], and signature verification [32].
In recent years, Bernoulli polynomials have been shown to be a powerful mathematical tool in dealing with various problems of a dynamical nature, e.g., for solving numerically high-order Fredholm integrodifferential equations [33], pantograph equations [34], partial differential equations [35], linear Volterra and nonlinear Volterra-Fredholm-Hammerstein integral equations [36], as well as optimal control problems [37]. Here, we employ a spectral method based on Bernoulli polynomials in order to obtain numerical solutions to the problem (1) and (2). Our method consists of reducing the problem to a system of nonlinear algebraic equations. To do this, we introduced an accurate operational matrix of variable-order fractional integration for the Bernoulli polynomials' basis vector. To the best of our knowledge, this is the first time in the literature that such a method for solving a general class of multiterm variable-order FDEs based on the Riemann-Liouville fractional integral of the basis vector has been introduced.
The rest of this paper is organized as follows. In Section 2, some preliminaries of variableorder fractional calculus are reviewed and some properties of the Bernoulli polynomials are recalled. Section 3 is devoted to introducing the operational matrix of variable-order fractional integration for Bernoulli polynomials. In Section 4, we present a new numerical method for solving the problem (1) and (2) by using the operational matrix technique and collocation points. Section 5 is concerned with presenting an error estimate for the numerical solution obtained by this new scheme. In Section 6, several multiterm variable-order FDEs are considered and solved, using the introduced method. Finally, concluding remarks are given in Section 7, where some possible future directions of research are also pointed out.

Preliminaries
In this section, a brief review of the necessary definitions and properties of the variableorder fractional calculus is presented. Furthermore, Bernoulli polynomials, and some of their properties, are recalled.

Some Preliminaries of Variable-Order Fractional Calculus
We followed the notations of [1].

Bernoulli Polynomials
Bernoulli polynomials build a family of independent polynomials that form a complete basis for the space L 2 [0, 1], which is the space of all square integrable functions on the interval [0, 1]. The Bernoulli polynomial of degree m, β m (t), is defined as follows [38]: where b k , k = 0, 1, . . . , m, are the Bernoulli numbers that appear in the series expansion of trigonometric functions [39] and can be defined by the following identity: The first four Bernoulli polynomials are: The Bernoulli polynomials satisfy the following property [39]: Any arbitrary function y ∈ L 2 [0, 1] can be approximated using the Bernoulli polynomials as where and The coefficient vector A in (7) is calculated by the following formula (see [37]): where ·, · is the inner product, defined for two arbitrary functions f , g ∈ L 2 [0, 1] as and D = B(t), B(t) , which is calculated using (6).

Operational Matrix of Variable-Order Fractional Integration
The aim of this section is to introduce an accurate operational matrix of variable-order fractional integration for Bernoulli functions. To do this, we display the Bernoulli basis vector B(t), given by (8), in terms of the Taylor basis functions, as follows: where T is the Taylor basis vector and Q is the change-of-basis matrix, which is obtained using (5) as It is obvious that the matrix Q is a nonsingular matrix. Therefore, we can write: Taking (10) into account and by applying the left Riemann-Liouville fractional integral operator of order α(t) to the vector B(t), we obtain where S is a diagonal matrix, which is obtained using Lemma 1 as follows: Finally, by using (11) into (12), we have: where P , which we call the operational matrix of variable-order fractional integration of order α(t) of Bernoulli functions.
Since Q and Q −1 are lower triangular matrices and S is also a lower triangular matrix. For example, with M = 2, one has: where

Numerical Method
This section is devoted to presenting a new numerical method for solving the multiterm variable-order FDE (1) with initial conditions (2). To this aim, we set: Then, by assuming y ∈ C n [0, 1], we considered an approximation of the n-th order derivative of the unknown function y using the Bernoulli functions as follows: where A is an (M + 1) × 1 vector with unknown elements and B(t) is the Bernoulli basis vector given by (8). Taking into account the initial conditions (2) and using (3), (13) and (14), we obtain: In a similar way, using (4), (13), and (14), we obtain and By substituting the approximations given in (15)- (17) into Equation (1), we have: . (18) By using (18) at M + 1 collocation points, which are chosen as t j = j+1 M+2 , with j = 0, 1, . . . , M, we obtain the following system of nonlinear algebraic equations: System (19) includes M + 1 nonlinear algebraic equations in terms of the unknown parameters of vector A. After solving this system, an approximation of the solution of the problem (1) and (2) is given by (15).

Error Estimate
The purpose of this section is to obtain an estimate of the error norm for the approximate solution obtained by the proposed method in Section 4. We assume that f (t) = y (n) (t) is a sufficiently smooth function on [0, 1] and q M (t) is the interpolating polynomial to f at points t s , where t s , s = 0, 1, . . . , M, are the roots of the (M + 1)-degree shifted Chebyshev polynomial in [0, 1]. Then, according to the interpolation error, we have: Therefore, We assume that there is a real number κ such that: By using (21) in (20) and taking into consideration the estimates for Chebyshev interpolation nodes [40], we obtain: As a consequence of (22), we obtain the following result.  (8) and (9). Then, there exists a real number κ such that: Proof. Let Π M be the space of all polynomials of degree ≤ M on t ∈ [0, 1]. By definition, f M is the best approximation of f in Π M . Therefore, we have: where g is any arbitrary polynomial in Π M . Therefore, we can write: where q M is the interpolating polynomial of f , as discussed before. As a result, taking into consideration (22) in (24), we easily obtain (23).
With the help of Theorem 1 and by assuming f (t) = y (n) (t), we obtain the following result.

Theorem 2.
Let the exact solution y of the problem (1) and (2) be a real-valued sufficiently smooth function and y M be the approximate solution of this problem obtained by the method proposed in Section 4. Then, we have: where y is the exact solution and κ = max t∈(0,1) y (nM+n) (t) , n = max 0<t≤1 { α(t) }.
Proof. Let X be the space of all real-valued functions defined on [0, 1] and 0 I n t : X → X be the Riemann-Liouville integral operator. We use the following definition for the norm of the operator 0 I n t : 0 I n t 2 = sup To continue the proof, first, we introduce an upper bound for 0 I n t 2 . To this end, using the definition of the left Riemann-Liouville integral operator and Schwarz's inequality, we obtain: Therefore, we have: On the other hand, from Theorem 1, we have the following error bound: Finally, using (3), (15), (25) and (26), we obtain that: which completes the proof.

Illustrative Examples
In this section, we apply our method to some variable-order FDEs and, moreover, to one variable-order fractional pantograph differential equation (Example 4), comparing the results with the ones obtained by existing methods in the literature. We implemented our method and performed our numerical simulations with Mathematica 12. Example 1. In our first example, we considered the following multiterm variable-order FDE [19,22,23]: where with initial conditions y(0) = 2 and y (0) = 0. The exact solution to this problem is y(t) = 2 − t 2 2 . As can be seen, we have α(t) = 2t. Therefore, to implement the proposed method, we introduce: We set M = 1 and suppose The operational matrices of variable-order fractional integration are given as follows: Furthermore, we have: . Now, using the initial conditions and the aforementioned operational matrices, we obtain the following approximations for y(t) and its variable-order derivatives: .
By substituting these approximations into (27), collocating the resulting equation at t 0 = 1 3 , t 1 = 2 3 , and finally, solving the obtained system, we obtain: Therefore, we have: which is the exact solution. In this case, since the exact solution is a second-order polynomial, we can obtain it by applying the numerical method with just two basis functions.

Example 2.
In our second example, we considered the following nonlinear variable-order FDE borrowed from [41]: The exact solution of this problem is y(t) = t 7 2 . By considering α(t) = 1 − 0.5 exp(−t), we solved this problem with different values of M. The numerical results are displayed in Figure 1 and Table 1. In Figure 1, the approximate solutions obtained with M = 1, 2, 3, together with the exact solution of this problem, are plotted. Furthermore, by considering M = 2, 6, 10, the absolute errors at some selected points are reported in Table 1. From these results, the convergence of the numerical solutions to the exact one can be easily seen.  y(t) + y(t) + e t y(t 5 ) = g(t), 0 < t ≤ 1, where The exact solution of this problem is y(t) = t 3 + t 2 . To solve Problem (28), we applied the method with M = 1 and M = 2. The numerical solution obtained with M = 1, together with the exact solution are plotted in Figure 2. With M = 2, according to the method described in Section 4, we set n = max 0<t≤1 { sin(t) } = 1. Therefore, by assuming that and using the initial condition, we have: By substituting these approximations into (28) and using the collocation points t 0 = 1 4 , t 1 = 1 2 , and t 2 = 3 4 , we obtain a system of three nonlinear algebraic equations in terms of the elements of the vector A. By solving the resulting system, one obtains: Finally, using these values, we obtain: which is the exact solution.

Example 4.
Let us now consider the fractional pantograph differential equation: The exact solution of this problem, when α(t) = 1, is y(t) = e −t . By choosing M = 1 and α(t) = 1, we set: Then, by considering the initial condition, we have: By substituting these approximations into (29), and using the collocation points we obtain a 0 = −0.620328, a 1 = 0.621053, which gives This approximate solution and the exact solution to the problem, corresponding to α = 1, are displayed in Figure 3. By computing the L 2 -norm of the error for this approximation, we have: which shows that the method gives a high-accuracy approximate solution, even with a small number of basis functions. A comparison of the absolute errors, obtained by the proposed method with M = 6, 8, 10 at some selected points, with the results proposed in [43], using modified hat functions, and those of [44], using Bernoulli wavelets, are reported in Table 2. From this table, it is seen that our method gave more accurate results with a smaller number of basis functions when compared to previous methods. Moreover, the approximate solutions obtained with M = 2 and different α(t), along with the exact solution of corresponding first-order equation, are given in Figure 4. This figure shows that the numerical solution is close to the exact solution for the case α(t) = 1 when α(t) is close to one.  2.00 × 10 −8 9.30 × 10 −9 1.59 × 10 −11 2.42 × 10 −13 2 −5 5.34 × 10 −10 3.70 × 10 −9 6.47 × 10 −9 1.21 × 10 −11 1.29 × 10 −13 2 −6 2.27 × 10 −9 2.03 × 10 −8 3.83 × 10 −9 7.58 × 10 −12 6.72 × 10 −14 Exact  Example 5. As our last example, consider the following variable-order FDE found in [45]: where α(t) = 0.25(1 + cos 2 (t)) and the initial condition is y(0) = 1. The exact solution is y(t) = e t . We solved this problem using the proposed method in this paper with different values of M. A comparison of the absolute errors at some selected points, obtained by our method and the one of [45], based on a class of Lagrange polynomials, is given in Table 3.

Concluding Remarks
Many researchers have employed fractional differential equations (FDEs) in order to model and analyze various scientific phenomena. Typically, such FDEs do not have known analytical solutions, and approximate and numerical approaches have to be applied [46]. Here, a new numerical method, based on Bernoulli polynomials, was presented for solving multiterm variable-order fractional differential equations. The operational matrix of variableorder fractional integration for the Bernoulli basis functions was introduced, which is a lower triangular matrix and helps to reduce the computational effort of the method. Our scheme uses this matrix to give some approximations of the unknown solution of the problem and its variable-order fractional derivatives in terms of the Bernoulli functions. Substituting these approximations into the equation and using some collocation points allowed us to reduce the problem to a system of nonlinear algebraic equations, which greatly simplifies the problem. An error estimate of the method was proven, and the applicability of our method was illustrated by solving five illustrative examples. The obtained results confirmed the efficiency, accuracy, and high performance of our technique, when compared with the state-of-the-art numerical schemes available in the literature. We emphasize that the accuracy of our method is preserved even when the solution of the problem is not infinitely differentiable. This can be observed in Example 2, whose exact solution is t 7 2 . We used a variable-order definition where the operator has no order memory, a so-called type-I operator. As future work, it would be interesting to extend the proposed numerical method to approximate variable-order derivatives with weak (type-II) and strong (type-III) variable-order definitions [18,47]. Other interesting lines of research include the stability analysis of the proposed numerical method and its application to different areas in science and engineering, for example in structural mechanics.