Trapezoidal methods for fractional differential equations: theoretical and computational aspects

12/20/2019 ∙ by Roberto Garrappa, et al. ∙ University of Bari Aldo Moro 0

The paper describes different approaches to generalize the trapezoidal method to fractional differential equations. We analyze the main theoretical properties and we discuss computational aspects to implement efficient algorithms. Numerical experiments are provided to illustrate potential and limitations of the different methods under investigation.



There are no comments yet.


page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

It is nowadays well established that several real–life phenomena are better described by fractional differential equations (FDEs), where the term fractional, used for historical reasons, refers to derivative operators of any real positive order. Applications of FDEs are commonly found in bioengineering, chemistry, control theory, electronic circuit theory, mechanics, physics, seismology, signal processing and so on (e.g., [3, 7, 5, 33, 22, 41, 50]). We refer to [40] for an historical perspective on fractional calculus.

As a consequence of the growing interest for fractional order models, the analysis and the development of numerical methods for FDEs and partial differential equations with time or space fractional derivatives have become an active area of research (see, for instance,

[4, 18, 17, 13, 25, 26, 32, 35, 38, 44, 45, 47, 48, 46, 53, 54]).

Some of the methods for FDEs are directly derived from methods for integral equations; this is the case, for instance, of product integration (PI) rules. A completely different approach, named as fractional linear multistep methods (FLMMs) [38], has instead been specifically tailored to provide a solid theoretical background for the numerical treatment of FDEs.

Interestingly, these two approaches (PI and FLMM) lead to the same methods when the fractional order

converges to the nearest integer; indeed, different methods for FDEs are often generalizations of the same corresponding method for ordinary differential equations (ODEs).

Some important differences are however observed in the pure fractional case. For instance, the PI and FLMM generalizations of the implicit Euler method behave in a different way when applied to problems with discontinuous right–hand side [10, 11, 12] and only the method devised in the framework of FLMMs succeeds in stabilizing the solution on the switching surface, as expected from theoretical considerations, without generating unwanted spurious oscillations [24].

The presence of distinguishing characteristics in methods for FDEs developed from the same method for ODEs deserves to be more thoroughly investigated. Some differences are of theoretical interest (convergence, stability, etc.); nevertheless, the two approaches also present distinctive features of computational nature that affect the efficiency of the solution process.

The main aim of this paper is to describe and compare some methods for FDEs, highlight strengths and weaknesses and investigate the circumstances under which a certain method is preferable to some others.

The analysis will be confined to those PI rules and FLMMs that generalize, in some respects, the implicit second–order trapezoidal rule. This –stable one–step method has the smallest possible error constant [34] and is appreciated since provides a satisfactory balance between accuracy, stability properties and computational complexity. Our interest is in exploring when these strengths are inherited by the corresponding methods for FDEs.

There are also some more specific motivations for restricting the attention to second–order methods: a) PI rules very seldom converge with an order equal or greater than 2; b) the implementation of higher order FLMMs is usually computationally demanding; c) as for multistep methods for ODEs, order 2 represents a barrier for stability.

Despite this apparently limited extent, we will however be able to include in our analysis a sufficiently large number of methods.

The choice of focusing on implicit methods is also motivated by the fact that they have been so far very seldom investigated (very often, implicit methods are taken into account in the framework of predictor–corrector algorithms in which the implicit scheme is actually implemented as an explicit one) and explicit methods are unduly preferred in most of the works on FDEs, regardless of stability issues which discourage their use.

The paper is organized as follows. In Section 2 we present some basic definitions and properties concerning FDEs. Section 3 is devoted to describe PI rules and discuss different implementations on uniform and graded grids. In Section 4 we present FLMMs and we describe three different generalizations of the implicit trapezoidal rule; the use of an effective formula for the evaluation of the weights is also proposed. In Section 5 linear stability is investigated and some implementation details are discussed in Section 6. Finally, in Section 7 the results of some numerical tests are presented to compare the methods under investigation and some concluding remarks are presented at the end of the paper.

2 Fractional differential equations

Let us consider the initial value problem for a system of FDEs in the form


where is the smallest integer such that , is assumed to be sufficiently smooth, is the unknown solution and denotes classical derivative of integer order .

The fractional derivative operator is introduced in this paper according to the Caputo’s definition [14, 42, 49]

