Numerical solution for multi-term fractional (arbitrary) orders differential equations

Numerical solution for multi-term fractional (arbitrary) orders differential equations

A. M. A. El-Sayed^I; A. E. M. El-Mesiry^II; H. A. A. El-Saka^II

^IFaculty of Science, Alexandria University, Alexandria, Egypt, E-mail: [e-mail]

^IIDepartment of Mathematics, Damietta Faculty of Science, New Damietta, Egypt, E-mail: [e-mail]/ [e-mail]

ABSTRACT

Our main concern here is to give a numerical scheme to solve a nonlinear multi-term fractional (arbitrary) orders differential equation. Some results concerning the existence and uniqueness have been also obtained.

Mathematical subject classification: 34A99.

Key words: multi-term fractional differential equation, Caputo derivative, existence, uniqueness, predictor-corrector method.

1 Introduction

The use of fractional orders differential and integral operators in mathematical models has become increasingly widespread in recent years (see [11], [22] and [26]). Several forms of fractional differential equations have been proposed in standard models, and there has been significant interest in developing numerical schemes for their solution (see [11], [14], [22] and [26]). However, much of the work published to date has been concerned with linear single term equations and, of these, equations of order less than unity have been most often investigated (see [1] and [3]-[6]).

Let a Î (n,n+1],a_kÎ (k-1,k],k = 1,2,...,n and a_o = 0.

Consider the initial value problem

The existence of at least one solution of this initial value problem has been proved (see [13]) where the function f(t,U) = f(t,u₁(t),...,u_n(t)) satisfies Caratheodory conditions, i.e., t ® f(t,U) is measurable for every U Î Rⁿ⁺¹ and U ® f(t,U) is continuous for every t Î I. f(t,U) is nondecreasing for all variables, and there exist a function a(t) Î L₁ and constants b_k > 0, such that

|f(t,U)| < a(t)+ b_k|u_k(t)| for all (t,U) Î I×Rⁿ⁺¹.

In this paper we focus on providing a numerical solution to the nonlinear multi-term fractional (arbitrary) orders differential equation (see [15]).

subject to the initial values

where a_i are real numbers (i = 1,2,...,m), such that

0 < a₁ < a₂ < ... < a_m < n

and n is any positive integer number.

A theorem proving the existence and uniqueness of the solution will be proved.

Applications for such equations arise, e.g., in various areas of mechanics [26], the Bagley-Torvik equation [11] and the Basset equation [22].

Now we give the definition and some properties of the fractional order differential and integral operators.

Let L₁ = L₁[a, b] be the class of Lebesgue integrable functions on [a, b],a < b < ¥.

Definition 1.1. Let f(t) Î L₁,b Î R⁺. The fractional (arbitrary) order integral of the function f(t) of order b is defined by (see [18], [21], [24] and [27])

when a = 0 we can write I^bf(t) = f(t) = f(t)*f_b(t), where f_b(t) = for t > 0,f_b(t) = 0 for t < 0 and f_b® d(t) (the delta function) as b ® 0 (see [16]).

Definition 1.2. The fractional derivative D^a of order a Î (0,1] of the absolutely continuous function g(t) is defined as (see [2],[13], [18], [24] and [27])

2 Formulation of the problem

Definition 2.1. By a solution of the initial value problem (1.1) and (1.2) we mean a function x(t) Î C(I) and all its derivative up to order (n-1) are vanishing at t = 0.

Now equation (1.1) can be written in the form

where

and

From Eq. (2.2) we get

and

Lemma 2.1. The fractional order differential equation (2.1) can be transformed to the system

where

