Quasilinearized Semi-Orthogonal B-Spline Wavelet Method for Solving Multi-Term Non-Linear Fractional Order Equations

: In the present article, we implement a new numerical scheme, the quasilinearized semi-orthogonal B-spline wavelet method, combining the semi-orthogonal B-spline wavelet collocation method with the quasilinearization method, for a class of multi-term non-linear fractional order equations that contain both the Riemann–Liouville fractional integral operator and the Caputo fractional differential operator. The quasilinearization method is utilized to convert the multi-term non-linear fractional order equation into a multi-term linear fractional order equation which, subsequently, is solved by means of semi-orthogonal B-spline wavelets. Herein, we investigate the operational matrix and the convergence of the proposed scheme. Several numerical results are delivered to conﬁrm the accuracy and efﬁciency of our scheme.


Introduction
Fractional calculus, a generalization of conventional integer calculus, has been found to be more appropriate for describing important phenomena in the fields of dynamics [1], physics [2,3], medicine [4], chemistry [5], and other scientific areas [6]. There have been plenty of works, especially those considering real-world physical systems, where the Caputo-type fractional derivative and the Riemann-Liouville fractional integral have been used to describe complex dynamical systems. As most fractional order equations cannot be resolved analytically, numerical methods have been developed to give their numerical solutions [7][8][9][10]. In recent years, multi-term fractional order equations have received increasing attention due to their wider applicability. Many papers have been devoted to numerical methods for approximating multi-term fractional differential equations, such as the Bernstein polynomials method [11] for deriving the operation matrix of the Caputo derivative based on the collocation method, and the B-spline method [12] for constructing an operation matrix based on linear B-spline functions. In contrast, only a few papers have focused on the development of numerical techniques for multi-term fractional order equations with fractional integral and derivative. In [13], Kojabad and Rezapour theoretically discussed the existence of solutions of multi-term fractional order equations, then numerically verified that the differences between the results of the Legendre method and the Chebyshev method are negligible; however, they did not provide a method with greater accuracy and speed for multi-term fractional order equations. Zheng et al. [14] studied linear multi-term fractional initial problems using the discontinuous Galerkin finite element method; however, the study of numerical methods for non-linear multi-term fractional order equations, which can describe material transport behaviors in complex physical systems, is scarce. This motivated us to propose an efficient and accurate numerical method.
In this paper, we consider the following Problem I about multi-term fractional order equations: with p + 1 boundary conditions in [0, b], g k 0, u(0), · · · , u (p) (0) = 0, k = 1, · · · , j − 1 (2) and g k b, u(b) · · · , u (p) (b) = 0, k = j, · · · , p + 1, where 0 < α 1 < α 2 < · · · < α q , β ≥ 0, p < α q ≤ p + 1, and p, q ∈ N. The functions f and g k (k = 1, · · · , p + 1) are non-linear functions of u(x) and its derivatives u k (x)(k = 1, · · · , p). The operator D α i is the Caputo-type fractional derivative operator of order α i and I β is the Riemann-Liouville fractional integral operator of order β. The definitions of the Caputo fractional derivative and the Riemann-Liouville fractional integral are introduced in the next section. Wavelet methods, which are relatively novel approaches, have been applied to solving problems in fractional calculus. Due to the structural characteristics that the basis of wavelets can achieve through the dilation and translation of mother wavelet functions, wavelet methods can enhance numerical efficiency. Some orthogonal wavelet methods have been developed to estimate approximate solutions for fractional order equations [15,16]. Compared with orthogonal wavelet methods, the semi-orthogonal B-spline wavelet method (SOBWM) has the following advantages: compact support, explicit analytical form, and finite basis functions in any wavelet subspaces [17]. Due to its accuracy and efficiency, SOBWM has been applied in several kinds of differential and integral equations. Maleknejad et al. [18] adopted this kind of method to solve non-linear Fredholm-Hammerstein integral equations of the second kind, Aram et al. [19] used the method to deal with integro-difference equations, while Liu et al. [20] approximated multi-term linear fractional differential equations by this method.
The quasilinearization method [21,22] and the homotopy method [23,24] are common approaches for linearizing non-linear functions. For multi-term fractional order equations, however, our simulation results show that the solutions obtained by the homotopy method easily diverge. By contrast, the quasilinearization method is more suitable for non-linear multi-term fractional order equations. In this paper, we utilize the quasilinearization method to linearize the non-linear fractional order equations in Problem I. Then, we solve the linearized fractional order equations using the semi-orthogonal B-spline wavelet collocation method (SOBWCM).
The remainder of this article is structured as follows: In Section 2, we introduce some basic notation, definitions, and lemmas in fractional calculus. The definition of the semi-orthogonal B-spline wavelets (SOBW) and related theorems and properties are given in Section 3. Section 4 presents the implementation process of the quasilinearized semi-orthogonal B-spline wavelet method (QSOBWM). The convergence of QSOBWM is analyzed in Section 5. In Section 6, the validity of the presented scheme is examined through illustrative examples. In Section 7, we draw a concise conclusion.