with the Riemann–Liouville (RL) integral of order on the interval


A classical result in fractional calculus [14] allows to compute, for any , the RL integral of the power function as


It is a well known result that (1) can be equivalently written in terms of the integral equation


where is the Taylor expansion of centered at

One of the major difficulties in the treatment of FDEs derives from the lack of smoothness of the true solution on the whole interval of integration. Indeed, as proved in [36], the true solution of (1) possesses an expansion in mixed (integer and real) powers of the form


with , and

with the consequence that the derivatives of are unbounded at .

3 Product integration rules

The main idea behind PI rules [52] is to evaluate the integral in (4

) by approximating the vector field

with suitable polynomials. For the integer–order case this is the way in which Adams multistep methods are devised.

To generalize the –step trapezoidal rule consider a (not necessarily uniform) grid on the whole interval of integration and, in the piecewise decomposition of (4)


in each subinterval by the first–degree polynomial interpolant

where , and is the numerical approximation to the solution . After exactly solving the integrals

where, for shortness, we denoted

it is elementary, after a simple change of the summation index, to obtain the numerical approximation of the solution at as



According to the analysis presented in [6], this scheme is convergent of second order whenever the true solution is assumed sufficiently smooth. Unfortunately, this is a rather optimistic result compared to the expansion (5) of : due to the presence of unbounded derivatives at the convergence rate of (6) is actually lower than the theoretical value , as it can be formally stated by a direct application of the main result by Dixon in [19].

Theorem 1.

Let be Lipschitz continuous with respect to the second variable and the numerical approximation obtained by applying the PI trapezoidal rule (6) on the interval . There exists a constant , which does not depend on , such that

Unless very restrictive (and usually unrealistic) assumptions are made to ensure the smoothness of also at the left endpoint of , the rate by which the error reduces to is therefore, away from , proportional to when . The same drop in the convergence order occurs also when higher degree polynomials are employed; for this reason PI rules of higher order are never taken into account for FDEs when . When order of convergence is instead obtained.

3.1 Uniform meshes

Using (6) on uniform grids , with a constant step–size, offers some advantages of computational type. The evaluation of the weights is indeed less expensive and just involves the computation of real powers of integer numbers. Furthermore, (6) presents a convolution structure, in which each weight depends just on , and it can be simplified as


where now

Formulation (7) of the trapezoidal PI rule is well–known since it is a building component of a widely used predictor–corrector algorithm [18].

The main advantage of the convolution structure is however that (7) can be implemented in a fast way by means of FFT algorithms [28].

3.2 Graded meshes

Unlike FLMMs, PI rules can be implemented on non uniform grids, with the aim of overcoming the order reduction stated in Theorem 1. This issue has been investigated, in the context of weakly singular integral equations, by several authors [1, 9] and graded grids are highly praised to this purpose.

By clustering the grid–points at the left boundary of the integration interval, graded grids are able to cope with the lack of regularity and counteract the inefficiency of polynomials to satisfactorily approximate real powers.

A graded mesh on the interval is defined according to

where is a suitable grading exponent. After denoting and observing that , a simple manipulation allows to reformulate (6) according to



The following result [1, 2] on the second order convergence of graded grids for linear FDEs can be readily extended to more general equations.

Theorem 2.

Let , . For it is

The use of PI rules on nonuniform grids is more demanding; not only the weights must be recomputed at each step but, since the absence of a convolution structure, it is not possible to employ fast algorithms. By means of some numerical experiments we will verify, at the end of the paper, whether or not this further amount of computation is, in some way, balanced by the enhancement in the accuracy.

4 Fractional linear multistep methods

In the pioneering work on discretized fractional calculus [38], Lubich proposed an elegant and effective strategy for classical linear multistep methods (LMMs) originally devised for integer–order ODEs.

Although the potentials of FLMMs have been recognized in several papers (see, for instance, the review article by Gorenflo [27]), their use for practical computation has not captured so far widespread attention, perhaps due to some difficulties in explicitly formulating FLMMs.

The key aspect in FLMMs is the approximation of the RL integral (2) by means of the convolution quadrature


on uniform grids , , and where convolution and starting quadrature weights and do not depend on .

Starting weights play an important role, especially in the first part of the integration interval, in order to cope with the possibly singular character of the integrand function at ; we will discuss the problem of their evaluation later on.

