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

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 f(t) is defined as a contour integral

ft=12π<!-- p -->i∫<!-- ? -->γ<!-- ? -->−<!-- - -->i∞<!-- 8 -->γ<!-- ? -->+i∞<!-- 8 -->estFsds,$$\begin{array}{} \displaystyle f\left(t\right)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} e^{st}F\left(s\right)ds, \end{array}$$

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.

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.

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

er=1−<!-- - -->vcve.$$\begin{array}{} \displaystyle e_{r}=\left|1-\frac{v_{c}}{v_{e}}\right|. \end{array}$$

The Laplace transform and its inversion (1) is

f^<!-- ^ -->t=L−<!-- - -->1π<!-- p -->e−<!-- - -->1/ss=cos2tt,t∈<!-- ? -->0,6.$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{\sqrt{\pi}\;e^{-1/s}}{\sqrt{s}}\right\}=\frac{\cos\left(2\sqrt{t}\right)}{\sqrt{t}},\; t \in \left(0,6\right \gt . \end{array}$$(1)

Fig. 1

The Laplace transform and its inversion (2) is

f^<!-- ^ -->t=L−<!-- - -->11s+1s+2=e−<!-- - -->t−<!-- - -->e−<!-- - -->2t,t∈<!-- ? -->0,10.$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{1}{\left(s+1\right)\left(s+2\right)}\right\}=e^{-t}-e^{-2t},\; t \in \left(0,10\right \gt . \end{array}$$(2)

Fig. 2

The Laplace transform and its inversion (3) is

f^<!-- ^ -->t=L−<!-- - -->1sinb⋅<!-- ? -->s+a⋅<!-- ? -->cosbs2+a2=sina⋅<!-- ? -->t+b,a=1,b=1,t∈<!-- ? -->0,10.$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{\sin\left(b\right)\cdot s + a \cdot \cos\left(b\right)}{\left(s^2+a^2\right)}\right\}=\sin\left(a\cdot t + b\right),\; a=1,b=1,\; t \in \left(0,10\right \gt . \end{array}$$(3)

Fig. 3

The Laplace transform and its inversion (4) is

f^<!-- ^ -->t=L−<!-- - -->11s2+a2=J0at,a=1,t∈<!-- ? -->0,50,$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathscr{L}^{-1}\left\{\frac{1}{\sqrt{s^2+a^2}}\right\}=\mathscr{J}_{0}\left(at\right),\; a=1,\; t \in \left(0,50\right \gt , \end{array}$$(4)

in which 𝓙₀(t) is the Bessel function of the first kind of order 0.

Fig. 4

The Laplace transform and its inversion (5) is

f^<!-- ^ -->t=L−<!-- - -->11s+as+b=e−<!-- - -->a+bt/2I0a−<!-- - -->b2t,a=2,b=1,t∈<!-- ? -->0,30,$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathscr{L}^{-1}\left\{\frac{1}{\sqrt{s+a} \sqrt{s+b}}\right\}=e^{-\left(a+b\right)t/2} \; \mathscr{I}_{0}\left(\frac{a-b}{2}\; t\right),\;a=2,b=1,\; t \in \left(0,30\right \gt , \end{array}$$(5)

in which 𝓘₀(t) is modified Bessel function of order 0.

Fig. 5

The Laplace transform and its inversion (6) is

f^<!-- ^ -->t=L−<!-- - -->11s⋅<!-- ? -->s+a=ea2t⋅<!-- ? -->erfcat,a=1,t∈<!-- ? -->0,30,$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{1}{\sqrt{s}\cdot\left(\sqrt{s}+a\right)}\right\}=e^{a^2t}\cdot erfc\left(a\sqrt{t}\right),\;a=1,\; t \in \left(0,30\right \gt , \end{array}$$(6)

in which erfc is the complementary error function, 1 − erf(x).

Fig. 6

The Laplace transform and its inversion (7) is

f^<!-- ^ -->t=L−<!-- - -->1sα<!-- a -->−<!-- - -->1sα<!-- a -->+λ<!-- ? -->=Eα<!-- a -->−<!-- - -->λ<!-- ? -->tα<!-- a -->,α<!-- a -->=0.5,λ<!-- ? -->=1,t∈<!-- ? -->0,20$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{s^{\alpha-1}}{s^{\alpha}+\lambda}\right\}=\mathcal{E}_{\alpha}\left(-\lambda t^{\alpha}\right),\; \alpha=0.5,\lambda=1,\; t \in \left(0,20\right \gt \end{array}$$(7)

in which

Eα<!-- a -->z=∑<!-- ? -->k=0∞<!-- 8 -->zkΓ<!-- G -->α<!-- a -->k+1,Rα<!-- a -->><!-- > -->0,z∈<!-- ? -->C$$\begin{array}{} \displaystyle \mathcal{E}_{\alpha}\left(z\right)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma\left(\alpha k+1\right)} ,\,\mathcal{R}\left(\alpha\right) \gt 0\,,z\in \mathcal{C} \end{array}$$

is the Mittag-Leffler function with one parameter [43].

Fig. 7

The Laplace transform and its inversion (8) is

f^<!-- ^ -->t=L−<!-- - -->1sα<!-- a -->−<!-- - -->β<!-- ? -->sα<!-- a -->+λ<!-- ? -->=tβ<!-- ? -->−<!-- - -->1Eα<!-- a -->,β<!-- ? -->−<!-- - -->λ<!-- ? -->tα<!-- a -->,α<!-- a -->=0.9,β<!-- ? -->=0.1,λ<!-- ? -->=1,t∈<!-- ? -->0,20,$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{s^{\alpha-\beta}}{s^{\alpha}+\lambda}\right\}=t^{\beta-1} \mathcal{E}_{\alpha,\beta}\left(-\lambda t^{\alpha}\right),\; \alpha=0.9,\beta=0.1,\lambda=1,\; t \in \left(0,20\right \gt , \end{array}$$(8)

in which

Eα<!-- a -->,β<!-- ? -->z=∑<!-- ? -->k=0∞<!-- 8 -->zkΓ<!-- G -->β<!-- ? -->+α<!-- a -->k,Rα<!-- a -->><!-- > -->0,Rβ<!-- ? -->><!-- > -->0,z∈<!-- ? -->C$$\begin{array}{} \displaystyle \mathcal{E}_{\alpha,\beta}\left(z\right)=\sum_{k=0}^{\infty}\frac{z^{k}}{\Gamma\left(\beta+\alpha k\right)},\, \mathcal{R}\left(\alpha\right) \gt 0,\, \mathcal{R}\left(\beta\right) \gt 0 \,,z\in \mathcal{C} \end{array}$$

is generalized the Mittag-Leffler function with two parameters [44].

Fig. 8

The Laplace transform and its inversion (9) is

f^<!-- ^ -->t=L−<!-- - -->1sα<!-- a -->−<!-- - -->1sα<!-- a -->+λ<!-- ? -->=tα<!-- a -->−<!-- - -->1Γ<!-- G -->α<!-- a -->e−<!-- - -->λ<!-- ? -->t,α<!-- a -->=0.5,λ<!-- ? -->=1,t∈<!-- ? -->0,20.$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{s^{\alpha-1}}{s^{\alpha}+\lambda}\right\}=\frac{t^{\alpha-1}}{\Gamma\left(\alpha\right)}e^{-\lambda t} ,\; \alpha=0.5,\lambda=1,\; t \in \left(0,20\right \gt . \end{array}$$(9)

Fig. 9

The Laplace transform and its inversion (10) is

f^<!-- ^ -->t=L−<!-- - -->11sα<!-- a -->=tα<!-- a -->−<!-- - -->1Γ<!-- G -->α<!-- a -->,α<!-- a -->=0.5,t∈<!-- ? -->0,20.$$\begin{array}{} \displaystyle \hat{\text{f}}\left(t\right)=\mathcal{L}^{-1}\left\{\frac{1}{s^{\alpha}}\right\}=\frac{t^{\alpha-1}}{\Gamma\left(\alpha\right)} ,\; \alpha=0.5,\; t \in \left(0,20\right \gt . \end{array}$$(10)

Fig. 10

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.