Numerical Investigation of Fractional-Order Differential Equations via 
 φ
 -Haar-Wavelet Method

In this paper, we propose a novel and efficient numerical technique for solving linear and nonlinear fractional differential equations (FDEs) with the 
 
 φ
 
 -Caputo fractional derivative. Our approach is based on a new operational matrix of integration, namely, the 
 
 φ
 
 -Haar-wavelet operational matrix of fractional integration. In this paper, we derived an explicit formula for the 
 
 φ
 
 -fractional integral of the Haar-wavelet by utilizing the 
 
 φ
 
 -fractional integral operator. We also extended our method to nonlinear 
 
 φ
 
 -FDEs. The nonlinear problems are first linearized by applying the technique of quasilinearization, and then, the proposed method is applied to get a numerical solution of the linearized problems. The current technique is an effective and simple mathematical tool for solving nonlinear 
 
 φ
 
 -FDEs. In the context of error analysis, an exact upper bound of the error for the suggested technique is given, which shows convergence of the proposed method. Finally, some numerical examples that demonstrate the efficiency of our technique are discussed.


Introduction
Fractional differential equations are used to describe a wide range of phenomena in natural science, and because of its numerous applications in physical, chemical, and biological sciences, fractional calculus has captivated the scientific community. Several researchers have recently focused their attention on the concept of the fractional derivative. The fractional derivative is introduced in fractional calculus through the fractional integral. Riemann, Liouville, Caputo, Hadamard, Grunwald, and Letinkow are the pioneers in this field, having contributed and published extensively on the subject. The nonlinear fractional Schrodinger equations with the Riesz space and the Caputo time-fractional derivatives are studied using the finite difference/spectral-Galerkin method in [1]. For the Higgs boson equation in the de Sitter spacetime, a finite difference/Galerkin spectral scheme was introduced in [2] which retains the discrete energy dissipation property. For the two-dimensional fractional wave equation with the Weyl space-fractional operators, Ref. [3] proposes a highorder compact difference method with fourth-order precision in space and second order in time. Explicit solutions to differential equations of complex fractional orders with respect to functions and continuous variable coefficients are determined in [4]. Different types of fractional derivatives have appeared in the literature that strengthen and generalize the classical fractional operators defined by the aforementioned authors [5,6]. Katugampola recently discovered a new type of fractional integral operator which encompasses the Riemann-Liouville and Hadamard operators in a single form [7,8]. Moreover, several other fractional operators are being introduced to date. Due to a wide range of definitions for fractional-order integrals and derivatives [9][10][11], the idea of a fractional derivative of one function with respect to another function emerged. This class of fractional operators depends on a kernel function and unifies many definitions of fractional operators. Almeida used the idea of fractional derivatives in the Caputo sense and introduced the φ-Caputo fractional derivative of one function with respect to some other function [12]. The proper choice of a trial function helps in the modeling of physical phenomenon and makes the approach more suitable from the application point of view [13,14].
Wavelet analysis is a well-known and widely used mathematical method in engineering and other sciences [15,16]. Wavelets are made up of function expressions that have been extended into a sum of basic functions. A mother wavelet function is translated and compressed to obtain these basic functions. As a result, it inherits properties of locality and smoothness, making it simple to research the properties of integer and locality during the process of expressing functions. Wavelets have sparked a lot of interest in using them to solve classical ordinary and partial differential equations numerically. Researchers have recently succeeded in extending several standard wavelet methods to numerical solutions for fractional differential equations. Numerical integration and numerical solutions of fractional ordinary and fractional partial differential equations are some of the other applications of wavelet methods in applied mathematics. So, for now, wavelets such as the Haar-wavelet, B-spline, Daubechies, and Legendre wavelet are used [17][18][19][20][21]. In Ref. [22], the Genocchi wavelet-like operational matrix was used together with the collocation method to solve nonlinear FDEs. For solving fractional integrodifferential equations, the Jacobi wavelet operational matrix of fractional integration is constructed and utilized in [23]. The Haar-wavelet is a simple form of orthonormal wavelets with compact support and has been used by many researchers. The Haar-wavelet family consisted of rectangular functions. It also includes the lower member of the Daubechies wavelet family, which is suitable for computer implementation. The Haarwavelets are used to transform a fractional differential equation into an algebraic structure of finite variables [24][25][26][27].
For modeling different physical problems, it is difficult to pick the right operator. Therefore, generalized operators of fractional order should be developed for which classical operators are special cases. An effective way to deal with such a variety is to merge these definitions into one by considering fractional derivatives of function f with respect to another function φ. The Riemann-Liouville operators of fractional order are generalized by introducing the fractional-order differentiation and integration of a function by another function [28,29]. In [12,30], Almeida defined the φ-Caputo fractional differential and integral operators and discussed its characteristics. The contribution made by Almeida et al. plays a pivotal role in putting together a wide range of fractional operators. Moreover, recent work on the φ-Caputo derivative indicates that φ-Caputo fractional differential-based mathematical models are more flexible and provide felicitous results in many situations. In order to evaluate the growth of the world population, Almeida [12] implemented the φ-Caputo derivative and illustrated that the appropriate selection of a fractional operator determines the model's precision. Using fixed-point theorems, Almeida et al. in [13] investigated the existence and uniqueness of a solution for nonlinear FDEs involving a φ-Caputo derivative. Almeida et al. in [31] introduced the φ-shifted Legendre polynomials for solving fractional oscillation equations containing the φ-Caputo derivative of fractional order. We therefore see the theory of φ-FDEs as a promising field for further study. In this paper, taking motivation by the work cited above, we developed a new numerical method for solving linear and nonlinear boundary value problems in φ-FDEs.
The rest of the paper is organized as follows: We start Section 2 with an overview of the fractional calculus followed by a discussion of the classical Haar-wavelet and an approximation of the functions by the Haar-wavelet. In Section 3, we developed the φ-Haar operational matrix of fractionalorder integration of the Haar-wavelet and then utilize it for a numerical solution of the φ-FDEs. In Section 3.1, the error estimate of the developed technique is discussed in depth. Section 4 is devoted to some numerical results and figures that show the precision and effectuality of the developed technique. Finally, a conclusion is given in the last section.

