## Abstract

In this paper, we further develop Podlubny’s matrix approach to discretization of integrals and derivatives of non-integer order. Numerical integration and differentiation on non-equidistant grids is introduced and illustrated by several examples of numerical solution of differential equations with fractional derivatives of constant orders and with distributed-order derivatives. In this paper, for the first time, we present a variable-step-length approach that we call ‘the method of large steps’, because it is applied in combination with the matrix approach for each ‘large step’. This new method is also illustrated by an easy-to-follow example. The presented approach allows fractional-order and distributed-order differentiation and integration of non-uniformly sampled signals, and opens the way to development of variable- and adaptive-step-length techniques for fractional- and distributed-order differential equations.

## 1. Introduction

The results that we present in this paper were motivated by two important challenges in applied numerical methods of fractional calculus.

First, until recent times, the fractional derivatives were discretized using equidistant nodes. This has roots in the famous Grünwald–Letnikov definition of fractional-order differentiation, which is based on generalization of finite differences defined on an equidistant grid, and which gives the simplest and very efficient approximation for numerical evaluation of fractional derivatives. This Grünwald–Letnikov-based approach to discretization of fractional derivatives had very strong impact on the way of thinking in fractional calculus, that even fractional integrals were routinely discretized on equidistant grids, too. However, it is clear that for fractional integrals it was not a necessity at all. On the other hand, it was unclear what one could do with approximation of fractional integrals on non-equidistant grids, if one wants to have a uniform approach to discretization of both fractional derivatives and fractional integrals.

Second, solution of fractional differential equations in large time intervals is a computational problem due to the fact that the number of points taken into account in computations grows with the growing value of the time variable. This is caused by the non-local character of fractional-order differentiation. So far, the only aid in this respect was the ‘short memory principle’ [1]. Methods known from classical numerical solutions of integer-order differential equations, especially variable-step-length techniques, were not available for fractional differential equations.

The systematic and continuous development of the ‘matrix approach’ [2] allowed us to find some answers to both challenges, and in this paper, we present them as two mutually related methods for solving problems of discrete fractional calculus on non-uniformly spaced discretization grids. Moreover, we extend this approach to distributed-order operators and distributed-order differential equations.

We start with demonstrating how the matrix approach can be extended to numerical evaluation of fractional-order integrals and derivatives on non-equidistant grids, and how fractional differential equations with constant-order derivatives can be solved on such grids. This finally unifies the discretization of fractional derivatives and fractional integrals on arbitrary (equidistant and non-equidistant) grids.

After that we make the next step and extend the matrix approach to discretization of distributed-order operators and to numerical solution of distributed-order differential equations.

Then, we move the focus onto using variable step length for solving fractional differential equations. In this paper, we, for the first time, present the method that we call ‘the method of large steps’. We provide the general framework and illustrate this method by a numerical example that, for verification purposes, allows easy exact solution as well.

Because each ‘large step’ consists of a set of ‘small steps’, it can be done using the matrix approach, and the ‘small steps’ can be equidistant or non-equidistant. This is illustrated by included little pieces of Matlab code using our publicly available toolbox [3,4].

The methods presented in this paper finally allow fractional-order differentiation and integration of non-uniformly sampled signals, and the development of variable-step-length techniques for solving fractional differential equations (ordinary and partial).

The standard basic notation and basic definitions of fractional derivatives and fractional integrals can be found in Podlubny [1], Samko *et al.* [5] and Kilbas *et al.* [6].

## 2. Fractional-order integration and differentiation on non-equidistant grids

Although equidistant grids are frequently used in application, in many situations, the use of non-equidistant grids brings notable advantages. For example, many numerical methods for solving differential equations use a variable-time-step technique, so the time step can increase or decrease depending on how rapidly the resulting solution is changing.

Up to now, the fractional derivatives were discretized using equidistant nodes. This was, of course, due to the famous Grünwald–Letnikov definition of fractional-order differentiation, which is based on generalization of finite differences defined on an equidistant grid.

The matrix approach to discretization of integrals and derivatives of arbitrary real order, developed by Podlubny [2,7], allows us to generalize discretization of fractional-order integrals and derivatives to non-equidistant grids.

The idea is to create first a discretization matrix *I*^{α} for integration of order *α*. After the matrix *I*^{α} for discrete fractional integration on a non-equidistant grid is obtained, we can easily derive the matrix *D*^{α} for discretization of fractional-order derivatives by matrix inversion:

