Green–Haar wavelets method for generalized fractional differential equations

The objective of this paper is to present two numerical techniques for solving generalized fractional differential equations. We develop Haar wavelets operational matrices to approximate the solution of generalized Caputo–Katugampola fractional differential equations. Moreover, we introduce Green–Haar approach for a family of generalized fractional boundary value problems and compare the method with the classical Haar wavelets technique. In the context of error analysis, an upper bound for error is established to show the convergence of the method. Results of numerical experiments have been documented in a tabular and graphical format to elaborate the accuracy and efficiency of addressed methods. Further, we conclude that accuracy-wise Green–Haar approach is better than the conventional Haar wavelets approach as it takes less computational time compared to the Haar wavelet method.


Introduction
The mathematical theory of fractional calculus can be described as that of derivatives as well as integrals of any order possible. Primarily, fractional calculus is a generalized form of integer order calculus. Fractional calculus has been exploited as a crucial tool for applications that concern science and engineering. These applications of fractional calculus have been elaborated previously by several authors. In essence, fractional calculus has been deployed for modeling the transfer of heat in heterogeneous media [1], nonlinear oscillation of earthquakes [2], signal processing [3], neural networks [4][5][6], fluid dynamic traffic flow [7], electromagnetism [8], bioengineering [9], economics [10], anomalous diffusions and fractal-like nature [11,12]. For the qualitative analysis of fractional differential equations, we refer the reader to [1,2,[13][14][15] and the references therein.
Several books have been written on the philosophy and development of fractional calculus [16][17][18][19]. In fractional calculus the fractional derivative is introduced via fractional integral. Riemann, Liouville, Caputo, Hadamard, Grunwald and Letinkow are the pioneering researchers who have been contributing and publishing extensively about these applications. Meanwhile, the literature has witnessed the appearance of different types of fractional derivatives that improve and generalize the classical fractional operators defined by the above listed authors [20][21][22]. Recently, Katugampola rediscovered a new type of fractional integral operator which covers both Riemann-Liouville and Hadamard operators and represents them in a single form [23,24].
The study of wavelet theory dates back to mid-20th century. Once it had been introduced, the theory has had prominent contributions in mathematical studies [25]. It is a significant tool for science and engineering. Wavelets are being used for analyzing signals, for representation of waveform and segmentation, optimal control, numerical analysis, fast algorithm for easy implementation, and time-frequency analysis [26]. There are many kinds of wavelets, for example, Haar [27][28][29][30], Daubechies [31], B-spline [32], Battle-Lemarie [33], Legender [34], as well as Green-CAS [30]. A naive form of orthonormal wavelets which employ compact support has been used by many researchers and is called the Haar wavelet. Mathematically, Haar wavelet family consists of rectangular functions.
Further, it contains the lower member of Daubechies family of wavelets which is appropriate for computer implementations. Primarily, Haar wavelets convert a fractional differential equation into an algebraic system of equations with finite variables. The Haar wavelets approximation for tackling linear and nonlinear systems has been discussed in [35][36][37][38][39].
The prospective study has been established to solve the generalized fractional differential equations numerically. The numerical computation utilizes Haar wavelets, as well as Green-Haar approach for this purpose. Operational matrices have been developed using Haar wavelet approach. These matrices are thus employed to solve generalized fractional differential equations. The integral operator used for the purpose of computing operational matrix is the generalized Reinmann-Liouville fractional integral operator. An error analysis for convergence of the proposed technique has been undertaken through a generalized Caputo-type fractional differential operator. The method has been further elaborated in terms of efficiency and accuracy by considering a number of documented examples. A comparison of these results has also been presented against previous studies [40] to further emphasize the accuracy and efficiency of the proposed technique.
One important feature of the method is that it does not require an operational matrix at all. Stability and convergence of this method have also been derived for further applications. The undertaken study shows that the method is even more computationally efficient against the standard Haar wavelet technique discussed in the same study. Interestingly, the accuracy is not compromised, but rather enhanced by using Green-Haar method for solving generalized fractional boundary-value problems.
The present study is structured as follows: In Sect. 2, we review basic mathematical expressions of fractional calculus with their respective definitions. Furthermore, we reflect on Haar wavelets which is an essential preliminary topic for subsequent sections in this paper. In Sect. 3, we develop operational matrices using a generalized integral operator and Haar wavelets which help in estimating the numerical solution of a generalized Caputo-type fractional differential equation. Moreover, we establish an upper bound for the proposed technique through Haar wavelets for the generalized fractional differential equation. Further, numerical solutions are given to elaborate the accuracy and efficiency of the numerical scheme. We propose a new method called Green-Haar method for the boundary value problems and compare our results against the Haar wavelet approach in Sect. 4. We summarize the outcomes of this paper in the last section.

Preliminaries
In this section, for the sake of convenience, we will review some necessary definitions. These definitions serve as essential preliminaries for fractional calculus and Haar wavelets. These definitions are going to assist in upcoming sections. ([23, 24]) Consider α, ρ ∈ R + such that α > 0. The generalized fractional integral (I α,ρ a f ) (in the sense of Katugampola) is given by

Fractional calculus
Now we introduce the Caputo-type generalized fractional derivative such that, at two convenient limits, this generalized Caputo-type fractional derivative recovers the wellknown Caputo-Hadamard and Caputo fractional derivatives. ([23, 24]) Consider α, ρ ∈ R + such that α > 0 and n = α + 1. The generalized Caputo-type fractional derivative is defined by

Definition 2.2
where δ n ρ = (t 1-ρ d dt ) n and α represents the order of the fractional operators.

Haar wavelets and function approximation
The domain of Haar wavelets is an essential component of a set of those wavelets which employ compact support. The functions forming the family of Haar wavelets consists of step functions over the real line. These are the functions which are restrained to only the values -1, 0, and 1. These functions have two characteristics. Firstly, they are discontinuous in their nature. Secondly, their derivative vanishes. Each function that falls into the category of Haar wavelets is essentially defined over the interval t ∈ [a, b) except for the scaling function conveyed as [39] where, ξ 1 (i) = a + (ba) k m , ξ 2 (i) = a + (ba) 2k+1 2m , ξ 3 (i) = a + (ba) k+1 m . We define the quantity m = 2 j , where j = 0, 1, 2, 3, . . . , J and k = 0, 1, 2, 3, . . . , m -1. Here parameter j is used as a representation for the level of wavelet or dilation parameter, translation is represented by t, while the maximal level of resolution for the Haar wavelet is represented by J. The relation between the parameter m, k, and i is as i = m + k + 1. Equation (7) is valid for i ≥ 3. It is presumed that the values i = 1 and i = 2 correspond to the following scaling functions, respectively: and If u(t) is a function defined on the interval [0, 1], it should decompose as where c i = u(t), h i (t) . In particular, the first m terms are considered, such that m is a power of 2, Lemma 2.7 ([43]) Suppose that a function v(t) is differentiable and has bounded first derivative over (0, 1), that is, there exits M > 0 such that |v (t)| ≤ M for all t ∈ (0, 1), and also assume that v k (t) is an approximation of v(t), then we have

Operational matrices method
Operational matrices have been widely used to deal with fractional order systems. Several authors established Haar wavelets operational matrices to deal with various problems, such as to find the numerical solutions of linear and nonlinear initial as well as boundary value problems of fractional order [29,30,44,45]. Hsiao and Chen [46] established an operational matrix to study lumped dynamical systems with distributed-parameters. Wang and Hsiao [47] have solved an optimal control system by linearly changing through Haar wavelets. Dai and Cochran Jr. [48] have considered a Haar wavelet technique to transform an optimal control system in the direction to nonlinear programming (NLP) parameters using collocation points. This NLP can be solved using a nonlinear programming solver such as SNOPT and exploiting Haar wavelet operational matrices for the purpose of analyzing the optimal control system [49]. Now our aim is to integrate the Haar wavelets. The generalized fractional integration of Haar vector H = [h 0 , h 1 , h 2 , . . . , h m-1 ] is given as where P α,ρ (t) is a square m-dimensional operational matrix of generalized integrals. In general these generalized fractional integrals can be calculated analytically as This formula holds for i > 1. For i = 1, we obtain The generalized fractional order integration matrix P α,ρ (t) can be obtained by using collocation points in equations (13) and (14). In particular, the Haar wavelet operational matrix for the fixed variables α = 0.75, m = 8, and ρ = 1.5 is Furthermore, in Sect. 4 we will present a new approach to solve certain classes of linear or nonlinear boundary value problems of generalized fractional differential equations numerically, called Green-Haar wavelet method. This technique will be free of operational matrices

Error analysis
An error analysis for a function approximation by Haar wavelets is carried out in [43]. Here we derive an inequality in the context of an upper bound for the Caputo-Katugampola fractional differential operator which shows the convergence of the Haar wavelet technique. The proof of the following theorem is similar to that in [50].
where a, b ∈ R + , and D α,ρ a u m is an approximation of D α,ρ a u, then we have Proof The function D α,ρ a u defined over [a, b] can be approximated as where Consider the first m terms of the sum, denoted by D α,ρ a u m , which approximate D α,ρ a u(t), that is, where m = 2 β+1 , β = 1, 2, 3, . . . , then where I m is the identity matrix of order m. Therefore, From equation (17) we have Recall the mean value theorem of integration: Therefore, Together with the definition of Caputo-Katugampola fractional derivative and By the mean value theorem, ∃ξ ∈ [t 1 , Substituting equation (21) into (20), we get Putting equation (22) into equation (19), we have Let m = 2 β+1 , then we have Hence, one can achieve the error bound for the given partial sum, provided a numerical value of M is given.
To get an estimate for M, we proceed as follows. Since u (n) (t) is bounded and continuous on the interval [a, b], The integral of Haar wavelets of order α is given as Integrating equation (23) yields Similarly, equation (24) is integrated as Therefore, Equation (25) can be written as Writing equation (26) in matrix form, we get By solving the linear system in (27), we can determine the vector C t , and putting this vector into equation (23), u (n) (t) can be calculated for each t ∈ [a, b]. Now assuming x i ∈ [a, b] and u (n) (x i ) are calculated for i = 1, 2, 3, . . . , , where the points are equidistant, an estimate of M may be considered as + max |u (n) (x i )| 1≤i≤ . Clearly, the estimate shall be relatively more precise if increases and is chosen, for example, equal to b.
From equation (28) it is evident that

Numerical examples
In this part we present a few numerical examples which can help us compare the solutions obtained by the numerical methods with exact solutions and solutions by other methods.
Example 3.3 Consider the αth order Cauchy-type generalized fractional differential equation with the initial condition u(1) = 0, where 0 < α ≤ 1 and ρ > 0. It is easy to check that the analytic solution of system (29) is u(t) = (t ρ -1) 2 . To find the approximate solution, we apply Haar wavelets technique to equation (29). Let then computing the α order integral of (30) along with initial condition leads to Putting the values from (30) and (31) into equation (29), we have where F m H m (t) = 2ρ n (3-n) (t ρ-1 ) 2-n -(t ρ -1) 2 -1. After solving (32) for Haar coefficient vector C t m , and using the result in equation (31), we get the required approximate solution. This problem is also solved in [41] by a decomposition formula. The maximum absolute difference of the numerical and exact solutions of equation (29) for distinct values of α, ρ is documented in Table 1. The numerical results are in good agreement with the exact solutions.
Example 3.5 Consider a fractional differential equation with variable coefficient defined as with the initial condition u(1) = 1, where 0 < α ≤ 1. For a(t) = (1 + t) and g(t) = 2 α (2-α) (t ρ -1) 1-α + (1 + t)t ρ , it can be verified that the analytic solution for equation (37) is u(t) = t ρ . We apply the Haar wavelets technique with the aid of Haar matrices while seeking an approximate solution. Let After applying integral I α,ρ 1 + on both sides of equation (38), we obtain Using equations (38) and (39) in equation (37), we obtain where f (t) = g(t) -(1 + t) is approximated as f (t) = F t m H m (t) and a(t)P α,ρ m×m H m (t) = P α,ρ m×m H m (t). For ρ = 2 and several fixed values of α and m, Table 3 contains the maximum absolute error obtained from the exact and approximate results achieved through the Haar wavelets along with the method discussed in [40]. The tabulated results show that the presented method is nearly as accurate as the discretization method but with comparatively fewer nodes.
For a = b = c = 1 and f (t) = ρ α (α + 1) + ρ 3 2 (α+1) (α-3 2 +1) (t ρ ) α-3 2 + (t ρ ) α , the exact solution is u(t) = (t ρ ) α , and for ρ = 1 and n = 2, equation (41) becomes Bagley-Torvik equation considered in [51]. To find an approximate solution, we use Haar wavelets technique as follows. Letting and performing integration I α,ρ 0 on both sides, as well as using initial conditions, we have and In the same way, the input function g(t) can be approximated by Haar functions as Putting equations (42), (43), (44), and (45) into equation (41), we have By solving (46), we can get the Haar coefficients C t m . Then using equation (43), we can obtain the required output u(t). The absolute error is shown in Table 4 for ρ = 1.5, m = 64, and distinct values of α. Example 3.7 Consider a generalized fractional differential equation of inhomogeneous type with boundary conditions: with u(0) = u 0 , u(1) = u 1 , and 1 < α ≤ 2. For a(t) = 1, f (t) = ρ α t ρ + (t ρ ) α+1 (α+2) , u 0 = 0, and u 1 = 1 (α+2) , the analytic solution of the differential equation is u(t) = (t ρ ) α+1 (α+2) . To find a numerical solution, the integral form of equation (47) is given by where Integrating on both sides of equation (49), we have Using (49) and (50) in (48) yields where the approximations g(t) = G t m H m (t) and t ρ I α,ρ m×m H m (t) are used for convenience. To obtain the value of C t m , we have to solve the algebraic linear system in equation (51) and, putting value of C t m into equation (49), we have an approximate solution. The maximum absolute error of the exact and numerical solutions for ρ = 1.5 and different values of α and m is given in Table 6; see Sect. 4.

Nonlinear problems
A nonlinear differential equation can be transformed to a sequence of linear differential equations. One of the possible way to achieve this goal is the application of quasilinearization method. The quasilinearization technique was presented by Kalabas and Bellman [52] as a generalization of a specific method (Newton-Raphson) [53] which assists in solving nonlinear functional equations. Further, Haar wavelets with the quasilinearization technique have been applied for the numerical solution of the individual or system of nonlinear fractional differential equations [44]. Here we apply the quasilinearization technique to solve generalized nonlinear fractional differential equations.

Green-Haar method
The fractional Green's function was defined by K.S. Miller and B. Ross [54] who applied it to fractional differential equations consisting of derivatives of order kα only, where k ∈ Z. In this section we present a numerical method which is based on standard Haar wavelets and the Green function. This method applies to boundary value problems of a certain type. An interesting feature of the method is that in does not require fractional operational matrices or specific metrics reserved for solving boundary value problems. The study undertaken shows that the method is computationally efficient against the standard Haar wavelet technique discussed in the previous section. Thus the efficiency of the method is found to be considerably higher than that of some relevant numerical methods. Interestingly, the accuracy is not compromised, but rather enhanced.
At this stage, we shall consider following class of fractional boundary value problems: where Proof Suppose u(t) is the solution of (60), then the integral form of equation (60) is given by Using the boundary condition, we obtain or where g(t) = u 0 + t ρ (u 1u 0 ) and namely where Conversely, assume that u(t) satisfies the Fredholm integral equation. Taking generalized Caputo-type fractional derivative on both sides of equation (64), we end up with equation (60).
The graph for the function in equation (66) at the values of α = 2, ρ = 1, and m = 64 is shown in Fig. 5.

Linear case
In this part we present numerical solution of a class of generalized linear fractional differential equations with boundary conditions. Example 4.2 We consider the fractional boundary value problem in equation (47) with Dirichlet boundary conditions u(0) = u 0 , u(1) = u 1 .
The integral representation for this boundary value problem is given by or where and (α+2) , u 0 = 0, and u 1 = 1 (α+2) , the analytic solution of the boundary value problem is u(t) = (t ρ ) α+1 (α+2) . Let Using (69) and (70) in (67), we have where , and using the orthonormality of the sequence {h i (t)} on [0, 1), we obtain where I m×m represents an identity matrix of dimension m × m. We can solve the algebraic equation (71) for the Green-Haar coefficient vector C t m and, by equation (69), we get the required numerical solution. For ρ = 1.5 and different values of α and m, the tabular data in Table 6 presents the maximum absolute error given by Haar wavelet method in Example 3.7 and Green-Haar wavelet method, respectively. Green-Haar wavelet technique provides significantly more accurate numerical results in comparison with Haar wavelet technique. Moreover, it is also computationally less intensive and takes less time compared to the Haar wavelet method.  (47) with Dirichlet boundary conditions u(0) = u 0 , u(1) = u 1 , but this time with different data. Particularly, we take g(t) = ρ α ( (α + 1) -(α + 2)t ρ ) + (t ρ ) α -(t ρ ) α+1 . The exact solution of fractional differential equation is u(t) = (t ρ ) α -(t ρ ) α+1 . The corresponding integral representation for the boundary value problem is given by where and Using (76) and (77) in (74), and using the orthonormality of the sequence {h i (t)} on [0, 1) in equation (72), we obtain where f (t) is approximated as f (t) = F t m H m (t). The vector C t m can be obtained by solving the algebraic linear system in equation (78), which leads to the numerical solution when inserted into (76). The numerical solution is in good agreement with the exact solution as shown in Fig. 6. The absolute error for several fixed values of m is given in Table 7. Furthermore, the numerical solutions for ρ = 1.45, m = 64, and different values of α are given in Fig. 7. Figure 8 shows the numerical solutions for m = 64, α = 1.55, and several values of ρ.
where f (t) = F t m H m (t), G t m H m (t) = 1 0 G 2 (t, τ ) dτ , which are calculated by using trapezoidal and Simpson's rules together.
If we choose a 3 (t) = -1 and f (t) = (2α+1) (α+1) (t ρ ) α -(t ρ ) 6α , then the exact solution of equation (79) is u(t) = (t ρ ) 2α . The maximum absolute error of the numerical and exact solutions is given in Table 8 for various fixed values of α, m, and ρ = 1. Also exact and numerical solutions for equation (79) are shown in Fig. 9. Graphical results show that the numerical and exact solutions match with each other.

Error analysis
In this section we establish a bound of the absolute error for the Green-Harr wavelets. (0, 1), that is, there exits M > 0 such that |u (t)| ≤ M for all t ∈ (0, 1), and also assume that u k (t) is an approximation of u(t), Then we have u(t)u k (t) E ≤ ρ 1-α 3 (α + 1)

Theorem 4.5 Suppose that function u (t) is continuous and bounded on
• Green-Haar method has been proposed for numerical solutions of linear and nonlinear fractional boundary value problems. • A comparison has been conducted for the proposed method with conventional Haar wavelet technique. We conclude that Green-Haar method is relatively more efficient and accurate. • The convergence and stability of Green-Haar method have also been discussed.