# An algorithm for the numerical solution of differential equations of fractional order

Electronic Transactions on Numerical Analysis. Volume 5, pp. 1-6, March 1997. Copyright ? 1997, Kent State University. ISSN 1068-9613.

ETNA

Kent State University etna@mcs.kent.edu

AN ALGORITHM FOR THE NUMERICAL SOLUTION OF DIFFERENTIAL EQUATIONS OF FRACTIONAL ORDER

KAI DIETHELMy Abstract. Differential equations involving derivatives of non-integer order have shown to be adequate models for various physical phenomena in areas like damping laws, diffusion processes, etc. A small number of algorithms for the numerical solution of these equations has been suggested, but mainly without any error estimates. In this paper, we propose an implicit algorithm for the approximate solution of an important class of these equations. The algorithm is based on a quadrature formula approach. Error estimates and numerical examples are given. Key words. Fractional derivative, Riemann-Liouville derivative, differential equation, numerical solution, quadrature formula, implicit method. AMS subject classi?cations. 26A33, 65L70, 65L05.

1. Introduction and Main Results. 1.1. The differential equation. In this paper, we discuss a numerical method for the solution of the fractional differential equation (1.1) (1.2)

?

Dq x ? x0] (t) = x(t) + f (t); x(0) = x0;

0

t

1;

where 0 < q < 1, f is a given function on the interval 0; 1], 0, and x is the unknown function. Here, Dq x denotes the Riemann-Liouville fractional derivative of order q of the function x, de?ned by [6] (1.3)

(

dZt Dq x)(t) := Γ(1 1 q) dt (tx(u))q du: ? ?u 0

Following the common practice in the theory of these differential equations [1], we have incorporated the initial condition (1.2) into the differential equation (1.1). Existence and uniqueness of the solution have been shown in [5]. Since the complete initial value problem may easily be transformed to an arbitrary interval, our choice of the interval 0; 1] does not mean an essential restriction. For q = 1=2, such an equation describes, e. g., the behaviour of a damping model in mechanics [4] where x denotes the displacement, = ?1= ( is the viscosity), and f (t) = lN (t)=(EA ). Here, l is the length of the object under consideration, A is its volume, E is Young’s modulus, and N (t) the external force. From [6, pp. 201ff.], we may see that equation (1.1) can also be used to describe diffusion processes. Other applications are in areas like electromagnetics, electrochemistry, material science, the theory of ultra-slow processes, and special functions, see [1] and the references cited therein. 1.2. The approximation method. A numerical method for the solution of this equation has been proposed in [1]. The method is based on collocation using spline functions, but no error analysis is given.

Received November 28, 1996. Accepted for publication February 5, 1997. Communicated by R. S. Varga. Institut f¨ r Mathematik, Universit¨ t Hildesheim, Marienburger Platz 22, D-31141 Hildesheim, Germany u a (diethelm@informatik.uni-hildesheim.de).

y

1

ETNA

Kent State University etna@mcs.kent.edu

2

K. Diethelm

Our algorithm is based on the observation [3] that we may interchange differentiation and integration in (1.3) to obtain (1.4)

(

1 Dq x)(t) = Γ(?q)

Z

t

0 (

x(u) t ? u)q+1 du;

where now the integral must be interpreted as a Hadamard ?nite-part integral. Then, for a given n, we introduce an equispaced grid tj = j=n on the interval where the solution of eq. (1.1) is sought. Discretizing with this grid and applying (1.4), we obtain for j = 1; 2; : : :; n

Z t 1 x(u) ? x(0) du f (tj ) + x(tj ) = Γ(?q) (tj ? u)q +1 0 Z 1 t?q x(tj ? tj w) ? x(0) dw: j = Γ(?q) 0 wq+1

j

Now, for every j , we replace the integral by a ?rst-degree compound quadrature formula with the equispaced nodes 0; 1=j; 2=j; : ::; 1,

Qj g] :=

with remainder term

j X k=0

