Review of numerical methods for NumILPT with computational accuracy assessment for fractional calculus

Abstract In the paper we present results of accuracy evaluation of numerous numerical algorithms for the numerical approximation of the Inverse Laplace Transform. The selected algorithms represent diverse lines of approach to this problem and include methods by Stehfest, Abate and Whitt, Vlach and Singhai, De Hoog, Talbot, Zakian and a one in which the FFT is applied for the Fourier series convergence acceleration. We use C++ and Python languages with arbitrary precision mathematical libraries to address some crucial issues of numerical implementation. The test set includes Laplace transforms considered as difficult to compute as well as some others commonly applied in fractional calculus. Evaluation results enable to conclude that the Talbot method which involves deformed Bromwich contour integration, the De Hoog and the Abate and Whitt methods using Fourier series expansion with accelerated convergence can be assumed as general purpose high-accuracy algorithms. They can be applied to a wide variety of inversion problems.


Introduction
Fractional calculus (FC) can be interpreted as an extension of the concept of derivative operator from integer order n to arbitrary order ν, where n is an integer and ν can be a real number. English term fractional calculus is misleading because it suggests that differentiation and integration order may assume non-integer orders only. The more appropriate description of this branch of mathematics would be differentiation and integration of any order. In practice, both integration as well as differentiation orders may assume real or complex values. In this context derivatives and integrals of integer order can be considered as one of special cases of fractional calculus [1][2][3][4].
FC can be applied for finding solutions of differential equations of fractional order (abbreviation FODE). For this reason, the research for a numerical method for finding their accurate and stable solutions has been set for the last 70 years within the scope of many scientific projects as for example [5][6][7][8][9][10].
The Laplace Transform (abbreviation LPT) is commonly applied for solving linear differential equations of integer order by converting them into algebraic equations, which, in this form, are simpler to solve [11][12][13]. However, a found solution has to be inverted by applying the Inverse Laplace Transform (abbreviation ILPT). This inversion is commonly done by using tables with so called the Laplace pairs, consisting of a Laplace transform and an corresponding original. If an original is not found, the numerical approximation of ILPT (abbreviation NumILPT) can be applied to find it.
NumILPT is an ill-posed inversion problem due to the inclusion of the multiplication of an exponential function of time in the inversion formula, in which algorithmic and finite precision errors can lead to exponential divergence of numerical solutions [14]. Due to this difficulty, one has been trying to find a method to mitigate its influence on final accuracy of numerically inverted ODE solutions for the last 70 years (as it is stated in [15] the first recognized algorithm of NumILPT is the algorithm by D.V. Widder dated from 1946 [16]). As a result, currently there can be distinguished not only several main numerical approaches to the inversion problem, but also within a single approach there can be noted numerous different approximation methods to the various parts of the inversion formula.
Main approaches to the NumILPT include the use of: Bromwich contour integration, Fourier series expansion, finite differences and Laguerre and Legendre polynomials.
Bromwich contour can have a fixed shape or, to obtain better results, it can be deformed. The deformation can be done by using cotangent function or any other trigonometric function depending on which deformation is then more accurately integrated by a selected method of integration, e.g. the trapezoidal rule.
The exponential part of the inversion formula can be approximated by applying Taylor series expansion or Padée approximation. For the same purpose there can also be used Fourier series expansion, within which the exponential part is approximated with the use of a Padée approximant.
In this context, the main task within the scope of this paper is present the results of the evaluation of some of these algorithms for the inversion accuracy of Laplace transforms frequently used in FC.
However, it must be stressed that the accuracy evaluation is not the most difficult part of this task, but the right selection of the algorithms according to the aim of the scope of the evaluation (numerical inversion difficulty level) and their efficient programming which must include solutions to the implementation issues associated, among the others, with the limitations of the double precision computer arithmetic.
In the paper [17] we gave an initial evaluation for the inversion accuracy of several methods representing the main approaches to the NumILPT. Testbed included Laplace transforms of FODE solutions we have found in the most cited papers at those times. Currently, we would like to present results of a new evaluation inspired by the lecture of [18] and [19] publishing an impressive list of Laplace transforms used in the very practical applications of FC. The new results were obtained by using improved computer implementations of selected algorithms of NumILPT using C++ and Python programming languages.
The paper is divided in the following sections. The first section includes short information about the Inverse Laplace Transform, its numerical counterpart -the numerical approximation of the Inverse Laplace Transform. It list all the methods used in the evaluation and gives briefly insights into their mathematical algorithms as well. The third subsection discuses the most important problems of theor programming and presents suggested solutions to the most severe of them. Sections four and five presents Laplace transforms criteria selection and how computational accuracy evaluation was conducted. Section six contained the presentation of the results in form of two plots with each evaluated method for each selected Laplace transform. The paper end with usual conclusions and extended list of references.