In the simplest case, the function under differentiation can be approximated by a piecewise constant function, and for the non-equidistant discretization nodes *t*_{k} (*k*=1,…,*N*), the coefficients of the lower triangular matrix *I*^{α} can be evaluated as
2.1
with
2.2
Other formulae for numerical integration will give other expressions for the coefficients *I*_{k,j}; however, as we demonstrate below, even this simple approximation works well.

In the case of non-equidistant nodes, however, the matrices *I*^{α} and *D*^{α} are not strip matrices, although for one-sided fractional integrals and derivatives they are still triangular.

In the examples given below, we use non-equidistant nodes generated with random steps. We generate *N* random points between 0 and 1, sort them in ascending order, and then scale to the considered interval of length *L*. After that, we replace the first and the last randomly generated nodes with the exact left and right bounds of the considered interval.

### (a) Example 1: evaluation of integer-order integrals and derivatives

Let us consider the function . In figure 1, its exact first-order derivative (dashed line) and its exact (dotted line) onefold integral (dash-dotted line) are depicted.

The results of numerical differentiation (solid line) and numerical integration (dash-dotted line) of the same function using the matrix approach on a non-equidistant grid of 200 randomly generated nodes in the interval [0,2] are also shown in figure 1.

### (b) Example 2: evaluation of fractional-order integrals and derivatives

The proposed approach is suitable also for evaluating fractional-order integrals and derivatives. In figures 2 and 3, fractional-order integrals and derivatives of orders *α*=0.1, *α*=0.3, *α*=0.5 and *α*=0.7 of function are plotted. Each derivative obtained using non-equidistant steps is compared with the solution obtained using the ‘matrix approach’ with equidistant steps. The match shows good agreement of the results.

## 3. Solution of fractional differential equations on non-equidistant grids

Using discrete analogues of fractional-order derivatives on a non-equidistant grid, we can easily and conveniently perform discretization and numerical solution of fractional differential equations on such grids. We illustrate the developed approach on two classical benchmark problems: the two-term fractional relaxation equation and the Bagley–Torvik equation.

### (a) Example 3: fractional relaxation equation

In the first work on the matrix approach to discrete fractional calculus [2], the following sample two-term fractional differential equation in terms of the Caputo derivatives [1] under zero initial conditions was considered: 3.1 The exact analytical solution of this problem can be expressed using the Mittag–Leffler function: 3.2

In figure 4, the comparison of the exact analytical solution (dotted line) and numerical solution obtained with the help of developed approach (solid line) using non-equidistant nodes (with *N*=500) for the case of *α*=1.8 is shown.

### (b) Example 4: fractional oscillator equation

As mentioned earlier, the proposed approach allows easy solution of ordinary differential equations with derivatives of arbitrary real order (integer and non-integer). Let us consider the following classical initial value problem for the Bagley–Torvik equation (also known as damped fractional oscillator equations): 3.3

The solutions of the Bagley–Torvik equation for *A*=1, *B*=1, *C*=1 in the time interval [0;30], obtained using two approaches, are shown in figure 5, and the solutions of the same problem for *A*=1, *B*=0.5, *C*=0.5 are depicted in figure 6. In both cases, the number of discretizaton steps is 400. The dotted lines represent the numerical solutions obtained using equidistant steps (with *h*=0.075), and the solid lines represent the numerical solutions with the same number of randomly generated non-equidistant steps.

## 4. Solution of distributed-order differential equations on non-equidistant grids

The presented extension of the matrix approach to discretization of non-integer order integrals and derivatives and to numerical solution of equations with such operators on non-equidistant grids can be used for solving distributed-order differential equations, too. In this paper, we use the following definition of distributed-order differential/integral operators [8]:
4.1
where *w*(*α*) denotes the weight function of distribution of order *α*∈[*γ*_{1},*γ*_{2}]. The weight function *w*(*α*) is normalized, so
4.2

The idea of distributed-order differential equations was first introduced most probably by Caputo [9,10].

As Jiao *et al.* [8] showed recently, distributed-order derivatives can be discretized as
4.3
and
4.4
where matrices *B*^{αk} are discrete analogues of fractional derivatives of orders *α*_{k} on a given grid—in our case, on a non-equidistant grid. The matrices *B*^{αk} for discrete differentiation of order *α*_{k} on a non-equidistant grid are obtained as described earlier, and then the matrix *B*^{w(α)} for discrete distributed-order differentiation on the same non-equidistant grid is computed using relationship (4.4).

### (a) Example 5: distributed-order relaxation equation

Let us consider the following initial value problem for the distributed-order relaxation equation:
4.5
and
4.6
where the distribution of the orders *α* is given by the function *w*(*α*)=6*α*(1−*α*) (0≤*α*≤1). To be able to use the matrix approach, we need a zero initial condition. Introducing an auxiliary function *u*(*t*),
gives the following initial value problem for the new unknown *u*(*t*):
4.7
and
4.8

