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

[Tolga Uzuner]

Tolga Uzuner wrote:

> 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)
>
>
By the way, I am told "This significance of this specification is that 
Zakian’s Algorithm is very accurate for overdamped and slightly 
underdamped systems. But it is not accurate for systems with prolonged 
oscillations."

for those considering using it.