Convolution quadrature weights are the main components of the quadrature rule and characterize the specific FLMM. They are obtained starting from any LMM for ODEs

being and the characteristic polynomials. As shown in [51], LMMs can be equivalently reformulated as (8) (but with for ODEs) with weights obtained as the coefficients of the formal power series (FPS)

and the function goes therefore under the name of generating function of the LMM. The idea proposed in [38] is to generate a quadrature rule for fractional problems (1) by evaluating the convolution weights as the coefficients in the FPS of the fractional–order power of the generating function


Methods of this kind, named as FLMMs, when applied to (1) read as


and their convergence properties are stated in the following result [37, 38].

Theorem 3.

Let be a stable and consistent of order implicit LMM with the zeros of having absolute value . The FLMM (10) is convergent of order .

One of the major difficulties in FLMMs (10) is in evaluating the weights as the coefficients in the FPS (9). Although some sophisticated algorithms are available for manipulating FPS [31], for most of the methods an efficient tool is the J.C.P. Miller formula stated by the following theorem [30, Theorem 1.6c].

Theorem 4.

Let be a FPS. Then for any ,

where coefficients can be recursively evaluated as

The Miller formula allows to evaluate, in the most general case, the first coefficients of with a number of operations proportional to . In most cases of practical interest, this computational effort is actually reduced by an order of magnitude. This is the case of the evaluation of the coefficients in the FPS of ; indeed, since and , it is elementary to see that the application of Theorem 4 leads to


involving just multiplications and additions.

4.1 Fractional trapezoidal rule

In order to generalize to FDEs the classical trapezoidal rule, let us first consider its classical formulation for ODEs


with characteristic polynomials and and generating function

To practically compute the first weights in the FPS of we advise against the direct application of the Miller formula to the infinite series in since it would involve a number of operations proportional to . It is instead preferable to first apply twice Theorem 4 to and , by means of (11), and hence evaluate the coefficient of their product by an FFT algorithm, with a number of operations proportional to when is a power of [31]. This last task can be efficiently performed by means of the very few lines of Matlab code

x = fft([omega1,zeros(size(omega1))]) ;
y = fft([omega2,zeros(size(omega2))]) ;
omega = ifft(x.*y) ; omega = omega(1:N) ;

where omega1 and omega2 are the row arrays with coefficients of and respectively, evaluated in a fast way thanks to (11).

4.2 Newton–Gregory formula

The trapezoidal rule (12) belongs to the more general family of –step Adams–Moulton methods

with the classical backward differences (i.e., and ); the coefficients are the first terms in the truncation

of the power series expansion, with respect to , of and the resulting generating function is .

It is easy to see that and , thus immediately leading to (12) when . Therefore, a straightforward generalization to fractional–order problems by means of the –power of leads, in the case of the –step method, to the trapezoidal rule discussed in the previous subsection.

It has been shown in [36] that alternative methods, named as Newton–Gregory formulas, can be devised by first evaluating the power of and successively truncating to the first terms as

The interchange between the truncation series symbol and the –power is important and generates a different family of methods. By operating for simplicity the change of variable ,

The first coefficients can be easily evaluated thanks again to Theorem 4 (see [21] for the explicit formulation of the first coefficients).

Since and , the generating function of the FLLM obtained by the –step Adams method is therefore

and, after a simple manipulation, it is possible to evaluate the first coefficients in , with a cost proportional to , as

with the coefficients in the FPS evaluated by means of (11).

4.3 Fractional BDF formula

The second order BDF formula for ODEs is given by


with characteristic polynomials and and generating function

BDF method (13) can be considered in no way as coming from the trapezoidal rule; it is indeed obtained by differentiating a second–degree interpolant polynomial. However, it is equally of interest in our analysis due to its second order of convergence as the trapezoidal rule and the excellent stability properties that we are interested in verifying also for FDEs.

The coefficients in the FPS of can be profitably evaluated again by the Miller formula applied to , with negative power , to first produce

and hence by computing with an overall number of operations that still remains proportional to .

5 Linear stability analysis

The analysis of stability is of primary importance to understand the qualitative behavior of a numerical method. By assuming the existence of a stable steady state, stability analysis aims to verify whether or not approximated trajectories move towards the steady state.

As usual, stability issues for FDEs can be addressed by studying the scalar linear test equation


