[R] Re: Inverse of the Laplace Transform/Gaver Stehfest algorithm

[Tolga Uzuner]

Tolga Uzuner wrote:

> Hi there,
>
> Is there an implementation of the Gaveh Stehfest algorithm in R 
> somewhere ? Or some other inversion ?
>
> Thanks,
> Tolga
>
Well, at least here is Zakian's algorithm, for anyone who needs it:

Zakian<-function(Fs,t){
# Fs is the function to be inverted and evaluated at t
    a = c(12.83767675+1.666063445i, 
12.22613209+5.012718792i,10.93430308+8.409673116i, 
8.776434715+11.92185389i,5.225453361+15.72952905i)
    K = c(-36902.08210+196990.4257i, 
61277.02524-95408.62551i,-28916.56288+18169.18531i, 
+4655.361138-1.901528642i,-118.7414011-141.3036911i)
    ssum = 0.0
#   Zakian's method does not work for t=0. Check that out.
    if(t == 0){
        print("ERROR:   Inverse transform can not be calculated for t=0")
        print("WARNING: Routine zakian() exiting.")
        return("Error")}
#   The for-loop below is the heart of Zakian's Inversion Algorithm.
    for(j in 1:5){ssum = ssum + Re(K[j]*Fs(a[j]/t))}
 
    return (2.0*ssum/t)
}

#   InvLap(1/(s-1))=exp(t)
#   check if Zakian(function(s){1/(s-1)},1)==exp(1)
lapfunc<-function(s){1.0/(s-1.0)}

#   Function Zakian(functobeinverted,t) is invoked.
Zakian(lapfunc,1.0)