Message-ID: <42614329.90204@coubros.com>
Date: 2005-04-16T16:54:01Z
From: Tolga Uzuner
Subject: Inverse of the Laplace Transform/Gaver Stehfest algorithm
In-Reply-To: <425EFDE6.5040204@coubros.com>
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)