The Inverse Laplace Transform
The Inverse Laplace Transform f (t) is defined as a contour integral where all the singularities of F (s) are to the left of the vertical line γ − i∞, γ + i∞.
The analytic Inverse Laplace Transform involves application of the well known theory of complex variables. For isolated singularities, the Bromwich contour is the standard approach.

The Selection of Algorithms for Numerical Approximation of the Inverse Laplace Transform
If there is no analytical inversion formula, numerical approach to the problem is required. This is done by using numerical Inverse Laplace Transform, denoted asf (t).
Numerical inversion of Laplace transforms is not an easy task [20]. There neither exists an universal method of accurate inversion nor an inversion can be always conducted with satisfactory accuracy. Numerical inversion can be done using numerous distinct algorithms and only a few of them engage numerical integration of the contour. In the paper we evaluate for computational accuracy numerous methods representing all major lines of approach to the numerical inversion.
The Widder formula, mentioned in the introduction, required high-order analytic image function derivatives and it was thus impractical for computational purposes. The method of Stehfest is a discrete version of this formula. Another popular method used for numerical inversion is Fourier series expansion with accelerated convergence. There also have to be noted a modification of Stehfest finite differences evaluation method known as Gaver-Stehfest algorithm [21]. Bilinear Inverse Laplace Transform is possible by the use of Laguerre [22] and Legendre polynomials [23][24][25]. Delta function approximation by a finite linear combination of exponential functions was also used, as well as the approximation of a part of the Inverse Laplace Transform by a rational function. The method which follows closely the analytical approach to the numerical inversion is the use of deformed or fixed Bromwich contour integration using the Trapezoidal rule.
Next we present a briefly describe some of the methods.

Stehfest Method
The Stehfest inversion formula [26] is based on computing a sample of the time function by using finite differences and Salzer summation algorithm. The approximation assumes the following form: where the weight coefficients K n are given by The available parameters are the choice of N (which must be even) and the values of s at which F (s) tdepended must be sampled.

Abate and Whitt Method
An algorithm involves Fourier series was originally proposed by Dubner and Abate [27] and improved by Abate and Whitt [28]. The algorithm employs Euler summation formula [29] (van Wijngaarden algorithm), which is shown to accelerate convergence of a slowly converging truncated Fourier series [30].
It is essentially the Trapezoidal rule applied to Laplace inversion formula. By applying the Trapezoidal rule with step size h and setting h = π/ (2t) and a = A/ (2t), we get The choice of A has to be made in such a way that a falls at the left of the real part of all singularities of the inverted Laplace transform.

De Hoog Method
It is an improved procedure for numerical inversion of Laplace transforms. It is based on accelerating the convergence of the Fourier series obtained from the inversion of integral by using non-linear double acceleration with Padé approximation and an analytic expression for the remainder in the series [31]. Non-linear acceleration techniques drastically reduce the required number of functions evaluations [32].
Speed of algorithm is primarily a function of the order of the Taylor series expansion, which is controlled manually.

FFT Method
The non-accelerated Fourier series inverse algorithm is almost useless because it requires thousands of F(s) evaluations. Practical approaches accelerate the convergence of the sum by using for example well known FFT (Fast Fourier Transform) algorithm [33].
The program inverts numerically a Laplace transform f(s) intof(t) using the Fast Fourier Transform (FFT) algorithm for a specific time t, an upper frequency limit ω, a real parameter σ and the number of integration intervals N.
Parameter σ contains singularities coordinates of F(s) and must be input by a user.

Vlach and Singhai Method
The basic inversion formula which works for real time functions assumes the following form [34] where K i are coefficients obtained as a result of Padée approximation of exponential part of the Laplace inversion formula.

Zakian Method
Zakian developed a method [35], in which approximation for f (t) is given bŷ where sets of α i K i can be obtained as solution to The only free parameter of choice is N. Its value should be 10 for best results.

