Monte Carlo method for fractional-order differentiation

In this work the Monte Carlo method is introduced for numerical evaluation of fractional-order derivatives. A general framework for using this method is presented and illustrated by several examples. The proposed method can be used for numerical evaluation of the Grünwald-Letnikov fractional derivatives, the Riemann-Liouville fractional derivatives, and also of the Caputo fractional derivatives, when they are equivalent to the Riemann-Liouville derivatives. The proposed method can be enhanced using standard approaches for the classic Monte Carlo method, and it also allows easy parallelization, which means that it is of high potential for applications of the fractional calculus.


Introduction
It is surprising that, to our best knowledge, although there exists a huge amount of literature on Monte Carlo method for integration, there are no available works on using Monte Carlo method for differentiation. This paper is devoted to introducing the numerical evaluation of the Grünwald-Letnikov fractional derivatives by the Monte Carlo method. It is well known and documented that the Grünwald-Letnikov fractional derivatives play an important role in various numerical methods and approximations of the Riemann-Liouville and Caputo fractional derivatives, based on fractional-order differences, see [11] or [1,2,9,10].
In this work we, however, develop a totally different approach. After recalling the general framework and necessary notions, we introduce a stochastic interpretation of the Grünwald-Letnikov definition of fractional-order differentiation, and demonstrate, that the evaluation of a fractional-order derivative at a given point can be replaced by the study of a certain stochastic process.
Then we outline the basic scheme of the proposed Monte Carlo approach, and provide the algorithm for computations. Close look at the separate stages of this algorithm discovers an unexpected link between our Monte Carlo method for fractional differentiation, on one side, and the known finite difference methods on non-equidistant or nested grids, on the other side. In particular, in both cases the function values are evaluated at the nodes that are more dense near the current point, and less dense towards the beginning of the considered interval.
Implementation in the form of a toolbox for MATLAB allowed experiments with the functions, that are most frequently used in applications of fractional-order modeling, like Heaviside's unit step function, power function, exponential function, trigonometrical functions, and the Mittag-Leffler function. The provided examples demonstrated excellent agreement of the exact fractional-order derivatives of the considered functions and the numerical results by the proposed method.
In the concluding remarks we emphasize two important aspects. First, the proposed method can be enhanced using standard approaches for improving the classical Monte Carlo method, such as reduction of variance, importance sampling, stratified sampling, control variates or antithetic sampling. Second, the proposed method allows parallelization, which means that finally parallel algorithms can be used in the applications of the fractional calculus.

