Fractional differential equations (FDEs) are generalizations of differential equations that replace integral order derivatives by fractional order derivatives. In general, ordinary differential equations are applied on describing dynamic phenomena in various fields such as physics, biology and chemistry. However, for some complicated systems the common simple differential equations cannot provide agreeable results. Therefore, in order to obtain better models, FDEs are employed instead of integer order ones, see [3, 8, 16]. On the other hand, the FDEs are too complicated to solve by analytical methods and theoretical background for this problem is not well developed. Hence, in recent years mathematicians have discovered new methods of numerical solution. There are several methods to solve FDEs, such as variational iteration method [14, 15], Adomain decomposition method [2], fractional differential transformation method [1], fractional finite difference method [12], and wavelet method [9, 17].

Taylor series development can be very straightforward for representing a function as power series, and Taylor series are frequently employed to approximate solutions of complicated problems. There are several studies in which Taylor series methods are applied to solve problems in control theory [10] or in partial differential equations [4]. In the recent years, Taylor series have also been employed to solve FDEs numerically, with the Bagley-Torvik (B-T) equation as an example. However, the B-T equation is an FDE with constant coefficients [7]. Nevertheless, in this paper, Taylor series is the chief tool used to solve FDEs with variable coefficients numerically, by assuming that the solution can be expanded as a Taylor series.

In this paper, we introduce a novel method to approximate the solutions of FDEs with given initial values. In this technique, the solution is approximated by Taylor vectors. Moreover, an operational matrix to integrate the FDEs is provided, and is utilized to project a continuous space to a discrete space, and to form a system of algebraic equations. Some basic definitions and effective theorems of fractional calculus are introduced in Section 1. Section 2 presents an error analysis of our method. Due to the fact that solutions are approximated by Taylor series, a Taylor basis operational matrix of fractional integration will be provided in Section 3. In Section 4, numerical examples demonstrate obtaining approximate solutions. Finally, the Conclusion will be in Section 5.

Due to the fact that there are several definitions of fractional derivatives and integrals, this section is necessary to define the fractional calculus used. The Caputo derivative and the Riemann-Liouville integral are well-known and widely applied. In this paper, we focus on not only the Caputo's fractional derivative, but also on the Riemann-Liouville fractional integral.

Definition 2.1.([13])

The Caputo's fractional derivative of order α is defined as

where α is the order of the derivative and n is the smallest integer which is greater than α.

Definition 2.2.([13])

The Riemann-Liouville fractional integral operator of order α ,I^α, is given by

An important property of the Riemann-Liouville fractional integral that it is a linear operator [5], i.e.

where λ is a constant. Moreover, a well-known Riemann-Liouville fractional integral formula regards its effects on powers of the integrand

Furthermore, Caputo's fractional derivative of order α also obeys

where α∈R and n−1<α≤n with n∈N. The relation between the Caputo fractional derivative and Riemann-Liouville integral operator is that they are "almost" inverses, except for integration constants that necessarily emerge:

Since the purpose of this paper is to find numerical solutions of FDEs, the error bounds need to be considered, which is now addressed. The Taylor basis vector is given by

where m is a positive integer. Clearly Tm⊂H , where H=L²[0,1]. Let S=span {1,t,t2,...,tm} and y be an arbitrary element in H. Since S is a finite dimensional vector subspace of H, there exists a unique y₀∈ S which is the best approximation of y, i.e.

Since y₀∈ S, there exits the unique coefficients a₀,a₁,....,a_m such that

where AT=[a0,a1,a2,...,am].

By Taylor's theorem, see [11], for smooth enough y∈Cm+1[0,1]. Then

where y0(t)=∑ k=0my (k)(0)tkk! and Rm(t)=y(m+1)(c)tm+1(m+1)! for some c∈R.

Lemma 3.1. Let y₀(t) be the best approximation of y∈ S and y(t)∈Cm+1[0,1] then

where M=supt∈[0,1]∥y(m+1)(t)∥.

Proof. The proof is similar to those in [6, 7].

The next theorem is the main theorem in this paper. When a solution of an FDE is approximated by a truncated Taylor series, the error bounds for fractional integrals must be considered.

Theorem 3.2. Suppose that the conditions of Lemma 3.1 are satisfied. Then, for t∈[0,1],

Proof. Since the fractional integral is a linear operator, we have

Substituting inequality (3.4) into (3.5), we can achieve the desired result.

Now we apply Definition 1.2 of the Riemann-Liouville fractional integral I^α. The Taylor basis operational matrix representing the linear operation of fractional order integration Iα is defined by

In this section, we apply the Taylor basis operational matrix of fractional order integration to solve four examples of fractional differential equations with variable coefficients. These problems are based on initial value problems. The examples demonstrate that our method is practically applicable. To obtain better results, we should increase the number of basis functions used in the approximation, or essentially truncate the Taylor series later. For convenience, the programs for computational solutions were written in MatLab.

Example 5.1. Consider the following fractional differential equation

with the initial condition is u(0)=0 and f(t)=32Γ(23)t2/3+t4/3. The exact solution is u(t)=t. This example is provided by [17].

We provide two alternative approaches to approximate solution. First, we assume that

Applying the equation (2.6) with an initial condition, we have

Substituting equations (5.2) and (5.3) into equation (5.1), we obtain

In order to determine A, we collocation at the points ti=t0+ih with h=1/m and t₀=0. For convenience, the coefficient in (5.1) can be set as

For collocation the right-hand side of (5.1) is similarly discretized to

Now the continuous problem has been projected to discrete space. Hence, we have algebraic equations with At to be solved from