Talbot Method
Talbot developed a method for the numerical inversion of the Laplace transform, in which the inversion is approximated by the Trapezoidal rule along a special deformed contour.
The Inverse Laplace Transform formula can be rewritten aŝ where B is the Bromwich contour from r − i∞ to r + i∞ and t > 0, so that B is to the right of all singularities of F (s).
Talbot provided [36] a recipe for a contour, in which the contour is moved to the left so as to reduce in magnitude the factor e st in the integrand, but the contour must not be moved too close to singularities of F (s), as to do so will result in peaks in the integrand function. This requires to be known the locations of the singularities of F (s) [37].
Weideman [38] proposed an improved formula involving the cotangent function with additional parameters.
3 Some Aspects of Programming

Numerical Inversion Issues
The numerical Inverse Laplace Transform is an ill-posed problem by inherent sensitivity due to the multiplication by an exponential function of time: • Finite precision and algorithmic errors can lead to the exponential divergence of numerical solutions, • High precision variables are required for most Laplace Transform inversion methods: double precision variable can hold the values from the range 10 −308 − 10 308 . Value 10 308 ∼ e 709 is large enough for only modest t.
To counteract possible issues we apply some modern programming tools and techniques briefly described in the next two subsections.

Modern Programming Languages
C++ and Python are modern languages which have similar capabilities. However, Python achieves a higher abstract level than C++, i.e. an individual programmer can achieve the same results that in a much shorter time and with far fewer lines of code using it. It also has especially clean and straightforward syntax. It can lead to programs' shorter executing time and lower computational complexity.
Still, both languages' most respectful inventory in context of the goals of the following research are the arbitrary precision libraries: GNU GMP/MPFR for C++ and mpmath for Python.

Advantages of Arbitrary Precision Libraries
The GNU Multiple Precision Floating-Point Reliable Library (MPFR) [39] is an arbitrary precision package for C language and it is based on GNU Multiple-Precision Library (GMP) [40]. MPFR supports arbitrary precision floating point variables. It also provides proper rounding of all implemented operations and mathematical functions.
The mpmath [41] library is a free Python library for real and complex floating-point arithmetic with arbitrary precision.
The main advantages of the libraries are: • Almost any calculation can be performed just as well at 10-digit or 1000-digit precision, with either real or complex numbers, and in many cases they implement efficient algorithms that scale well for extremely high precision work, • Proper rounding in compliance with IEEE 754-2008 standard, • A high number of special functions, with arbitrary precision and full support for complex numbers, • Rudimentary support for interval arithmetic.
The libraries enable a user to set the precision of the arbitrary precision variables by specifying the number of bits to use in the mantissa of the floating point number. Due to the design of the libraries it is possible to work with any precision between 2 bits and maximum bits allowed for a computer.
The most common errors in numerical calculations are caused by wrong rounding. The libraries support proper rounding in compliance with the IEEE 754-2008 standard, as well as an additional one not included in it, i.e. round away from 0.
An increased computational time complexity for arbitrary precision calculations, which is their main disadvantage, is caused not only by applying variables of increased precision (more bits for mantissa in the floating point number), but also by applying the corresponding functions and different syntax, which is required. However, if an application is speed-depended -so called wrapper, as for example MPFR C++ [42] can be used. It enables using arbitrary precision in connection with standard C++ syntax. This simplifies programming and does not increases computational time complexity. For example, computational time complexity using this wrapper for up to 100-digits precision computing does not increase it more than it is in case of double precision application done on mid-2000 laptop.

Laplace Transforms Criteria Selection
We spent a large amount of time studying papers devoted to numerical Inverse Laplace Transform problem and solutions of fractional order differential equations. To be able to select aptly Laplace transforms which are useful in both areas, we consulted numerous noble mathematicians actively involved in FODE about the subject. As a result, we formed a list with Laplace transforms, which are the most frequently used and which represent various kinds of inversion problems. The list also included transforms which represent generally difficult computational problems.