Definitions and Properties of Fractional Calculus
In this section, the definitions of the Riemann-Liouville fractional integral and the Caputo fractional derivative are presented. Then, the definition and some properties of the Laplace transform are provided, which is applied to the proof in the next section, due to its ability to simplify the solving of differential equations by reducing them to algebraic equations.

Definition 3 ([25]
). The Caputo fractional derivative of order α > 0 for a function u(x) ∈ C n −1 is defined as where n ∈ N + , u (n) (x) is the conventional integer derivative of order n for u(t).

Definition 4 ([26]
). The functionû(s) of the variable s is defined bŷ which is called the Laplace transform of the function u(x). where c 0 lies in the right half-plane of the absolute convergence of the Laplace integral Equation (6).

Semi-Orthogonal B-Spline Wavelets in [0, b]
In this section, we describe the construction of SOBW in [0, b]. The integral and derivative for SOBW in [0, b] are also presented.

Construction of Semi-Orthogonal B-Spline Wavelets
The m-th order B-spline function is defined by Here, [·, · · · , ·] ξ is the m-th divided difference of (x − ξ) m−1 + with respect to the variable ξ and Let j 0 be the smallest integer satisfying which is the minimum to contain one complete wavelet function in [0, b].

Fractional Integral and Derivative of Semi-Orthogonal B-Spline Wavelets
From Definitions 6-8, the fractional integral and derivative of SOBW can be reduced to solving the fractional integral and derivative of function (ax − b) k + , where k > 0, a > 0, and b > 0.

Function Approximation
The quasilinearization approach [29] is applied to approximate u(x) in Problem I. Concisely, we note with linearized non-linear boundary conditions and where u (0) The initial value, u 0 (x), can be selected from mathematical or physical conditions, and u r+1 (x) is further obtained by iteration.
Equations (24)- (26) are linear equations in u r+1 (x) and, so, Problem II consists of multi-term linear fractional order equations that can be addressed efficiently by SOBWCM. For simplicity, Problem II is reformulated into the equivalent Problem III as with the boundary conditions , and b(x) are related to the functions and values of r-th iteration by the quasilinearization approach. A can be expanded by semi-orthogonal B-spline scaling functions and wavelets [27] as To meet the needs of practical applications, the higher frequency components are truncated at M, such that the infinite series in Equation (30) is approximated as where C and Ψ are the (2 M+1 + m − 1) × 1 vectors Substituting u r+1 (x) from Equation (31) into Equations (25)- (29) of Problem III, we obtain: with the boundary conditions In order to increase the computational efficiency, Equation (32) can be rewritten, in matrix form, as: where and The terms D α i Ψ and I β Ψ can be implemented on a computer, based on Lemmas 3 and 4. Therefore, Problem III is converted into the solution of a system of linear equations, which consists of Equations (33)-(35).

Convergence Analysis
The QSOBWM is a hybrid numerical method that combines the quasilinearization method and SOBWCM. Thus, we analyzed the convergence of the two basic methods, followed by the convergence of the method proposed in this paper, obtained in the process of function approximation.
, l, n = 0, 1, · · · , p, and the difference function δu r+1 (x) = u r+1 (x) − u r (x) is considered in the quasilinearization method for Problem II. Then, there exists a positive constant C M such that where δu r is the maximum value of any of δu Proof. By iteration of Equations (24)-(26), we have with the corresponding boundary conditions: According to the mean value theorem in [30], r−1 (x) and u In view of the n-term linear fractional Green's function properties in [26], there exists a Green function G(x, y) such that Therefore, r (y) f u (l) u (n) y,ū r−1 (y), . . . ,ū r−1 (y) , and p such that Hence, The convergence of the quasilinearization method is shown in Theorem 1. Then, to estimate the error of SOBWCM, we use the following theorem from [20]: is approximated by SOBWCM of order m in Problem III. Then, the truncation error for j = M is Theorem 2 implies that |u r+1 (x) −ũ r+1 (x)| → 0 when M → ∞, reflecting the convergence of SOBWCM. As each iteration step (i.e., from u r (x) to u r+1 (x)) and the process of approximating u r+1 (x) are convergent, QSOBWM proposed in this paper is convergent.

Numerical Examples
This section presents some numerical examples to demonstrate the validity of the proposed scheme. Computation was carried out on a personal computer with an Intel(R) Core(TM) i5-7500 CPU @ 3.40 GHz with 8.00 GB RAM and the codes were written in Matlab 2014a. We denote u as the exact solution and u as the numerical solution. To assess the performance of the method, we calculated the absolute error, L 2 error, and L ∞ error.
The absolute error in [0, b] is the L 2 error is defined as and the L ∞ error is defined as where N x is the number of collocation points in [0, b].