The discretization of equation (4.7) gives, as usual in the matrix approach, the following system of algebraic equations in matrix form:
4.9
where *U*_{n} is the vector of the values of *u*(*t*) at the discretization nodes, and *F*_{n} is the vector of the values of the right-hand side, *f*(*t*)−*b*, at the same nodes; *E*_{n} is the identity matrix.

The results of computations for the case of 500 randomly generated nodes in the interval [0,5] for *b*=0.1 are shown in figure 7 by the solid line. They are in agreement with the numerical solution obtained for the same number of uniformly spaced nodes (dotted line in figure 7).

### (b) Example 6: distributed-order oscillator

Let us consider an initial value problem for the Bagley–Torvik equation, where the damping term is expressed in terms of distributed-order derivatives: 4.10 and 4.11

Similar to example 5, we just replace continuous operators with their corresponding matrix-form discrete analogues on the considered non-equidistant grids, and the known and unknown function by the vectors of their values in the discretization nodes. This gives the following algebraic system in matrix form: 4.12

The results of computations for the case of 400 randomly generated nodes in the interval [0,30] for *A*=1, *B*=1, *C*=1, *φ*(*α*)=6*α*(1−*α*), *α*∈[0,1], are shown in figure 8 by the solid line. They are in obvious agreement with the numerical solution obtained for the same number of uniformly spaced nodes (dotted line in figure 8).

The Matlab implementation of the matrix approach to discretization of distributed-order operators and numerical solution of distributed-order differential equations can be found in the update to the available toolbox [11].

## 5. Method of ‘large steps’

We will illustrate the idea of the proposed ‘method of large steps’ on an easy-to-follow and sample problem, which allows exact analytical solution. In this section, we prefer using small snippets of the Matlab code in order to illustrate the simplicity of the procedure, and the idea of how ‘large steps’ is performed.

### (a) Sample problem in the interval (0, 2)

It is worth remembering that we use the Caputo derivatives [1]. Let us consider the following sample problem for the large-steps method:
5.1
and
5.2
One can easily verify that the exact solution of this problem is *y*(*t*)=*t*.

### (b) First ‘large step’: interval (0,1)

We can solve the problem (5.1)–(5.2) numerically in the interval (0,1) using the recently developed matrix approach [2,7] and the corresponding MATLAB toolbox [2]:

`clear all`

`h = 0.01;`

`t = 0:h:1;`

`N = 1/h + 1;`

`M = zeros(N,N);`

`M = ban(0.5, N, h) + eye(N,N);`

`F = (t.^(0.5)/gamma(1.5) + t)’;`

`M = eliminator(N,[1])*M*eliminator(N,[1])’;`

`F = eliminator(N,[1])*F;`

`Y = M\F;`

`Y0 = [0; Y];`

`plot (t,Y0,’b’)`

`set(gca, ’xlim’, [0 2], ’ylim’, [0 2] )`

`grid on, hold on`

So, we have solved (see figure 9) the previous problem in (0,1) and we know *y*(*t*) for *t* in (0,1). How can we continue from the point at which we have arrived?

### (c) Second ‘large step’: interval (1,2)

Taking into account that for *t*>1 (we recall that we use the Caputo derivatives here)
5.3
and that we already have *y*(*t*)=*t* in the interval (0,1), the problem (5.1)–(5.2) can be written as
5.4

The integral in the last term can be easily evaluated as 5.5

Now, we are ready to make the second ‘large step,’ i.e. solution in the interval (1,2). In the interval (1,2) (second step), we have to solve the following problem: 5.6 and 5.7

To solve this problem using the matrix approach [2,7], we need to obtain zero initial conditions. For this, we substitute
5.8
and for the auxiliary function *u*(*t*) we have the desired initial value problem with zero initial condition:
5.9
and
5.10

Now we solve the problem for *u*(*t*) using the same matrix approach toolbox, and plot the solution:

`clear all`

`h = 0.01;`

`t = 1:h:2;`

`N = 1/h + 1;`

`M = zeros(N,N);`

`M = ban(0.5, N, h) + eye(N,N);`

`F = (t.^(0.5)/gamma(1.5) + t - 2*t.^(0.5)/gamma(0.5) ... + 2*(t-1).^(0.5)/gamma(0.5) - 1)’;`

`M = eliminator(N,[1])*M*eliminator(N,[1])’;`

`F = eliminator(N,[1])*F;`

