Chebyshev Wavelet Finite Difference Method: A New Approach for Solving Initial and Boundary Value Problems of Fractional Order

A new method based on a hybrid of Chebyshev wavelets and finite difference methods is introduced for solving linear and nonlinear fractional differential equations. The useful properties of the Chebyshev wavelets and finite difference method are utilized to reduce the computation of the problem to a set of linear or nonlinear algebraic equations. This method can be considered as a nonuniform finite difference method. Some examples are given to verify and illustrate the efficiency and simplicity of the proposed method.

In recent years, wavelets have received considerable attention by researchers in different fields of science and engineering. One advantage of wavelet analysis is the ability to perform local analysis [32]. Wavelet analysis is able to reveal signal aspects that other analysis methods miss, such as trends, breakdown points, and discontinuities. In comparison with other orthogonal functions, multiresolution analysis aspect of wavelets permits the accurate representation of a variety of functions and operators. In other words, we can change and simultaneously to get more accurate solution. Another benefit of wavelet method for solving equations is that after discreting the coefficients matrix of algebraic equations is sparse. So the use of wavelet methods for solving equations is computationally efficient. In addition, the solution is convergent. The operational matrix of fractional order integration for the Chebyshev wavelet, Legendre wavelet, and Haar wavelet has been introduced in [33][34][35] to solve the differential equations of fractional order. A CAS wavelet operational matrix of fractional order integration has been developed by Saeedi et al., to solve fractional nonlinear integrodifferential equations [36,37].
The paper is organized as follows. Section 2 included some necessary definitions and mathematical preliminaries of fractional calculus, Chebyshev polynomials, and Chebyshev wavelets. In Section 3, we introduce the Chebyshev finite difference method. In Section 4, Chebyshev wavelet finite difference method (CWFD) is presented. Section 5 included convergence analysis of the proposed method. In Section 6, the proposed approach is used to approximate the fractional differential equations. As a result the fractional differential equation is converted to a system of algebraic equations which is simply solved. Some illustrative examples of different types are given to demonstrate the efficiency and accuracy of the method in Section 7. In Section 8, concluding remarks are given.

Preliminaries and Notations
In this section, we present some notations, definitions, and preliminary facts that will be used further in this work.

Fractional Integral and Derivative.
There are several definitions of fractional integral and derivatives [2,38,39], including Riemann-Liouville, Caputo, Weyl, Hadamard, Marchaud, Riesz, Grunwald-Letnikov, and Erdelyi-Kober. The most commonly used definition is of Riemann-Liouville and Caputo. One of the drawbacks of Riemann-Liouville method is that it cannot incorporate the nonzero initial condition at lower limit. The Caputo fractional derivative allows the utilization of initial and boundary conditions involving integer order derivatives, which have clear physical interpretations. Therefore, in this study we will use the Caputo fractional derivative proposed by Podlubny [14].

Definition 1.
A real function ( ), > 0, is said to be in the space , ∈ R if there exists a real number > such that ( ) = 1 ( ), where 1 ( ) ∈ [0, ∞), and it is said to be in the space as if and only if ( ) ∈ , ∈ N.

Chebyshev Polynomials.
Chebyshev polynomials of the first kind of degree can be defined as follows: which are orthogonal with respect to the weight function where is the Kronecker delta function and Chebyshev polynomials also satisfy the following recursive formula: The set of Chebyshev polynomials is a complete orthogonal set in the Hilbert space £ 2 [−1, 1]. A function ∈ £ 2 [−1, 1] can be written in terms of Chebyshev polynomials as

Wavelets and Chebyshev Wavelets.
Wavelets have been very successfully used in many scientific and engineering fields. They constitute a family of functions constructed from dilation and transformation of a single function called the mother wavelet ( ); we have the following family of continuous wavelets as [42,43]: If we set = − 0 , = 0 − 0 ; 0 > 1, 0 > 0, we will get the following family of discrete wavelet which forms a wavelet basis for 2 (R): Abstract and Applied Analysis 3 In particular, when 0 = 2 and 0 = 1, then , ( ) forms an orthonormal basis.