whose steady–state is at . The generalization of a well–known result in the theory of ODEs can be done by means of the following result [37, 43].

Theorem 5.

Let . The steady–state of (14) is stable if and only if , where .

In Figure 1 we present the sector of stability for two exemplifying cases. It is immediate to verify that coincides with the left complex semi–plane when , thus generalizing the classical result for ODEs.

Figure 1: Stability sectors for (left) and (right)

When PI rules or FLMMs are applied to (14), the numerical solution is evaluated by means of a convolution quadrature


The problem of determining the stability region such that whenever has been discussed in [37, 39] and the main result can be summarized in the following theorem.

Theorem 6.

Let . Under the assumption that the sequence is convergent and that the quadrature weights satisfy

the stability region of the convolution quadrature (15) is

We refer to [23] and [38] to verify that the PI rule and the FLMMs investigated in this paper fulfill the assumptions of Theorem 6.

In the numerical treatment of ODEs –stability means that the whole complex semi–plane is included in the stability region of the method. –stable methods are therefore particularly attractive since they ensure that the numerical solution of stable linear problems converges to regardless of the size of the step . Because of the result in Theorem 5, for FDEs it is more useful to use the concept of –stability, which means that the whole sector must be included in the stability region.

While is a priori known for FLMMs, the generating function of PI rules does not possess an analytical representation and, therefore, it can be just numerically approximated.

In Figure 2 we present a comparison of the stability regions (the gray areas) of the methods discussed in the previous sections for various values (dotted lines denote the corresponding sectors of stability ). To identify each method, here and throughout the following section we make use of the acronyms listed in Table 1.

Acronym Formula Subsection
PI U PI rule on uniform grid 3.1
PI G PI rule on graded grid 3.2
FT fractional trapezoidal rule 4.1
NG Newton–Gregory formula 4.2
FBDF fractional BDF formula 4.3
Table 1: Acronyms used to denote methods

We can distinctly observe that, as increases from to , the stability region of each method decreases but always includes ; we can draw the conclusion that all the methods are –stable for .

Figure 2: Stability regions (gray areas) for some orders

Within ODEs, the stability region of the trapezoidal rule is the negative semi–plane and exactly corresponds to the region in which the true solution converges to the steady–state. Among the methods presented in this paper, just FT preserves this feature and indeed its stability regions coincide with for any value of . NG and PI have regions of almost the same size but larger with respect to FT; as in the ODE case, FBDF possesses the largest regions of stability.

Interestingly, the situation is completely different for as shown in Figure 3. Whilst the stability regions of FT still coincide with the stability sectors and the stability properties of FBDF still remain the most remarkable ones, we observe that NG and PI completely lose the –stability. This is an important aspect to be taken into account: despite the implicit nature of these methods, numerical instability can arise when a not sufficiently small step–size is used and, therefore, they seem not suited to solve stiff problems when .

Figure 3: Stability regions (gray areas) for some orders

6 Implementation details and Matlab codes

In this section we discuss some implementation issues that are essential for the development of efficient codes.

6.1 Evaluation of the lag term

The evaluation of the (usually long) lag term is one of the most expensive tasks in convolution quadratures like (7) or (10). The integration along a large number of steps is usually a critical bottleneck and, if not well designed, it can be an impracticable task involving a number of operations proportional to .

Several authors have discussed this problem (see, for instance, [8]); in [15] a strategy based on parallel algorithms has been described.

In our codes we implement the FFT algorithm described in [28] by which the computational effort downgrades from to . We do not further discuss technical details of this algorithm, for which we just refer to the above mentioned paper; we simply observe that, without efficient algorithms of this kind, the numerical solution of FDEs very quickly degenerates in an unworkable task as the number of grid–points increases.

6.2 Starting weights

Starting weights are introduced in (10) to cope with the singular behavior of the solution close to the left endpoint of the integration interval. They are evaluated after imposing that (8) is exact for with and the order of convergence of the FLMM [38].

Because of (3), and since the methods investigated in this paper have order , it is immediate to see that weights must satisfy the relation


and, hence, they are evaluated by solving, at each step , a system of linear equations, with the cardinality of .

This problem has been in–depth discussed in [16], where the authors stressed the attention on the mildly ill–conditioned nature of the coefficient matrix of (16). Nevertheless, for the methods under investigation this is not a serious issue since the size of is small, being (in most cases will be 3 or 4 and the inversion of can be done analytically). The only exception is for very small values of which, however, are extremely uncommon in applications.