Example 1.
Consider the fractional integro-differential equation with weakly singular kernel [8,9]: with the initial condition u(0) = 0, The exact solution of the problem is We adopted the QSOBWM with m = 4 for several truncation values M, and the absolute errors for each case are exhibited in Table 1. The numerical results of the proposed scheme in Table 1 illustrate that the absolute errors decreased when the value of truncation M increased. More intuitively, we describe the absolute errors of the present scheme for various values of M in Figure 1. The results were in accordance with the convergence analysis of the present scheme.
In order to compare with the second-kind Chebyshev polynomials method (SKCPM) in [8] and the fractional order Euler functions method (FEFsM) in [9], we computed the L 2 errors of the present scheme with various values of M and list the results in Table 2, where N denotes the maximal degree of the polynomials in the space spanned by all polynomials for SKCPM and FEFsM, and M indicates the truncation in the present scheme. The degree of all polynomials was no more than three when m = 4 in the present scheme. Table 2 shows that the L 2 errors of the present scheme with m = 4 were smaller. Example 2. Consider the following fractional Langevin equation [14]: with the initial condition u(0) = 0, Γ(2.0+β) and α 1 , α 2 , β ∈ (0, 1). The exact solution is given by u(x) = x 2 − x.   We solved this problem by using the QSOBWM with m = 4 and M = 3 for fixed step size h = 1 20 . Table 3 exhibits the absolute errors, which demonstrates that the absolute errors of various values of α 1 , α 2 , and β were less than 10 −13 , and that the computation took only 4.10 s. Furthermore, the results confirmed the effectiveness of this method. Figures 2 and 3 depict the absolute errors of the present scheme for different α 1 and β, respectively. Figure 2 presents the absolute errors of the present scheme with α 2 = 0.8 and β = 0.3 for various values of α 1 , while Figure 3 displays the absolute errors of the present scheme with α 1 = 0.5 and α 2 = 0.8 for several values of β. These graphs demonstrate that the approximate solutions agreed closely with the analytical result for various values of α 1 and β, with absolute errors less than 10 −13 .    Table 4. Comparison of the errors of QSOBWM and the method in [14] for Example 2. In order to verify the performance of QSOBWM, as shown in Table 4, we calculated the L 2 error and L ∞ error of the present scheme and the method in [14] for α 1 = 0.5, α 2 = 0.8, and β = 0.3 with different grid sizes. The results demonstrated that the approximate results for the present scheme were closer to the analytical solutions than the method in [14].

Errors
Some interesting phenomena are shown in Figures 1-3. There may be two reasons: First, when x = 0, the error is much closer to zero as the value at x = 0 satisfies both the original equation and the initial value; however, the values with x ∈ (0, 1] only satisfy the original equation, which is more likely to have bigger errors. Second, it might be related to the reduced algebraic equations, which needs more exploration. Example 3. Consider the following fractional Duffing-Holmes model for a non-linear oscillator [31]: with the initial value: From [31], u(x) = x µ+2 sin x is the given exact solution.
We approximated the solution using QSOBWM with m = 4 and different M (M = 3, 4, 5, 6). This problem was also solved by the spectral collocation method (SCM) in [31] with µ = 0.5 and N = 5, 6, 7, 8, where N denotes the maximal degree of polynomials in the space formed by the basis of the method. In Table 5, we compare the L 2 errors and L ∞ errors of the two methods mentioned above. The degree N of the polynomials in the present scheme was three. The L 2 and L ∞ errors of the present scheme were much smaller than those of the SCM from Table 5. Moreover, the L 2 errors and L ∞ errors of the present scheme decreased when M increased. Thus, the results were in accordance with the convergence analysis discussed in the previous section. Example 4. Consider the following multi-term fractional non-linear boundary value problem [11,12]: with boundary conditions: The exact solution of this problem is u(x) = x 2 + 1. Table 5. Example 3: L 2 errors and L ∞ errors with µ = 0.5 for spectral collocation method (SCM) and QSOBWM.  The problem was resolved by the B-spline operational matrix method (BSOMM) [12] and the Bernstein operational matrix method (BOMM) [11], respectively. We also solved the example by applying the present scheme with m = 4 and list the values of the L 2 and L ∞ errors of the three methods in Table 6. The L 2 errors of the present scheme achieved 10 −14 within 7.98 s. As shown in Table 6, the present scheme was the most accurate. Table 6. The L 2 errors and L ∞ errors of QSOBWM, B-spline operational matrix method (BSOMM) (in [12]), and Bernstein operational matrix method (BOMM) (in [11]) for Example 4.

Conclusions
In this article, the quasilinearized semi-orthogonal B-spline wavelet method is proposed to approximate multi-term non-linear fractional order equations containing fractional integrals and derivatives. The proposed scheme significantly reduces the computational complexity of solving non-linear fractional order equations by using the semi-orthogonal B-spline wavelet collocation method. The solution procedure for approximating non-linear fractional equations is described. we also discuss the convergence properties of the present scheme. As the initial and boundary conditions are both considered during the process of function approximation, this method can be applied to solving both initial and boundary value problems of fractional order. Illustrative examples and comparison results were given to demonstrate the efficiency and accuracy of the proposed scheme.