Preliminaries
Here, we present some vital definitions of φ-fractional operators and their basic properties which will be used in the subsequent sections of the paper.

Existence and Uniqueness of Solution for Nonlinear φ-
FDEs. In this section, we provide existence and uniqueness theorems for nonlinear φ-FDEs.
Consider the nonlinear φ-FDE: We have the initial conditions, namely, yðαÞ = y α and y (9) if and only if y satisfies the following fractional integral equation: Theorem 7. Let f be a Lipschitz continuous function with respect to the second variable, that is, ∃ is a positive constant L such that Then, there is a constant h ∈ ℝ + , such that there exists a unique solution to problem (9) on the interval ½α, α + h ⊆ ½α , β.
Proof of Theorems 6 and 7 can be seen in [13].

Approximation of Function by the Haar-Wavelet. The ith
Haar-wavelet defined on the interval ½α, β is given by where Any function yðζÞ defined and square integrable over the interval ½0, 1Þ can be expressed in terms of the Haar-wavelet as follows: where the coefficients c i of the Haar-wavelet are defined by In practice, only the first m terms of the series in equation (14) are considered, where m is a power of 2, that is, with vector form as where

The φ-Haar-Wavelet Operational Matrix
In this section, our endeavor is to construct the φ-Haarwavelet operational matrix P ρ,φ of fractional-order ρ and use it to solve φ-FDEs numerically. The φ-fractional integration of the Haar-wavelet is performed using equation (2). Mathematically, the generalized fractional-order integration of the Haar-wavelet, Analytically, these generalized φ-fractional integrals can be approximated as follows: where Equation (19) holds for i > 1; for i = 1, we have The fractional-order φ-Haar-wavelet operational matrix P ρ,φ for the function φðζÞ = ζ 2 and ρ = 0:75 is given by Also, the approximate and exact φ-RL fractional integration of φðζÞ = sin ð5ζÞ for J = 6 and various choices of ρ is plotted in Figure 1.

Convergence Analysis of the φ-Haar-Wavelet Method.
In [33], the Caputo-type FDEs were recently analyzed for error. Furthermore, utilizing the Haar wavelet, [34] proves convergence for the solution of the nonlinear Fredholm integral equations. In the present work, the upper limit for the error estimate is calculated using the φ-Caputo fractional differential operator. The φ-Haar-wavelet method for FDEs is shown to be convergent.
Proof. D ρ,φ α y can be approximated as follows: where Suppose that D ρ,φ α y m is the following approximation of D ρ,φ α y: where m = 2 κ+1 , κ = 1, 2, 3, ⋯. Then, which implies that By orthogonality of the sequence fh m ðζÞg, we have Ð β α h m ðζÞh m ðζÞdζ = I m , where I m is the identity matrix of order m. Therefore, From equation (25), we have By the mean value theorem for integration, we have ∃ζ 1 , ζ 2 ∈ ðα, βÞ, such that Hence,

Journal of Function Spaces
Employing the definition of the φ-Caputo fractional derivative, the fact that φ is increasing and the condition j y ½n φ ðζÞj ≤ K, we arrive at where Since ζ 1 > α, ζ 2 > α, and ζ 2 > ζ 1 and φðζÞ is an increasing function, so Therefore, By mean value theorem, ∃ϰ ∈ ½ζ, ζ 2 ⊆ ½α, β, such that φ Putting equation (38) into equation (32), we get Equations (29) and (39) give which implies that Let k = 2 κ+1 ; (41) can also be written as From the value of K, we can get an upper bound for the error.
We first estimate the value of K. Since y n ðζÞ is continuous and bounded on ½α, β, so y ½n φ ðζÞ is also continuous and bounded on ½α, β and is given by Similarly,
Theorem 9 can be proven simply via Theorem 8. From equation (49), we can understand that kyðζÞ − y m ðζÞk E tends to 0 as m tends to ∞, which shows that the φ-Haar-wavelet technique converges.
Example 10. To demonstrate optimality of the upper bound in equation (49), we consider the following φ-FDE: with initial condition yð0Þ = 0, having the exact solution yðζ Þ = ðφðζÞÞ 2 ρ. Table 1 shows the optimal values of the upper bound of error obtained for various options J and ρ = 0:25.

