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)
______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