Chebyshev Finite Difference Method
Clenshaw and Curtis [45] introduced the following approximation of the function ( ) denoted by ( ) ( ) as where the summation symbol with double primes denotes a sum with both the first and last terms halved. Moreover, are the extrema of the th-order Chebyshev polynomial ( ) and are defined as These well-known Chebyshev-Gauss-Lobatto interpolated points, are the zeros of (1 − 2 )( ( )/ ). Using (4) we have, so can be rewritten as The first two derivatives of the function ( ) at the points , = 0, 1, . . . , are given by [46] as where with 0 = = 1/2, = 1 for = 1, 2, . . . , − 1. As can be seen from (24), the first two derivatives of the function ( ) at any point of the Chebyshev-Gauss-Lobatto points are expanded as a linear combination of the values of the function as these points.

Convergence Analysis
Lemma 5. If the Chebyshev wavelet expansion of a continuous function ( ) converges uniformly, then the Chebyshev wavelet expansion converges to the function ( ) [47].
, with bounded second derivative, say | ( )| ≤ , can be expanded as an infinite sum of Chebyshev wavelets, and the series converges uniformly to ( ); that is, Proof. We havê Abstract and Applied Analysis 5 if > 1, by substituting 2 − 2 + 1 = cos , it yieldŝ Using integration by part, we get the first part is zero, therefore, Using integration by part again, it yieldŝ where Thus, we get However Since ≤ 2 −1 , we obtain Now, if = 1, by using (35), we havê It is mentioned in [42] that { 0 } ∞ =1 form an orthogonal system constructed by Haar scaling function with respect to the weight function ( ), so Therefore, with the aid of Lemma 5, the series Theorem 7. Suppose ( ) ∈ 2 [0, 1) with bounded second derivative, say | ( )| ≤ ; then its Chebyshev wavelet finite difference expansion converges uniformly to ( ); that is, where the summation symbol with prime denotes a sum with the first term halved.
Proof. From Theorem 6, we have wherê, This series converges to ( ) uniformly. We first show that converges tô. We havê 6 Abstract and Applied Analysis by substituting 2 − 2 + 1 = cos , it yieldŝ Using the trapezoidal rule for integration with equidistant subintervals giveŝ From (26) and (28), for 1 < < , we havê According to approximation error for the trapezoidal rule, we havê− In view of (35) and triangle inequality, we get where Because | ( )| ≤ 2 /2 , it is understandable that Therefore, in view of Lemma 5, series (43) is uniformly convergent to ( ).
Moreover, we should impose continuity condition on the approximate solution and its first derivatives at the interface between subintervals which results in (2 −1 − 1)( + 1) equations.
We will totally have a system of 2 −1 ( + 1) algebraic equations, which can be solved for the ( ). Consequently, we obtain the solution ( ) to the given (60) using (27) and (28).

