Abstract

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.

1. Introduction

The study of fractional calculus dates back to 17th century, starting by G. W. Leibnitz (1695, 1697) and L. Euler (1730) [1, 2] and then has been developed by many researchers in different disciplines. In the year 1823, Liouville and Abel introduced the theory of fractional derivatives and integrals; for more details, please refer to [3, 4]. Fractional calculus has received much more attention from scientists and engineers in recent years. Many researchers in various fields found that derivatives of noninteger order are useful for the description of some natural physics phenomena and dynamic system processes such as damping laws and diffusion process [5, 6]. In general, it is difficult to solve fractional differential equations analytically. Therefore, it is necessary to introduce some reliable and efficient numerical algorithms to solve them. During the past decades, an increasing number of numerical methods were being developed. These methods include homotopy analysis method [7], homotopy perturbation method [810], variational iteration method [913], finite difference method [5, 1418], Adomian decomposition method [1923], fractional differential transform method [24, 25], predictor-corrector method [26], fractional linear multistep method [27], extrapolation method [28], integral transform [29], and generalized block pulse operational matrix method [30, 31].

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 [3335] 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.

2. Preliminaries and Notations

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

2.1. 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 , , is said to be in the space , if there exists a real number such that , where , and it is said to be in the space as if and only if , .

Definition 2. The Riemann-Liouville fractional integration of order of a function , is defined as [14]

Definition 3. The fractional derivative of in the Caputo sense is defined as [14] for , , , .
Some basic properties of the fractional operator are as follows [14, 40], for , , and we have(1),(2),(3),(4),(5),(6)

where denotes the smallest integer greater than or equal to .

2.2. Chebyshev Polynomials

Chebyshev polynomials of the first kind of degree can be defined as follows: which are orthogonal with respect to the weight function [41], 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 . A function can be written in terms of Chebyshev polynomials as

2.3. 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 , ; , , we will get the following family of discrete wavelet which forms a wavelet basis for :

In particular, when and , then forms an orthonormal basis.

Chebyshev wavelets have four arguments; can assume any positive integer, is degree of Chebyshev polynomials of the first kind, and denotes the time. Consider where and and . In (11) the coefficients are used for orthonormality. We should note that in dealing with the Chebyshev wavelets, the weight function has to be dilated and translated to get orthogonal wavelets as follows:

In view of (5), Chebyshev wavelets are an orthonormal set with respect to the weight function because

Lemma 4. The family of Chebyshev wavelets forms an orthonormal basis for with respect to the weight function [44].
For and , Chebyshev wavelets are as follows:

2.4. Function Approximation

A function defined over may be expanded as where , in which denotes the inner product in . If we consider truncated series in (16), we obtain where and are matrices given by

3. 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 . Using (4) we have, so can be rewritten as

The first two derivatives of the function at the points are given by [46] as where with , for .

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.

4. Chebyshev Wavelet Finite Difference Method

In this section, we present the Chebyshev wavelet finite difference (CWFD) method. Consider , , , as the corresponding Chebyshev-Gauss-Lobatto collocation points at the th subinterval such that

A function can be written in terms of Chebyshev wavelet basis functions as follows: where , , , are the expansion coefficients of the function at the subinterval and , , , are defined in (11).

In view of (11) and (19), we can obtain the coefficients as

Using (25), the first two derivatives of the function at the points , , , can be obtained as where

5. 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].

Theorem 6. A function , 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 have if , by substituting , it yields Using integration by part, we get the first part is zero, therefore,
Using integration by part again, it yields where
Thus, we get
However
Since , we obtain
Now, if , by using (35), we have
It is mentioned in [42] that form an orthogonal system constructed by Haar scaling function with respect to the weight function , so is convergent. Hence, we will have
Therefore, with the aid of Lemma 5, the series converges to uniformly [47].

Theorem 7. Suppose 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 where
This series converges to uniformly. We first show that converges to . We have by substituting , it yields
Using the trapezoidal rule for integration with equidistant subintervals gives
From (26) and (28), for , we have
According to approximation error for the trapezoidal rule, we have where
In view of (35) and triangle inequality, we get where Because , it is understandable that
Therefore, in view of Lemma 5, series (43) is uniformly convergent to .