6.3 Starting the computation

The initialization of (10) requires the knowledge of the first approximations of the solution. Since just is given by the problem, the remaining values must be evaluated in some other way, for instance by using some different methods.

To avoid mixing methods of different nature, we prefer to evaluate all the approximations at the same time by setting the single implicit system


This nonlinear system has size , where is the size of the system of FDEs we want to solve. No particular problems arise when the size of the problem is small; otherwise some memory allocation problems could occur in presence of very small values of , involving a very large number of initial approximations to compute.

6.4 Solution of the nonlinear systems

All the implicit methods analyzed in Sections 3 and 4 require, at each step, the solution of an algebraic nonlinear system in the form

where and do not depend on . To preserve the good stability properties, it is convenient to use Newton iterations

where is solution of the linear system

and, usually, . We refer to [29] for more insights on this issue and for a discussion on the stopping criteria.

We just highlight here that the convergence can depend on , instead of as in ODEs. Whilst this is not an issue when , some unexpected problems could arise whenever is small. For this reason it is necessary to use a quite small step–size with low values of .

7 Numerical experiments

By means of numerical experiments we compare the methods discussed in the paper. We perform all the experiments in Matlab, version, on a computer equipped with the Intel Dual Core E5400 processor running at 2.70 GHz under the Windows XP operating system; all the Matlab codes have been optimized, to get the best performance, according to the discussion reported in the previous sections.

As a first test problem we consider the linear FDE (14) for . For each method, and for an increasing number of grid–points, we report the error at

, with respect to a reference solution, and we provide an estimation of the order of convergence (EOC) as


Table 2: Errors and EOC at for the linear FDE (14) with and

The numerical results in Table 2 for clearly confirm theoretical findings on the order of convergence. As expected from Theorem 1, the PI rule on uniform grids fails to obtain order 2 and converges with order ; on the other hand, the PI rule on graded grids converges with full order .

As in the ODE case, FT produces the smallest error. Although also NG and PI are generalizations of the trapezoidal rule, it seems that only FT inherits this feature; by taking into account, from the analysis of Section 5, that FT preserves the other peculiar character of the trapezoidal rule of having the region of stability exactly equal to the region in which the solution converges towards the steady state, it is tempting to conclude that FT is the most natural generalization of the trapezoidal rule.

Unfortunately, the situation is more complicate since when it is NG that provides the lowest error, as shown in Table 3 for .

Table 3: Errors and EOC at for the linear FDE (14) with and

This phenomenon can be more clearly observed in Figure 4, where we show the errors as a function of (we used the same values for and from the previous experiment and a fixed number of grid points).

Figure 4: Errors at for the linear FDE (14) with and .

We must highlight that this is a not completely unexpected result. Stability and accuracy are indeed conflicting requirements in most of numerical methods; thus, FT shows the smallest regions of stability but the lowest error when , exactly as it happens with NG for . What we observe is a kind of roles swap between FT and NG when moves from the left to the right of 1.

The very bad performance of PI on graded grids for low values of is due to the fact that the grids generated by the graded exponent are very coarse close to right endpoint of the integration interval. Furthermore, graded grids do not add any improvement for , since the maximum expected order is already achieved by using uniform grids. We can draw the conclusion that the utility of graded grids is confined only to .

It is also interesting to analyze the computational effort in relation to the accuracy. To this purpose, the plots in Figure 5, which are related to the experiments presented in Tables 2 and 3, show that FT, NG and FBDF behave in a more efficient way with respect to PI when ; indeed, PI U pays for the reduction in the accuracy whilst PI G pays for the absence of a convolution structure which precludes the application of FFT algorithms. The minor efficiency of FBDF with respect to FT and NG is due to the major complexity in the evaluation of weights. It is worthwhile to note that the preeminence of FLMMs with respect to PI rules holds despite the need of solving linear systems for evaluating the starting weights. For the absence of order reduction increases the efficiency of PI on uniform nodes but however it is below that of NG which instead has lower errors.

Figure 5: Errors versus time of execution for (left) and (right)

To observe the behavior of the methods under investigation for nonlinear problems, we consider the fractional reaction diffusion system of Brusselator type


It has been shown [20] that there exists a marginal value