kj g(k=j )

Z

0

1

g(u)u?q?1 du

Rj g ] =

Z

0

1

g(u)u?q?1 du ? Qj g]

as proposed in [2]. Since the quadrature formula uses both end points of the integration interval as nodes, we obtain an implicit scheme. Explicit expressions for the weights kj are given in Lemma 2.1 below. Ignoring the quadrature error, we may solve the resulting equation for the values xj which will be our approximations for x(tj ) (j = 1; 2; : : :; n). We obtain the following formulas: (1.5)

xj =

0

j ? (j=n)q Γ(?q)

1

! j j q Γ(?q)f (t ) ? X x ? 1 x : j kj j ?k q 0 n k=1

Here, it is evident that, in contrast to the usual integration methods for differential equations with integer-order derivatives, we cannot say that the method is an m-step method for a certain ?xed m (i. e. that the approximation xj can be determined solely on the basis of the m previous approximations xj ?m; xj ?m+1 ; : : :xj ?1 ). Instead, we observe that every xj depends on all the previous values x0 ; x1 ; : : :; xj ?1. This re?ects the fact that, unlike the classical derivatives of integer order, fractional differential operators are not local operators (i. e. we cannot determine (Dq g)(z ) using only values of g in a neighbourhood of z ). Of course, this also means that our error analysis needs different methods compared to the error analysis in the classical case. Our main results are the following local and global error bounds. THEOREM 1.1. Assuming that the functions involved are suf?ciently smooth, there exists a constant depending on q and x (and therefore on f and ) such that the error of the approximation method described above is bounded by

jx(tj ) ? xj j

j q n?2;

j = 0; 1; : : :; n:

ETNA

Kent State University etna@mcs.kent.edu

Fractional differential equations

3

COROLLARY 1.2. If the functions involved are suf?ciently smooth, we have the following global error estimate for the approximation method described above:

j =0;1;:::;n

max

jx(tj ) ? xj j = O(nq?2 ):

The proof will be given in x 2. The scheme has been tested on some numerical examples. The results are reported in x 3. It is easily seen that the algorithm may be generalized to handle equations of the form

?

Dq x ? x0] (t) = (t)x(t) + f (t)

?

with non-constant . It may also be combined with an explicit scheme to form a predictorcorrector method for the more general nonlinear equation

Dq x ? x0 ] (t) = g(t; x(t)):

However, since the equation stated in (1.1) seems to be the most important case as far as applications are concerned, we shall not go into details about these two generalizations here. 2. Proofs. 2.1. Preliminaries. Before we come to the proofs of the main results, we state some auxiliary lemmas. LEMMA 2.1. For the weights kj of the quadrature formula Qj , j 1, we have

q(1 ? q)j ?q kj = : 2k1?q ? (k ? 1)1?q ? (k + 1)1?q (q ? 1)k ?q ? (k ? 1)1?q + k 1?q

8 <

?1

for k = 0, for k = 1; 2; : : :; j ? 1, for k = j .

Proof. This follows after a simple calculation from the de?nition of the quadrature formula. The following result is taken directly from [2, Theorem 2.3 and the remark following its proof]. LEMMA 2.2. Let q 2 (0; 1). (i) There exists a constant q > 0 such that, for every f 2 C 2 0; 1],

Z

0 1

f (t)t?q?1 dt ? Qj f ]

Z

0 1

q j q?2 kf 00 k1 :

(ii) If f is convex, then

f (t)t?q?1 dt Qj f ]:

Upper and lower bounds for q can also be found in [2, Theorem 2.3]. LEMMA 2.3. For 0 < q < 1, let the sequence (dj ) be given by d1 = 1 and

j ?1 ?q X kj dj ?k ; dj = 1 + q(1 ? q)j k=1

where kj is as in Lemma 2.1. Then, 1

j = 2; 3; : : :;

dj

sin q q q(1 ? q) j ;

j = 1; 2; 3; : : ::

ETNA

