Matrix approach to discrete fractional calculus III: non-equidistant grids, variable step length and distributed orders

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.


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 fractionalorder 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].  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 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.   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 y = sin(t) 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.

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: and The exact analytical solution of this problem can be expressed using the Mittag-Leffler function:  Ay (t) + By (3/2) and 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.

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]: where w(α) denotes the weight function of distribution of order α ∈ [γ 1 , γ 2 ]. The weight function w(α) is normalized, so 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 and where matrices B α k are discrete analogues of fractional derivatives of orders α k on a given gridin 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: and 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): and u(t) = 0. (4.8) The discretization of equation (4.7) gives, as usual in the matrix approach, the following system of algebraic equations in matrix form: 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.

(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: 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].

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: 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]: 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) and that we already have y(t) = t in the interval (0,1), the problem (5.1)-(5.2) can be written as The integral in the last term can be easily evaluated as 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: and To solve this problem using the matrix approach [2,7], we need to obtain zero initial conditions. For this, we substitute y(t) = u(t) + 1, (5.8) and for the auxiliary function u(t) we have the desired initial value problem with zero initial condition: 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 C 1 D 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 0 P α a y(t) 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: H(t − a))y(t)). (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:

Linear fractional differential equations
Let us consider a linear fractional differential equation with constant coefficients in the interval If we assume that a k < n (k = 1, . . . , m) and n − 1 < max k a k < n, then we have to add n initial conditions, for example, The equation for the second 'large step' in the interval (a, b) will be The initial conditions should be, as usual, transformed to zero initial conditions. For this, we have to introduce the auxiliary function u(t): This process of making 'large steps' can be continued as long as necessary.

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

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 leftsided, 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.