Grünwald-Letnikov fractional derivatives
Let α > 0, and let The function f (α) , defined uniquely by the uniqueness of the Fourier transform, is the Riemann-Liouville fractional derivative of f .
Similarly, let us define In order to calculate the fractional-order derivative of f , the Grünwald-Letnikov fractional derivative can be used, that was introduced in [11] as a non-local operator on L 1 (R) given by and When (2.1) is applied to function f ∈ F + 0 , we extend f to L 1 (R) by setting f (t) = 0 for t < 0. Hence the operator (2.1) can be regarded as an operator on L 1 (R + ), see [1,2,9,11] for more details.
In particular, the operator ((2.1) is well defined for a class F of bounded functions f , such that f and its derivatives of order up to n > 1 + α exist and are absolutely integrable, and its Fourier transform is (ik) αf (k), see, e.g., [9,pp. 22-23].
Applying the Stirling approximation The binomial series converges for any complex |z| ≤ 1 and any α > 0. Thus for z = 1 we have ∞ k=0 w k = (1 − 1) α = 0, and hence This has been noticed by Machado [8]. Denoting Note that [11] Using (2.4), we can obtain the following relationship: which motivates the general form of fractional-order derivative as Integration by parts gives the Caputo from of fractional derivative: which is just the regularized form of the Riemann-Liouville fractional derivative: Thus, the Grünwald-Letnikov fractional derivative can be considered as approximation of the Caputo of the Riemann-Liouville fractional derivatives in the numerical analysis of fractional differential equations. In particular, if f ∈ F (or F + ) then the non-local operator A α h f defined in (2.2) converges to RL D α f or C D α f in L 1 (R) (or L 1 (R + )) norm as h → 0, see [9,11] for more details.

Monte Carlo approach to the Grünwald-Letnikov fractional derivatives
Let Y be a discrete random variable such that for any fixed t and h. Let Y 1 , Y 2 ,…, Y n , …, are independent copies of the random variable Y ; then by the strong law of large numbers with probability one for any fixed t and h, and hence where −→ D means convergence in distributions and N (0, 1) is the standard normal law. Here Thus, for a given ε > 0 where 1− ε 2 is the quantily of the standard normal law, and hence for a large N we have with probability 1 − ε the following asymptotic confidence interval: for example, if ε = 0.05, 1− ε 2 ≈ 1.96.

The basic scheme of the method
The above results can be used as the basis of the Monte Carlo method for numerical approximation and computation of the Grünwald-Letnikov fractional derivatives. Indeed, we can replace the samples Y 1 , Y 2 , …, Y N by their Monte Carlo simulations. For the simulation of the random variable Y with distribution (3.1), we introduce the cumulative distribution function If U is a random variable uniformly distributed on [0, 1], then and hence, to generate Y ∈ {1, 2, . . .}, we set (4.1) Each trial (draw) of the proposed Monte Carlo method includes the following steps.

Evaluate
3. Generate N independent uniformly distributed random points, and compute the values Y i using (4.1).

Evaluate the expression
After repeating steps 1-4 K times (trials), the mean of the obtained K values gives an approximation of the value of the fractional derivative of order α, 0 < α ≤ 1, at point t.

Close look at the stages of implementation
The proposed algorithm has been implemented in MATLAB [13]; this allows useful visualizations and numerical experiments with the functions that are frequently used in the fractional calculus and in fractional-order differential equations.
Obviously, the approximation (4. As a result, the points t k , at which the function f (t) is evaluated, are distributed over the interval [0, t] non-uniformly. In Fig. 2, Fig. 3, and Fig. 4 are shown examples of   The overall picture reminds some recent efforts in the finite difference methods based on non-equidistant grids or nested meshes [3][4][5][6][7]14]. The common feature is represented by the higher density of discretization nodes near the current point, and the lesser density far from the current point; see [4, Fig. 1] and [7,Fig 2]. However, in the aforementioned papers the function values at the discretization nodes are taken with different weights, while in the proposed Monte Carlo approach all weights are the same.
The dependence of the density of nodes on α indicates that the finite difference approaches using non-uniform grids or nested meshes can be improved, if this dependence on the order α of fractional differentiation would eventually be taken into account.

Examples
The following examples demonstrate the use of the proposed Monte Carlo method for fractional-order differentiation. In all examples, the results of the computations are compared with the exact fractional-order derivatives of the function under fractional-   , z ∈ C, α,β > 0, that appears in some of the provided examples, is computed using [12]. In all examples the considered interval is sufficiently large, namely t ∈ [0, 10]. The exact fractional derivatives are plotted using solid lines, the results of the proposed Monte Carlo method are shown by bold points, the results of K individual trials (draws) are shown by vertically oriented small points (in all examples, K = 100), and the confidence intervals are shown by short horizontal lines above and below the bold points.
The result for α = 0.2 is shown in Fig. 5, and for α = 0.5 in Fig. 6. Since H (t) = 1 for t > 0, these figures represent, in fact, the Grünwald-Letnikov (and Riemann-Liouville) fractional-order derivative of a constant. We see that in this case the sample variance of the values obtained in individual trials (draws) is very small.
The result for ν = 1.3 and α = 0.5 is shown in Fig. 7. In this case, the sample variance of the values obtained in individual trials (draws) increases, but the confidence intervals remain sufficiently small.
The result for λ = 0.1 and α = 0.5 is shown in Fig. 8. The sample variance of the values obtained in individual trials (draws) increases, but the confidence intervals remain sufficiently small.
The result for λ = 0.1 and α = 0.5 is shown in Fig. 9. We see that in this case the sample variance of the values obtained in individual trials (draws) is very small, and the same holds for the confidence intervals.
The result for λ = 0.1 and α = 0.7 is shown in Fig. 10. We see that in this case, while the sample variance of the values obtained in individual trials (draws) is large enough, the confidence intervals are still sufficiently small.

Example 5. The trigonometric functions
Taking into account that sin(t) = t E 2,2 (−t 2 ) and cos(t) = E 2,1 (−t 2 ), we have: The results for (A) and (B) are shown in Fig. 11 and Fig. 12, respectively. We see that in this case also the sample variance of the values obtained in individual trials (draws) is very small, and the same holds for the confidence intervals.
Overall, these examples show that the proposed Monte Carlo method for fractionalorder differentiation works well for various kinds of functions that are important for applications of the fractional calculus, and for functions of various kinds of behavior.

Concluding remarks: a way to parallelization
In this work, the Monte Carlo method is proposed for approximation and computation of fractional-order derivatives. It can be used for evaluation of all three types of fractional-order derivatives, that usually appear in applications: the Grünwald-Letnikov, the Riemann-Liouville, and also the Caputo fractional derivatives, when they are equivalent to the Riemann-Liouville derivatives [11,Section 3.1].
The proposed method is implemented in the form of a toolbox for MATLAB, and illustrated on several examples. This opens a way to development of a family of Monte Carlo methods for the fractional calculus, using standard methods for enhancements, such as reduction of variance, importance sampling, stratified sampling, control variates or antithetic sampling.
By its nature, the proposed Monte Carlo method for fractional-order differentiation allows parallelization of computations on multiple core processors, GPUs, computer grids, and on parallel computers, and therefore have high potential for applications. by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.