[R] Re: Inverse of the Laplace Transform/Gaver Stehfest algorithm
Tolga Uzuner
tolga at coubros.com
Sat Apr 16 18:57:09 CEST 2005
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.
More information about the R-help
mailing list