Kent State University etna@mcs.kent.edu

4

K. Diethelm

Remark 1. This lemma also holds in the limit cases q = 0 and q = 1, for then the recurrence relation reduces to dj = 1 and dj = 1 + dj ?1, respectively, which immediately implies dj = 1 or dj = j . Remark 2. A short calculation yields that 1 < (sin q)=( q(1 ? q)) 4= for every q 2 (0; 1). Proof. The inequality 1 dj is an easy consequence of the fact that kj > 0 for k 1 (cf. Lemma 2.1). We prove the upper bound for dj by induction. Since sin q q(1 ? q) 1 = d1 ;

the induction basis (j = 1) is presupposed. For the induction step, we de?ne a function Therefore, by Lemma 2.2 (ii), we have for every j

( ) = (1

x

? x)q .

Obviously,

00(x)

0.

j +1 X k=0

k;j +1

(

k=(j + 1)) = Qj +1

=

Z

]

0

1

( ) )(1) =

t t?q?1dt

Γ(?q)Γ(q + 1) = ? sin

Γ(?q)(Dq

q:

Using this result and the fact that kj

> 0 for k

1, we obtain

dj +1 = 1 + q(1 ? q)(j + 1)?q

1+ sin

j X k=1

k;j +1dj +1?k

j +1?k q k;j +1 j +1 k=1 ? sin q Qj +1 ] ? 0;j +1 (0) = 1+ ? sin q 0;j +1 (0) = qsin ?qq) (j + 1)q : (1

2.2. Proof of Theorem 1.1. We are now in a position to prove the main theorem. First of all, we note that we can actually evaluate the formula (1.5) because the denominator 0j ? (j=n)q Γ(?q) cannot vanish due to the fact that 0j < 0 (cf. Lemma 2.1) and our assumption 0. This implies that the denominator is strictly negative. Next, we recall that, for a ?xed j 2 f1; 2; : : :; ng,

j +1 qX

where (2.1)

t?q Z 1 t f (tj ) + x(tj ) = Γ(j q) x(tj ? ujqu)1 ? x(0) du + ? 0 t?q ? j = Qj j ] + Rj j ] ; Γ(?q) j (t) = x(tj ? tj t) ? x(0), and Rj is the quadrature error. Thus, f (tj ) +

j j X t?q X x(tj ) = Γ(j q) kj x(tj ?k) ? x(0) kj + Rj j ] : ? k=0 k=0

!

ETNA

Kent State University etna@mcs.kent.edu

Fractional differential equations

5

By construction,

j X k=0

kj = Qj 1] =

Z

0

1

u?q?1 du = ? 1 : q

Using this, we solve (2.1) for x(tj ) to obtain

x(tj ) =

0

j ? (j=n)q Γ(?q)

1

! j j q Γ(?q)f (t ) ? X x(t ) ? x(0) ? R ] j j j kj j ?k n q k=1

which, combined with (1.5),

j := x(tj ) ? xj

=

=

(

j=n)q Γ(?q)

1

? kj (x(tj ?k ) ? xj ?k ) ? Rj j ] j ? (j=n)q Γ(?q) k=1 ! j X ? 0j k=1 kj j ?k + Rj j ] :

1

0

j X

!

We now majorize this relation and ?nd, using Lemmas 2.1 and 2.2 (i), that

j jj

(

?

j=n)q Γ(?q) ? 0j k=1 kj j j ?kj + jRj j ]j !

1 1

0

j X

!

j X k=1

j

kj j j ?kj + q j q?2 j X k=1 kj j j ?kj +

j 1 q j q n?2 kx00k

!

00

q(1 ? q)j ?q

1 :

? Because of the initial condition, 0 = 0, and therefore j 1 j q(1 ? q) q nP2kx00k1 . Let ?q j ?1 kj dj ?k, us now de?ne a new sequence (dj ) by d1 = 1 and dj = 1 + q(1 ? q)j k=1 j = 2; 3; : : :. Then,

