Indicators

• Cited by SciELO
• Access statistics

Print version ISSN 2238-3603On-line version ISSN 1807-0302

Comput. Appl. Math. vol.23 no.1 Petrópolis  2004

http://dx.doi.org/10.1590/S0101-82052004000100002

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

A. M. A. El-SayedI; A. E. M. El-MesiryII; H. A. A. El-SakaII

IFaculty of Science, Alexandria University, Alexandria, Egypt, E-mail: amasayed@maktoob.com
IIDepartment of Mathematics, Damietta Faculty of Science, New Damietta, Egypt, E-mail: elmisiery@hotmail.com/ halaelsaka@hotmail.com

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 ,  and ). 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 , ,  and ). 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  and -).

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

Consider the initial value problem The existence of at least one solution of this initial value problem has been proved (see ) where the function f(t,U) = f(t,u1(t),...,un(t)) satisfies Caratheodory conditions, i.e., t ® f(t,U) is measurable for every U Î Rn+1 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) Î L1 and constants bk > 0, such that

|f(t,U)| < a(t)+ bk|uk(t)| for all (t,U) Î I×Rn+1.

In this paper we focus on providing a numerical solution to the nonlinear multi-term fractional (arbitrary) orders differential equation (see ). subject to the initial values where ai are real numbers (i = 1,2,...,m), such that

0 < a1 < a2 < ... < am < 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 , the Bagley-Torvik equation  and the Basset equation .

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

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

Definition 1.1. Let f(t) Î L1,b Î R+. The fractional (arbitrary) order integral of the function f(t) of order b is defined by (see , ,  and ) when a = 0 we can write Ibf(t) = f(t) = f(t)*fb(t), where fb(t) = for t > 0,fb(t) = 0 for t < 0 and fb® d(t) (the delta function) as b ® 0 (see ).

Definition 1.2. The fractional derivative Da of order a Î (0,1] of the absolutely continuous function g(t) is defined as (see ,, ,  and ) 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 Ai+1(t)xi+1(t) = (t),i = 0,1,2,...,m-1,a0 = 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,x0(t),x1(t),...,xm(t)) satisfies the Lipschitz condition Theorem 2.1. Let f(t,x0,x1,...,xm) Î 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,x0,x1,...,xm) and (t,y0,y1,...,ym) Î , we get Since for t > 0 and i = 0,1,2,...,m-1 further, < ai+1-ai+1 < ai+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 , ,  and ).

In  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 ). That may require a large amount of computational effort to calculate its weights.

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

In  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  and ). The algorithm is used for solution of the linear and non-linear equations.

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

In  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  and  and investigated further in -. 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

a0 = 0, x0(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  and . 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 x3(t) is introduced. This problem was solved in . 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  and . 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  and .

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 . 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 ), subject to

x(0) = 0.

Where a = ba,b = 9/(1+2l) and l = 0.5,2,10 and 100. Also this problem was solved in . 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) = t3.

For a = 1,b = 1,c = 1,e = 1,a1 = 0.555 and a2 = 1.455 we found the following results (see Table 3). For a = 1,b = 0.5,c = 0.5,e = 0.5,a1 = 0.276 and a2 = 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,a1 = 0.821 and a2 = 1.333 we found the following results (see Table 5). For a = 4,b = 3,c = 2,e = 1,a1 = 0.637 and a2 = 1.211 we found the following results (see Table 6). Example 8. Consider the equation and

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

subject to

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

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

x(t) = 1+ t2.

Set

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

then we have the following equation where

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

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,a1 = 0.0159 and a2 = 0.1379 we found the following results (see Table 7). For a = 0.2,b = 1,c = 1,e = 0.5,k = 2,a1 = 0.00196 and a2 = 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- t2.

Set

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

then we have the following equation where

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

subject to

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

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

For a = 1,b(t) = t1/2,c(t) = t1/3,e(t) = t1/4,k(t) = t1/5,a1 = 0.333 and a2 = 1.234 we found the following results (see Table 9). For a = 3,b(t) = t, c(t) = t+1,e(t) = t2,k(t) = (t+1)2,a1 = 0.445 and a2 = 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) = ,a1 = 0.999 and a2 = 1.222 we found the following results (see Table 11). For a = 0.5,b(t) = t2/5,c(t) = t2-t,e(t) = t3-t2,k(t) = t+1,a1 = 0.123 and a2 = 1.888 we found the following results (see Table 12). 5 Conclusion

In our method a1,a2,...,am take arbitrary values such that 0 < a1 < a2 < ... < am < n. On the other hand a specific conditions must be satisfied by as (see ,  and ).

Examples 6-10 have been solved (see ,  and ) 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, a1,a2, ..., am,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 a1 = 0.55 and a2 = 1.45 the order of the system becomes 0.05 and the dimension of the system becomes 40 further after approximations a1 = 0.5 and a2 = 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.

REFERENCES

 L. Blank, Numerical treatment of differential equations of fractional order, Numerical Analysis Report 287. Manchester Center for Numerical Computational Mathematics (1996).         [ Links ]

 M. Caputo, Linear model of dissipation whose Q is almost frequency independent II, Geophys. J.R. Astr. Soc. 13 (1967), 529-539.         [ Links ]

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

 K. Diethelm and G. Walz, Numerical solution of fractional order differential equations by extrapolation, Numerical Algrithms 16 (1997), 231-253.         [ Links ]

 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).         [ Links ]

 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).         [ Links ]

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

 K. Diethelm, Predictor-Corrector strategies for single and multi-term fractional differential equations, Proceedings Hercma 2001.         [ Links ]

 K. Diethelm, N.J. Ford and A.D. Freed, Detailed error analysis for a fractional Adams method , Numer. Algorithms, to appear.         [ Links ]

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

 K. Diethelm and N.J. Ford, Numerical solution of the Bagley-Torvik equation BIT, 42 (2002), 490-507.         [ Links ]

 K. Diethelm and Y. Luchko, Numerical solution of linear multi-term differential equations of fractional order, J. Comput. Anal. Appl., (2003), to appear.         [ Links ]

 A.M.A. El-Sayed, Nonlinear functional differential equations of arbitrary orders, Nonlinear Analysis: Theory, Methods and Applications, 33 (2) (1998), 181-186.         [ Links ]

 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).         [ Links ]

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

 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).         [ Links ]

 I.M. Gelfand and G.E. Shilov, Generalized Functions, 1 (1958).         [ Links ]

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

 C. Lubich, Discretized fractional calculus, SIAM Journal of Mathematical Analysis, 17 (3) (1986), 704-719.         [ Links ]

 J. Leszczynski and M. Ciesielski, A numerical method for solution of ordinary differential equations of fractional order, V 1 26 Feb 2002.         [ Links ]

 K.S. Miller and B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley & Sons, Inc., New York (1993).         [ Links ]

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

 K.B. Oldham and J. Spanier, The fractional Calculus, Academic Press, New York and London (1974).         [ Links ]

 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).         [ Links ]

 A. Palczewski, Ordinary differential equations, WNT, Warsaw (1999) (in Polish).         [ Links ]

 I. Podlubny, Fractional differential equations, Academic Press (1999).         [ Links ]

 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).         [ Links ]

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

#580/03. All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License