The eventual numerical solution isu(t)=AtI1/3Tm(t). The numerical implementations for m = 5, 10, and 15 are illustrated by their absolute errors in Fig.1(a), since the analytical solution in this case is known. In Fig.1(a), we can see that the numerical solutions had good accuracy. Furthermore, Table 1 compares these absolute errors to those of the numerical solution in [17]. Moreover, the absolute errors consistently diminished with the number of basis functions used.

Table 1 . The absolute error in Example 5.1 shown in comparison to the method of [17]..

t	Ours (m=15)	Ref.[17](m=64)
0.0625	8.745028e-004	2.082093e-003
0.1875	2.498965e-004	9.743606e-004
0.3125	1.480609e-004	6.67794e-004
0.4375	1.037922e-004	5.129034e-004
0.5625	7.883851e-005	4.175219e-004
0.6875	6.289737e-005	3.523496e-004
0.8125	5.238139e-005	3.048147e-004
0.9375	3.629887e-005	2.685403e-004

Figure 1. This absolute error in Example 5.1 with various numbers m of Taylor basis functions. (a)
u ( t ) = A t I 1 / 3 T m ( t ) , (b)
u ( t ) = A t I 1 T m ( t ) . }

An alternative choice to solve this case is by setting

Then, we obtain

The approach to determine A is essentially unchanged, by employing collocation. In order to avoid a singular matrix, the collocation points are set at t_i=t₁+ih with h=1/m and t₁=h for i=1,2,3,...,m. Then the coefficient matrix of (5.1) is

and the right-hand side of (5.1) can be restated as

The system of an algebraic equation is analogous to that in (5.5). Finally, the solution is approximated by u(t)=AtI1Tm(t). With this latter approach, the approximate solution is very close to the exact solution. The absolute errors are shown in Fig 1(b). One can see that the errors are better than with the previous approach. For good results, we should set Dnu(t)=AtTm(t) with n=α where α is the order of the fractional differential equation.

The next example is more complicated than the previous, since the initial values are nonzero. Still the method of solution remains quite similar.

Example 5.2. Consider the following initial value problem

where f(t)=2+2cos(t)t2/3Γ(3/2)−t3−t with initial condition u(0)=1 and u'(0)=0. The exact solution is u(t)=t²⁺¹.

The method to solve this problem is similar to Example 5.1. First, we set

Then the initial conditions are applied to get

Applying the collocation points defined as Example 5.1, we let

Substituting the equation (5.10)-(5.16) into (5.9), we obtain a system of algebraic equation

This system of algebraic equations is solved for vector A. Then, employing (5.12), we get the numerical approximation of u(t). The absolute errors with m = 5, 10 and 15 are exhibited in Fig.2. We can see that the absolute errors are very small, so this example corroborates that our method can work very well.

Figure 2. This absolute errors of Example 5.2 for various numbers m of the Taylor basis functions.

Example 5.3.([17]) Consider the following initial value problem

where f(t)=t9−t8+56t6−42t5+sin(t)(327686435t15/2−2048429t13/2) with initial condition u(0)=0 and u'(0)=0. The exact solution is u(t)=t8−t7.

The method to solve this problem is similar to Example 5.1 . First, we set

Then the initial conditions are applied, we have

Applying the collocation points defined as Example 5.1, we let

align the equation (5.19)-(5.24) into (5.18), we obtain a system of algebraic equation

This is solved for vector A, and employing (5.21) we obtain the numerical approximation of u(t). The absolute errors with m = 5, 10 and 15 are shown in Fig.3. Moreover, Table 2 demonstrates the absolute errors in comparison to those in [17]. We can see that the absolute errors are very small, confirming that our method can work very well.

Table 2 . The absolute error in Example 5.3 shown in comparison to the method of [17]..

t	Ours (m=15)	Ref.[17](m=32)
0.0625	8.9369e-013	4.4517e-012
0.1875	2.6704e-009	8.3173e-008
0.3125	2.9526e-007	1.5545e-006
0.4375	6.0751e-006	9.6917e-006
0.5625	5.3743e-005	4.4897e-005
0.6875	2.7960e-004	1.8488e-004
0.8125	9.7785e-004	5.2730e-004
0.9375	2.3574e-003	8.3468e-004

Figure 3. This absolute errors of Example 5.3 for various numbers m of the Taylor basis functions.

Example 5.4. Consider the following initial value problem

with initial condition u(0)=0.

The exact solution of this problem is in general unknown, but for α=1 it is u(t)=t^2. The process to obtain an approximate solution is the same in Examples 5.1-5.3. Fig. 4 demonstrates that the approximate solutions appear to converge to the exact one, for the case α=1. For more detail see [5].

Figure 4. A comparison of the approximate solutions for different α values, with 15 Taylor basis functions.

The objective of this paper is to demonstrate numerical solutions of fractional differential equations with variable coefficients. The technique employs Taylor series approximations. Moreover, a convergence analysis with the Taylor basis was also proved. The operational matrix of fractional integration for a Taylor vector was provided. Furthermore, this method converts initial value problems into linear systems of algebraic equations. The method is computationally very easy and provides a structured approach to numerical approximate solutions. Numerical implementations were presented illustrating accuracy, application and efficiency of the approach.

The authors greatly acknowledge the valuable comment of Prof. Mohsen Razzaghi (Mississippi State University) and Assoc. Prof. Varayu Boonpongkrong (Prince of Songkla University) to obtain better result. Moreover, the author also gratefully acknowledge the copy-editing service of the Research and Development Office (RDO/PSU), and the helpful edits and comments by Dr. Seppo Karrila.