Accuracy Criterion
We have considered various accuracy assessment criteria of numerically obtained results which are usually applied in this kind of research. They include: (a) absolute/relative error to obtain a real knowledge how close to an exact inversion an approximation really is, (b) quadratic mean and quadratic mean weighted by the factor e −t and the difference between them which give root-mean-square deviation between the analytical and the numerical calculated values for small and large t and between them respectively, (c) order of convergence which informs how a method performs in respect to increasing number of sampling points, (d) time and memory complexity of computations which shows how a method occupies a computer with increasing amount of input data.
Due to the main goal of the numerical experiment, we selected relative error measure (a) for accuracy evaluation. It is commonly applied in similar to our comparisons. It actually can also be considered as (b), because it informs well how accurate to an exact inversion an approximation at each t is. However, due to the form of results presentation we selected method (a) because advantages of (b) can be utilized for results presented in tables.
The criterion (c) is abandoned due to the irrelevancy to the goals of the experiment. Order of convergence for each method can be found in source papers. They are presented in the references section.
The only drawback of (b) is an increasing value of relative error in case of decreasing values of exact inversion. It gives a false conclusion about an with time deteriorating accuracy. In such situation there is advised comparison of inversion and relative error plots.
The next section presents accuracy comparison results in form of two plots for each selected Laplace transforms: the first plot of an actual approximation by applying a respective numerical inversion method and the second plot with relative error e r of an approximation at each t for each of the methods v c computed in respect to the values of an analytical inversion formula v e e r = 1 − v c v e .

Results
• The Laplace transform and its inversion (1) iŝ • The Laplace transform and its inversion (2) iŝ • The Laplace transform and its inversion (3) iŝ • The Laplace transform and its inversion (4) iŝ in which J 0 (t) is the Bessel function of the first kind of order 0. • The Laplace transform and its inversion (5) iŝ in which I 0 (t) is modified Bessel function of order 0. • The Laplace transform and its inversion (6) is in which er f c is the complementary error function, 1 − er f (x). • The Laplace transform and its inversion (7) iŝ in which is the Mittag-Leffler function with one parameter [43]. • The Laplace transform and its inversion (8) iŝ in which is generalized the Mittag-Leffler function with two parameters [44]. • The Laplace transform and its inversion (9) iŝ

Conclusions
FC is an extension of classical calculus of integer order. This fact enables finding solutions of a problem within the available numerical methods.
The solution of the problem selected in the paper focuses on the selection and the programming of seven algorithms of NumILPT. The methods represent four main approaches to the numerical approximation of the Inverse Laplace Transform. They were evaluated against a test set of inversion problems frequently used in FC and representing difficult computational issues.
Based on our previous experience and presently obtained results, we have selected three numerical algorithms of the numerical approximation of the Inverse Laplace Transform for the actual comparison: the Abate & Whitt, the DeHoog and the Talbot method for their versatility and high-accuracy inversions. Due to the general application requirement and clarity of presentation for the present evaluation, remaining methods were left aside. Their performance is either mediocre or they calculated satisfactory accurate inversions in individual cases only. It suggests their application only on an particular inversion problem basis.
NumILPT can be applied as a high-accuracy method for obtaining solutions of FODE. During the evaluation, several test problems have correctly been inverted by applying more than one algorithm. Here, the Talbot method, which uses the Trapezoidal Rule for the integration of the cotangent function deformed Bromwich contour, is the most accurate and versatile. However, it requires precise input of Laplace transform poles coordinates for high-accuracy inversions, which can be difficult and can not be automated. For the best results, manual calculation and input of the coordinates is advised.
By applying the DeHoog algorithm, which uses Fourier series expansion integration, one can obtain equally accurate inversions. Wherein, the poles coordinates input is not required.
The Abate and Whitt, another algorithm, which employs the use of Fourier series expansion approximation as well, is advised for inverting Laplace transforms, whose originals are the periodic functions. It is characterized by high correctness of inverted solutions even for t > 50.
Difficult implementation issues were successfully solved by the use of arbitrary precision mathematical libraries and concurrent of two programming languages: C++ and Python for programming.
The commonly used standard double precision computer arithmetic has many flaws and limitations, which include the limitations of the values, the double precision variables can hold, no influence of a user on the method of mathematical operations rounding or precision of the intermediate variables. The decision in these cases is made on an operating systems or compiler level.
The state of the art technology called the infinite precision computing [45] enables eliminating most of these problems by defining of arbitrary precision variables, allowing a user to access the intermediate mathematical operations rounding modes by breaking them down into basic ones and the capability of a user to explicitly choose one of five rounding modes for each mathematical operation.
For speed-depended applications, the so called wrapper MPFR C++ has been used. It enables the use of arbitrary precision variables for programming in connection with the standard C++ syntax. Therefore, its application has all of the advantages associated with the infinite precision computing and no speed punishment for the computations involving the precision up to 100 digits.
T h i s p a g e i s i n t e n t i o n a l l y l e f t b l a n k