Numerical Solutions of φ-FDEs
In this section, we provide the solution to some problems in linear and nonlinear φ-FDEs by employing the φ-Haarwavelet operational matrix technique. Example 11. Consider the composite oscillation equation of a fractional order with the φ-Caputo fractional derivative: with the initial condition yð0Þ = 0. The exact solution of equation (51) is given by yðζÞ = ðφðζÞÞ 2 . For numerical solutions, we approximate D Integrating equation (52) with the φ-Caputo integral operator, we have Substituting the initial conditions in equation (51), we get Substituting equations (52) and (54) for equation (51), we have where f ðζÞ = ðφðζÞÞ 2 + ð2/Γð3 − ρÞÞðφðζÞÞ 2−ρ − ðφð0ÞÞ 2 . Equation (55) can be expressed in the matrix form as follows: where F = f ðζÞ: The value of C T m can be obtained from equation (56). By using C T m in equation (54), we can obtain the numerical solutions. Table 2 where 0 < ρ ≤ 1, ζ ∈ ½0, 1, and the initial condition The exact solution of the problem (57) is given as follows: yðζÞ = ðφðζÞÞ 4 − ð1/2ÞðφðζÞÞ 3 . To get numerical solutions, the φ-Haar-wavelet method is employed as follows. Let Integrating equation (59) with the φ-Caputo integral operator, we have Substituting initial conditions in equation (60), we get c 1 = y 0 . Equation (60) becomes     (57), we get where FðζÞ = ðφðζÞÞ 4 − ð1/2ÞðφðζÞÞ 3 − ð3/Γð4 − ρÞÞ ðφðζÞÞ 3−ρ + ð24/Γð5 − ρÞÞðφðζÞÞ 4−ρ . Approximate solutions are obtained by solving equations (61) and (62). The exact solution, approximate solutions, and the maximum absolute error are plotted in Figure 3 for J = 6 and ρ = 0:6. Also, the maximum absolute errors obtained for various choices of ρ and J are given in Table 3. We noticed that the maximum absolute error decreases with an increase in J.

Nonlinear Case
Example 13. Consider the fractional-order Riccati differential equation with the φ-Caputo fractional derivative: subject to the initial condition yð0Þ = 0. For ρ = 1, the exact solution of equation (63) is given by yðζÞ = ðe φð2ζÞ − 1/e φð2ζÞ + 1Þ. For numerical solutions, we first utilize the quasilinearization techniques to make the nonlinear terms of equation (63) linear and then solve the linearized problem with the φ -Haar-wavelet method. The linearized form of (63) is with the initial condition y r+1 ð0Þ = 0: Now, we apply the φ-Haar-wavelet method to equation (64). Let Operating the φ-Caputo integral on equation (65), we get Putting the initial conditions in equation (66) gives   The matrix form of equation (68) is given by where F = 1 + y 2 r : By solving the algebraic system given by equation (69) for C T m and substituting this value into equation (67), we will have the required numerical solution. In Table 4, the maximum absolute error is given for φðζÞ = ζ 3 . Also, approximate solutions for different choices of the function φ are plotted in Figure 4.
Then, we subject this to the initial condition: Taking the φ-Caputo integral of (73), Substituting the initial condition in equation (74), we have c 1 = 0: Using c 1 = 0 in equation (74), we get Substituting equations (73) (75). Table 5 shows that the maximum absolute error decreases by increasing the values of J. Also, the approximate solutions are displayed in Figure 5 for various values of ρ.
Example 15. For comparison with another method, we consider the following problem: with the initial condition yð0Þ = 0. The exact solution of equation (77) is given by yðζÞ = ðφðζÞÞ 2 . This problem is studied in [31] by using the operational matrix of the φ -shifted Legendre polynomials. For φðζÞ = ðζ 2 /2Þ + ðζ/2Þ, a comparison of the results obtained in [31] and by the proposed method is given in Table 6.

Conclusion
In this article, the φ-FDEs are solved numerically by introducing the φ-Haar-wavelet operational matrix of integration of fractional order. This operational matrix has been used to solve both linear and nonlinear problems with success. In comparison to the other methods, this approach is simple and more convergent. The developed method is used to solve a number of linear and nonlinear problems, demonstrating its efficiency and accuracy. Furthermore, the method's error analysis is thoroughly examined. As a future work, the proposed method may be applied to different wavelets as well as other operators.

Data Availability
The numerical data used to support the findings of this study are included within the article.