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.
______________________________________________
[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