Illustrative Examples
In this section, we consider some numerical examples for the fractional equation to demonstrate the validity of the proposed method (CWFD) in solving fractional differential equation. These examples are considered because closed form solutions are available for them, and they have been solved using other numerical methods. This allows one to compare the results obtained using this method with the analytical solution and the solutions obtained using other methods.
such that the second initial condition is for > 1 only. The exact solution of (70) and (71) is given as [26] ( ) = 8 − 3 4+ /2 + 9 4 . (72) It should be noted that for < 1, the slope of the solution at = 0 goes to infinity. Therefore, one can expect a large numerical error near = 0.
We applied the method introduced in Section 6 and solved this problem for different values of and , . Figure 1 shows the analytical and numerical results for = 9, 6, 3 and = 1 and = 0.75, 1.5. It can be seen that increasing the values of and results in more accurate solution. In Figure 2 numerical results for = 12, = 4, and = 0.5, 0.75, 0.95 and = 1 (exact solution) and also = 1.5, 1.75, 1.95, and = 2 (exact solution) are plotted. We see that as approaches 1 and 2, the solution of the fractional differential equation approaches to that of the integer-order differential equation. In Table 1, we compare the results obtained by the present method using = 12, = 6 with those in [50]. As can be seen, our results are much more better than those obtained by Saadatmandi and Dehghan [50]. Furthermore, as it is expected, the absolute error is reduced as approaches an integer value.
Example 2. As the second example, we consider the following initial value problem in the case of the inhomogeneous Bagley-Torvik equation [51]: where ( ) = 1 + subject to the following initial states The exact solution of this problem is ( ) = 1 + . We solve this fractional initial value problem by applying the method described in Section 6 with = = 3. The absolute error is shown in Figure 3, which shows that the numerical solution is in perfect agreement with the exact solution. We set Digits = 50. However, we get the absolute error less or equal to 10 −18 when we set Digits = 20.
Abstract and Applied Analysis 9   For = = = = 1, 1 = 0.333, 2 = 1.234, in Table 2, we compare the absolute errors in the solution with those obtained using Chebyshev wavelet method [33]. It can be seen that our results are much more accurate. We also show absolute errors for different values of and in Table 3. It is observed that the absolute errors in solution are reduced by increasing the values of and .
where the exact solution is ( ) = 2 . It is worth mentioning that the method presented above only can be employed for solving this problem for ∈ [0, 1]. That is because the first kind Chebyshev wavelet is defined on interval [0, 1]. However, the variable in this example is defined on interval [0, 50], and we should replace Chebyshev   wavelet bases ( ) with ( /50) in the discrete procedure. We applied the method described in Section 6 and got almost exact solution for ∈ [0, 50], although other authors had solved the problem for ∈ [0, 1] and ∈ [0, 5]. It can be seen from Figure 4 that the approximate solution and exact solution are closely overlapped for any ∈ [0, 50].
Example 5. Consider the following boundary value problem for nonlinear fractional order differential equation: where 1 < ≤ 2, , ∈ R, ∈ N, and ( ) is a given problem. For = 1.5, = 2, = 0, = 1, and ( ) = ((15√ )/8Γ((7/2)− )) (5/2)− + −2 5 , it can be easily verified that the exact solution is ( ) = 5/2 . We use the introduced scheme in Section 6 to solve this example with = 20, = 5. Absolute error in solution for different values of is presented in Table 4 which confirm the accuracy and efficiency of the proposed method. As it is expected, the absolute error is reduced as approaches an integer value. Furthermore, we plot exact and numerical solutions for = 2 in Figure 5.
where 1 < ≤ 2, 0 < ≤ 1 and ( ) = −1 + 1. The exact solution of the problem is not known generally. It can be easily verified that for = 2, = 1, the exact solution is ( ) = (1 − −1 ). The problem (80) is solved numerically for integer order case in [55,56] using Haar wavelet method and combined homotopy perturbation method, respectively. We solve the problem using = 12, = 3. In Table 5, the absolute errors are presented which confirm that the proposed method is more accurate. The numerical results for = 1 and different values of plotted in Figure 6 show that as tends to 2, the solution of fractional differential equation approaches to that of the integer-order differential equation.  Example 7. Consider the following linear fractional differential equation [26,49,50,57]: such that The condition (0) = 0 is only for 1 < ≤ 2. The exact solution of (74) and (75) is given by where ( ) = ∑ ∞ =0 /Γ( + 1) is the Mittag-Leffler function of order . It can be easily verified that for = 1 and = 2, the exact solutions are ( ) = exp(− ) and ( ) = cos( ), respectively.
We solved the problem using the proposed method for different values of .
The absolute errors for = 0.2, 0.4, 0.6, 0.8 (with = 6, = 6) and = 1.2, 1.4, 1.6, 1.8 (with = 6, = 8) are shown in Table 6. It can be seen that our results are more accurate than the ones reported in [50]. The numerical solutions for = 0.5, 0.75, 0.95, and 1 (Figure 7(a)) and = 1.5, 1.75, 1.95, and 2 (Figure 7(b)) are plotted in Figure 7. It should be noted that as approaches 1 and 2, the numerical solutions converge to the analytical solution ( ) = exp(− ) and ( ) = cos( ), respectively.  The exact solution of this problem is ( ) = √ ( + 1). This problem was solved in [58] by operational matrix of fractional derivatives using B-spline functions. We compare maximum absolute errors obtained by the introduced method with the ones reported in [58] in Table 7. It should be noted that the algebraic system of equations obtained by the method [58] is of order 2 + 1, while the resulting one by the current method is of order 2 −1 ( +1). It can be seen from Table 7 that   our results are more accurate while we need to solve a system of algebraic equations of lower order.

Conclusion
An efficient and accurate method based on hybrid of Chebyshev wavelets and finite difference methods was introduced. The fractional derivative is described in the Caputo sense.   The useful properties of Chebyshev wavelets and finite difference method make it a computationally efficient method for solving the problems. The main advantage of the present method is the ability to represent smooth and especially piecewise smooth functions properly. Several examples are given to demonstrate the powerfulness of the proposed method. We note that the accuracy can be enhanced either by increasing the number of subintervals or by increasing the number of collocation points in subintervals properly. The validity and accuracy of the method were investigated for large intervals as well.