j j j q(1 ? q) q n?2 kx00k1 dj ; j jj

sin

j = 1; 2; : : :; n:

This is obvious for j = 1, and for j = 2; 3; : : :; n, it follows by a simple induction. Now, an application of Lemma 2.3 yields our ?nal result, viz.

q

q kx00k j q n?2 : 1

3. Numerical Examples. We have tried out the algorithm on some examples. In this ?nal section, we report the results. All the calculations were performed in standard doubleprecision arithmetic. For the ?rst example, we have chosen f (t) = t2 + 2t2?q =Γ(3 ? q), = ?1, and the initial condition x(0) = 0. The exact solution in this case is given by x(t) = t2 . For various choices of q 2 (0; 1), we have always obtained convergence orders close to O(nq?2 ). This is well in line with the prediction of Corollary 1.2. The resulting errors at t = 1 and the experimentally determined orders of convergence (“EOC") are shown in Table 3.1.

ETNA

Kent State University etna@mcs.kent.edu

6

Results for

=

K. Diethelm

?1 and f (t) = t2 + 2t2?q =Γ(3 ? q).

TABLE 3.1

n

5 10 20 40

q = 0:5 Error at t = 1 ?0:02087 ?0:00773 ?0:00282 ?0:00102

EOC

?1:43 ?1:46 ?1:47

q = 0:75 Error at t = 1 EOC ?0:05307 ?0:02312 ?1:20 ?0:00991 ?1:22 ?0:00421 ?1:24

q = 0:25 Error at t = 1 EOC ?0:00620 ?0:00199 ?1:64 ?0:00063 ?1:66 ?0:00020 ?1:68

The second example is f (t) = 2 cos t + t?q (1 F1 (1; 1 ? q; i t) + 1F1 (1; 1 ? q; ?i t) ? 2)=(2Γ(1 ? q)), = ?2, and the initial condition x(0) = 1. The exact solution in this case is given by x(t) = cos t. Again, we have always obtained convergence orders close to O(nq?2 ) for different values of q as expected according to Corollary 1.2. The results are shown in Table 3.2.

Results for

=

?2 and f (t) = 2 cos

t + t?q (1 F1 (1; 1 ? q; i t) + 1 F1 (1; 1 ? q; ?i t) ? 2)=(2Γ(1 ? q)).

TABLE 3.2

n

5 10 20 40

q = 0:5 Error at t = 1 ?0:03852 ?0:01572 ?0:00604 ?0:00225

EOC

?1:29 ?1:38 ?1:43

q = 0:75 Error at t = 1 EOC ?0:08303 ?0:03899 ?1:09 ?0:01746 ?1:16 ?0:00761 ?1:20

REFERENCES

q = 0:25 Error at t = 1 EOC ?0:01266 ?0:00448 ?1:50 ?0:00150 ?1:58 ?0:00048 ?1:63

[1] L. BLANK, Numerical treatment of differential equations of fractional order, Numerical Analysis Report 287, Manchester Centre for Computational Mathematics, 1996. [2] K. DIETHELM, Generalized compound quadrature formulae for ?nite-part integrals, IMA J. Numer. Anal., (to appear). [3] D. ELLIOTT, An asymptotic analysis of two algorithms for certain Hadamard ?nite-part integrals, IMA J. Numer. Anal., 13 (1993), pp. 445–462. [4] L. GAUL, P. KLEIN, AND S. KEMPFLE, Damping description involving fractional operators, Mech. Systems Signal Processing, 5 (1991), pp. 81–88. [5] R. GORENFLO AND R. RUTMAN, On ultraslow and intermediate processes, in Transform Methods and Special Functions, P. Rusev, I. Dimovski, and V. Kiryakova, eds., So?a, 1995, pp. 61–81. [6] K. B. OLDHAM AND J. SPANIER, The Fractional Calculus, vol. 111 of Mathematics in Science and Engineering, Academic Press, New York, London, 1974.