and A_i+1(t)x_i+1(t) = (t),i = 0,1,2,...,m-1,a₀ = 0 where (') denotes the transpose of the matrix.

Proof. Using the definition of the fractional derivative and equations (2.3), (2.4) and (2.5), we can easily obtain the result.

Set I = [0,T], say, T is a suitable positive number and = I×C^*(t) where C^*(t) is the class of all continuous column vectors X(t) (defined by (2.7)) with the norm

Definition 2.2. By a solution of the system (2.6), we mean a column vector X(t) Î C^*(t) and X(0) = 0.

Assuming that the function f(t,x₀(t),x₁(t),...,x_m(t)) satisfies the Lipschitz condition

Theorem 2.1. Let f(t,x₀,x₁,...,x_m) Î C() and satisfies the Lipschitz condition (2.10). If

then Eq. (1.1) has one and only one solution x(t) Î C(I) that satisfies x(t) Î C(I), i = 1,2,...,m.

Proof. If we write

Then for (t,x₀,x₁,...,x_m) and (t,y₀,y₁,...,y_m) Î , we get

Since for t > 0 and i = 0,1,2,...,m-1

further, < a_i+1-a_i+1 < a_i+1+1 < n+1 hence

since

then

Hence the map F:C^*(I) ® C^*(I) is a contraction (and then, it has a fixed point X = FX) providing

and hence, there exists a unique column vector X(t) Î C^*(I), which is the solution of the system (2.6). Therefore, from (2.7) and the definition of C^*(I), we deduce that there exists one and only one solution x(t) Î C(I), and this solution satisfies x(t) Î C(I),i = 1,2,...,m.

3 Numerical methods and results

Numerical methods for the solution of linear fractional differential equations involving only one fractional derivative are well established (see for example [1], [3], [4] and [16]).

In [12] Luchko and Diethelm discussed a new algorithm for the numerical solution of initial value problems for general linear multi-term differential equations of fractional order with constant coefficients and fractional derivatives in the Caputo sense. Which is obtained by applying the convolution quadrature and discretized operational calculus to the analytical solution of the problem given in terms of the Mittag-Leffer type function (see [19]). That may require a large amount of computational effort to calculate its weights.

In [14] Ford et al. showed how the numerical approximation of the solution of a linear multi-term fractional differential equation can be calculated by reduction of the problem to a system of ordinary and fractional differential equations each of order at most unity (see [3]).

In [20] Leszczynski and Ciesielski proposed an algorithm for the numerical solution of arbitrary differential equations of fractional order. Which is obtained by using the following decomposition of the differential equation into a system of differential equation of integer order connected with inverse forms of Abel-integral equations (see [23] and [25]). The algorithm is used for solution of the linear and non-linear equations.

In [11] Diethelm and Ford considered the reformulation of the Bagley-Torvik equation as a system of fractional differential equations of order 1/2.

In [15] El-Sayed et al. gave a numerical methods for solving a nonnlinear multi-term fractional (arbitrary) orders differential equations.

An Adams-type predictor-corrector method has been introduced in [5] and [6] and investigated further in [7]-[10]. In this paper we use an Adams-type predictor-corrector method for the numerical solution of fractional integral equations. We will apply the PECE (Predict, Evaluate, Correct, Evaluate) method in Eqs.

where

a₀ = 0, x₀(t) = x(t)

and

4 Numerical examples

As an example that arises in application, we solve the Bagley-Torvik equation which arises, for instance, in modelling the motion of a rigid plate immersed in a Newtonian fluid. Examples 1-3 are The Bagley-Torvik equation.

Example 1.

where

subject to

x(0) = x'(0) = 0.

This problem was solved in [26] and [20]. The approximate solutions displayed in Figs. 1 and 2 for step size h = 0.01.

Example 2. Consider a non-linear form of the fractional differential equation

where

subject to

x(0) = x'(0) = 0.

This is a Bagley-Torvik equation where non-linear term x³(t) is introduced. This problem was solved in [20]. Fig. 3. shows a behaviour of the numerical solution for step size h = 0.01.

Example 3.

where

f(t) = c(t+1)

subject to

x(0) = x'(0) = 1.

It is easily verified that the exact solution of this problem is

x(t) = t+1.

This problem was solved in [11] and [12]. Set

y(t) = x(t)-t-1,

then we have the following equation

subject to

y(0) = y' (0) = 0.

Now, we can apply our method in Eq. (4.4). We evaluate the maximal error (m. e.) and the results are compared to results obtained by the methods of [11] and [12].

For a = 1,b = 0.5 and c = 0.5 we found the following results (see Table 1).

For a = 1,b = 1 and c = 1 we found the following results (see Table 2).

Example 4. The equation

subject to

x(0) = 1,x'(0) = 0.

This problem was solved in [14]. Set

y(t) = x(t)-1,

then we have the following equation

subject to

y(0) = y'(0) = 0.

Now, we can apply our method in Eq. (4.6). We show the approximate solutions in Figs. 4-8 for step size h = 0.01.

Example 5. Consider the linear fractional differential equation (4.7) which arises, for instance, in the study of the generalised Basset force occuring when a sphere sinks in a (relatively less dense) viscous fluid (see Mainardi [22]),

subject to

x(0) = 0.

Where a = b^a,b = 9/(1+2l) and l = 0.5,2,10 and 100. Also this problem was solved in [14]. Our calculated approximate solutions (with h = 0.1) to the Basset equation are given in Figs. 9 and 10.

Example 6. Consider the nonlinear equation

and

subject to

x(0) = x'(0) = 0.

It is easily verified that the exact solution of this problem is

x(t) = t³.

For a = 1,b = 1,c = 1,e = 1,a₁ = 0.555 and a₂ = 1.455 we found the following results (see Table 3).

For a = 1,b = 0.5,c = 0.5,e = 0.5,a₁ = 0.276 and a₂ = 1.999 we found the following results (see Table 4).

Example 7. Consider the nonlinear equation

and

f(t) = +et

subject to

x(0) = 0, x'(0) = 1.

It is easily verified that the exact solution of this problem is

x(t) = t.

Set

y(t) = x(t)-t,

then we have the following equation

subject to

y(0) = y'(0) = 0

Now, we can apply our method in Eq. (4.10).

For a = 1,b = 1,c = 0.5,e = 0.5,a₁ = 0.821 and a₂ = 1.333 we found the following results (see Table 5).

For a = 4,b = 3,c = 2,e = 1,a₁ = 0.637 and a₂ = 1.211 we found the following results (see Table 6).

Example 8. Consider the equation

and

f(t) = a + bt + + + k (1 + t²)

subject to

x(0) = 1,x'(0) = 0.

It is easily verified that the exact solution of this problem is

x(t) = 1+ t².

Set

y(t) = x(t)-1,

then we have the following equation

where

g(t) = a+bt+ + + t²

subject to

y(0) = y'(0) = 0.

Now, we can apply our method in Eq. (4.12).

For a = 1,b = 3,c = 2,e = 1,k = 5,a₁ = 0.0159 and a₂ = 0.1379 we found the following results (see Table 7).

For a = 0.2,b = 1,c = 1,e = 0.5,k = 2,a₁ = 0.00196 and a₂ = 0.07621 we found the following results (see Table 8).

Example 9. Consider the equation

and

subject to

x(0) = 2,x'(0) = 0.

It is easily verified that the exact solution of this problem is

x(t) = 2-t².

Set

y(t) = x(t)-2,

then we have the following equation

where

g(t) = -a--c(t)t--t²

subject to

y(0) = y'(0) = 0.

Now, we can apply our method in Eq. (4.14).

For a = 1,b(t) = t^1/2,c(t) = t^1/3,e(t) = t^1/4,k(t) = t^1/5,a₁ = 0.333 and a₂ = 1.234 we found the following results (see Table 9).

For a = 3,b(t) = t, c(t) = t+1,e(t) = t²,k(t) = (t+1)²,a₁ = 0.445 and a₂ = 1.178 we found the following results (see Table 10).

Example 10. Consider the equation

and

subject to

x(0) = 0,x'(0) = 1.

It is easily verified that the exact solution of this problem is

x(t) = t.

Set

y(t) = x(t)-t,

then we have the following equation

subject to

y(0) = y'(0) = 0.

Now, we can apply our method in Eq. (4.16).

For a = 2,b(t) = t,c(t) = t,e(t) = ,k(t) = ,a₁ = 0.999 and a₂ = 1.222 we found the following results (see Table 11).

For a = 0.5,b(t) = t^2/5,c(t) = t²-t,e(t) = t³-t²,k(t) = t+1,a₁ = 0.123 and a₂ = 1.888 we found the following results (see Table 12).

5 Conclusion

In our method a₁,a₂,...,a_m take arbitrary values such that 0 < a₁ < a₂ < ... < a_m < n. On the other hand a specific conditions must be satisfied by as (see [7], [8] and [10]).

Examples 6-10 have been solved (see [7], [8] and [10]) by transforming the given initial value problem into a system of equations, all of which must have the same order and the dimension of the system is d = n/q (q = gcd(1, a₁,a₂, ..., a_m,n)) which in most cases is very large because the order of the system q is small and it is obvious that this leads to much larger requirements concerning computer memory and run time.

For instance in Example 6 the order of the system q = 0.005 and the dimension of the system d = 400 and after approximations a₁ = 0.55 and a₂ = 1.45 the order of the system becomes 0.05 and the dimension of the system becomes 40 further after approximations a₁ = 0.5 and a₂ = 1.5 the order of the system becomes 0.5 and the dimension of the system becomes 4.

This produce two sources of errors one due to the approximations of as to reduce the dimension of the system and the other in the numerical solution of the resulting system.

Received: 28/IV/03. Accepted: 11/III/04.

#580/03.

• [1] L. Blank, Numerical treatment of differential equations of fractional order, Numerical Analysis Report 287. Manchester Center for Numerical Computational Mathematics (1996).
• [2] M. Caputo, Linear model of dissipation whose Q is almost frequency independent II, Geophys. J.R. Astr. Soc. 13 (1967), 529-539.
• [3] K. Diethelm, An algorithm for the numerical solution of differential equations of fractional order, Electronic Transactions on Numerical Analysis, Kent State University, 5 (1997), 1-6.
• [4] K. Diethelm and G. Walz, Numerical solution of fractional order differential equations by extrapolation, Numerical Algrithms 16 (1997), 231-253.
• [5] K. Diethelm and A. Freed, On the solution of nonlinear fractional order differential equations used in the modelling of viscoplasticity, Scientific Computing in Chemical Engineering II-Computational Fluid Dynamics, Reaction Engineering, and Molecular properties, F. Keil, W. Mackens, H. Voß and J. Werther (eds), 217-224, Springer, Heidelberg (1999).
• [6] K. Diethelm and A. Freed, The FracPECE subroutine for the numerical solution of differential equations of fractional order, Forschung und wissenschaftliches Rechnen 1998, S. Heinzel and T. Plesser (eds.), Gesellschaft für Wisseschaftliche Datenverarbeitung, Göttingen, 57-71, (1999).
• [7] K. Diethelm and N.J. Ford, The numerical solution of linear and non-linear Fractional differential equations involving Fractional derivatives several of several orders, Numerical Analysis Report 379, Manchester Center for Numerical Computational Mathematics.
• [8] K. Diethelm, Predictor-Corrector strategies for single and multi-term fractional differential equations, Proceedings Hercma 2001.
• [9] K. Diethelm, N.J. Ford and A.D. Freed, Detailed error analysis for a fractional Adams method , Numer. Algorithms, to appear.
• [10] K. Diethelm, N.J. Ford and A.D. Freed, A predictor-corrector approach for the numerical solution of fractional differential equations, Nonlinear Dynamics, 29 (2002), 3-22.
• [11] K. Diethelm and N.J. Ford, Numerical solution of the Bagley-Torvik equation BIT, 42 (2002), 490-507.
• [12] K. Diethelm and Y. Luchko, Numerical solution of linear multi-term differential equations of fractional order, J. Comput. Anal. Appl., (2003), to appear.
• [13] A.M.A. El-Sayed, Nonlinear functional differential equations of arbitrary orders, Nonlinear Analysis: Theory, Methods and Applications, 33 (2) (1998), 181-186.
• [14] J.T. Edwards, N.J. Ford and A.C. Simpson, The Numerical solution of linear multi-term fractional differential equations: systems of equations, Manchester Center for Numerical Computational Mathematics (2002).
• [15] A.E.M. El-Mesiry, A.M.A. El-Sayed and H.A.A. El-Saka, Numerical methods for multi-term fractional (arbitrary) orders differential equations, Appl. Math. and Comput., (2004), to appear.
• [16] N.J. Ford and A.C. Simpson, The approximate solution of fractional differential equations of order greater than 1, Numerical Analysis Reoprt 386. Manchester Center for Numerical Computational Mathematics (2001).
• [17] I.M. Gelfand and G.E. Shilov, Generalized Functions, 1 (1958).
• [18] R. Gorenflo and F. Mainardi, Fractional Calculus: Integral and Differential Equations of Fractional Order, in A. Carpinteri and F. Mainardi (Eds), Fractals and Fractional Calculus in Continuum Mechanics, Springer, 223-276, (1997), Wien.
• [19] C. Lubich, Discretized fractional calculus, SIAM Journal of Mathematical Analysis, 17 (3) (1986), 704-719.
• [20] J. Leszczynski and M. Ciesielski, A numerical method for solution of ordinary differential equations of fractional order, V 1 26 Feb 2002.
• [21] K.S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley & Sons, Inc., New York (1993).
• [22] F. Mainardi, Some basic problems in continuum and statistical mechanics, in A. Carpinteri and F. Mainardi (Eds), Fractals and Fractional Calculus in Continuum Mechanics, Springer, 291-348, (1997), Wien.
• [23] K.B. Oldham and J. Spanier, The fractional Calculus, Academic Press, New York and London (1974).
• [24] I. Podlubny and A.M.A. El-Sayed, On two definitions of fractional calculus, Solvak Academy of science-institute of experimental phys. UEF-03-96 ISBN 80-7099-252-2. (1996).
• [25] A. Palczewski, Ordinary differential equations, WNT, Warsaw (1999) (in Polish).
• [26] I. Podlubny, Fractional differential equations, Academic Press (1999).
• [27] S.G. Samko, A.A. Kilbas and O.I. Marichev, Integrals and Derivatives of the Fractional Order and Some of Their Applications, Nauka i Tehnika, Minsk, (in Russian) (1987), (English trans. 1993).

Publication Dates

• Publication in this collection
  26 Nov 2004
• Date of issue
  2004

History

• Accepted
  11 Mar 2004
• Received
  28 Apr 2003