`U = M \F;`

`U0 = [0; U];`

`Y0 = U0 + 1;`

`plot(t, Y0, ’g’)`

We see (figure 10) that finally we obtain the solution of the original problem in the interval (0,2) using two ‘large steps’: the first step was numerical solution in (0,1), and the second step was numerical solution in (1,2). In the right-hand side of the equation for the interval (1,2) two additional terms appeared as the result of considering fractional differentiation with a different lower terminal, namely .

## 6. Method of ‘large steps’: general scheme

In general, if we have considered the problem (0<*α*<1)
6.1
and
6.2
and obtained its solution in the interval (0,*a*) (and the final value *y*_{a} at *t*=*a*), then we can use this for transforming the problem to
6.3
and
6.4
where
6.5
is the contribution of the ‘past’ of the process *y*(*t*) in the interval [0,*a*] to the differential equations describing its current state in the interval [*a*,*b*].

It is useful to note here that can be evaluated as a fractional derivative of the function *y**(*t*)=(1−*H*(*t*−*a*))*y*(*t*), where *H*(*t*) is the Heaviside unit step function:
6.6

Introducing an auxiliary function *y*(*t*)=*u*(*t*)+*y*_{a}, we arrive at the problem with zero initial condition for the function *u*(*t*), which can be solved numerically:
6.7
and
6.8

This process of making ‘large steps’ can be continued as long as necessary.

## 7. Linear fractional differential equations

Let us consider a linear fractional differential equation with constant coefficients in the interval (0,*b*),
7.1
If we assume that *a*_{k}<*n* (*k*=1,…,*m*) and , then we have to add *n* initial conditions, for example,
7.2
The equation for the second ‘large step’ in the interval (*a*,*b*) will be
7.3
and the initial conditions for the second ‘large step’ will have the values of the final values of the solution in the first interval:
7.4

The initial conditions should be, as usual, transformed to zero initial conditions. For this, we have to introduce the auxiliary function *u*(*t*):
7.5

This process of making ‘large steps’ can be continued as long as necessary.

## 8. Method of ‘large steps’ and the problem of initialization of fractional derivatives

Lorenzo & Hartley [12,13] raised the question about initialization of fractional derivatives. Their motivation was to use or recover the information about the process *y*(*t*) in the interval (0,*a*), if we consider fractional derivatives of *y*(*t*) in (*a*,*b*).

It is worth noting that in the second ‘large step’ in the considered sample problem we used, in fact, the proper initialization of the fractional derivative in the interval (1,2) based on the known behaviour of *y*(*t*) in (0,1).

In other words, proper initialization in the interval (*a*,*b*) can be done only when we know all values of *y*(*t*) in the preceding interval (0,*a*).

## 9. Concluding remarks

The presented extension of the matrix approach to discretization of non-integer order derivatives and integrals of constant and distributed order allows numerical solution of differential equations with such derivatives on non-equidistant grids.

The matrix approach proves to be a very easy, algorithmic, modular and convenient method that unifies numerical solution of many types of problems in various settings. In the examples in this paper we, for simplicity, used only ordinary differential equations with generalized derivatives; partial differential equations on non-equidistant grids can be solved using the technique published earlier [7].

Of course, the proposed method can be used for equations containing any mixture of left-sided, right-sided and two-sided derivatives of integer orders, constant non-integer orders, and distributed orders, on equidistant (uniform) or non-equidistant (non-uniform) grids. In all cases, those derivatives are simply replaced with their discrete analogues in the form of easily computable matrices.

The methods presented in this paper finally allow fractional-order differentiation and integration, and also distributed-order differentiation and integration, of non-uniformly sampled signals.

The proposed method of ‘large steps’ allows one to avoid in many situations the limitations of the ‘short memory principle’ and to obtain numerical solutions in large intervals with higher accuracy. At the same time, the method of ‘large steps’ finally allows the development of variable- and adaptive-step-length techniques for solving differential equations of non-integer order (ordinary and partial).

We would like to mention that other kinds of efforts towards using non-equidistant grids for numerical solution of fractional differential equations can be found in recent studies [14,15], and that some existing methods, such as the collocation method [16,17] or explicit numerical methods [18], can be reconsidered in terms of non-equidistant grids as well.

## Acknowledgements

This work was supported in part by grants (nos APVV 0040-07, SK-UA-0042-09, and VEGA 1/0497/11, 1/0729/12 and 1/0746/11).

## Footnotes

One contribution of 14 to a Theme Issue ‘Fractional calculus and its applications’.

© 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/3.0/, which permits unrestricted use, provided the original author and source are credited.