Theorem 8 (accuracy estimation). Suppose with bounded second derivative, say , then one has the following accuracy estimation: where

Proof . We have
We know that the family forms orthonormal basis for , so . Therefore, in view of (52), we will have where

6. Discretization of Problem

In this section, the Chebyshev wavelet finite difference method (CWFD) is used for solving the following general form: subject to the conditions where , , , are located in and H can be linear or nonlinear while are linear functions.

We suppose the interval is divided into subintervals , . We also consider the shifted Chebyshev-Gauss-Lobatto collocation points on the th subinterval , , where is defined as follows:

In order to obtain the solution in (60), we first approximate according to (27) and rewrite , as fractional derivatives of , in the Caputo sense using (2). Collocating (60) at the shifted Chebyshev-Gauss-Lobatto points , , , we get We then continue as follows:

To calculate the first integral in the above summation, we use the Clenshaw-Curtis quadrature formula [48]: where are Chebyshev-Gauss-Lobatto nodes and the weights are given by where and , for .

In this paper we set . We obtain

Using integration by part and in view of (62), we convert second singular integral to nonsingular one as follows: get the second integral

Replacing (67) and (68) into (64), we obtain , and then replace them into (64) to get equations.

Furthermore, substituting (27) and (29) into (61), we get equations.

Moreover, we should impose continuity condition on the approximate solution and its first derivatives at the interface between subintervals which results in equations.

We will totally have a system of algebraic equations, which can be solved for the . Consequently, we obtain the solution to the given (60) using (27) and (28).

7. 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.

Example 1. As the first example, we consider a nonlinear equation defined as follows [49, 50]: such that the second initial condition is for only. The exact solution of (70) and (71) is given as [26]
It should be noted that for , the slope of the solution at goes to infinity. Therefore, one can expect a large numerical error near .
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 and and . It can be seen that increasing the values of and results in more accurate solution. In Figure 2 numerical results for , , and and (exact solution) and also , and (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 , 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 subject to the following initial states
The exact solution of this problem is .
We solve this fractional initial value problem by applying the method described in Section 6 with . The absolute error is shown in Figure 3, which shows that the numerical solution is in perfect agreement with the exact solution. We set . However, we get the absolute error less or equal to when we set .

Example 3. Following El-Mesiry et al. [52] and Li [33], we consider the following nonlinear fractional differential equation: where subject to
The exact solution of this problem is .
For , , , 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 .

Example 4. Consider the following boundary value problem in the case of the inhomogeneous Bagley-Torvik equation [53, 54] where the exact solution is .
It is worth mentioning that the method presented above only can be employed for solving this problem for . That is because the first kind Chebyshev wavelet is defined on interval . However, the variable in this example is defined on interval , and we should replace Chebyshev wavelet bases with in the discrete procedure. We applied the method described in Section 6 and got almost exact solution for , although other authors had solved the problem for and . It can be seen from Figure 4 that the approximate solution and exact solution are closely overlapped for any .

Example 5. Consider the following boundary value problem for nonlinear fractional order differential equation: where , , , and is a given problem. For , , , , and , it can be easily verified that the exact solution is . We use the introduced scheme in Section 6 to solve this example with , . 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 in Figure 5.

Example 6. Consider the following boundary value problem: where , and . The exact solution of the problem is not known generally. It can be easily verified that for , , the exact solution is . 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 , . In Table 5, the absolute errors are presented which confirm that the proposed method is more accurate. The numerical results for 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 is only for . The exact solution of (74) and (75) is given by where is the Mittag-Leffler function of order . It can be easily verified that for and , the exact solutions are and , respectively.
We solved the problem using the proposed method for different values of .
The absolute errors for (with , ) and (with , ) 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 , and 1 (Figure 7(a)) and , 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 and , respectively.

Example 8. As the last example, consider the following linear multiterm fractional boundary value problem: subject to
The exact solution of this problem is . 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 , while the resulting one by the current method is of order . 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.

8. 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.

Acknowledgment

The authors gratefully acknowledge that this research was partially supported by the University Putra Malaysia under the ERGS Grant Scheme having project no. 5527068.