################################################################# # Function Zakian(t), numerically inverts a Laplace transform # # F(s) into f(t) for a specific time "t" using the Zakian # # algorithm. # # # # Function F(s) is defined in separate as Fs(s) (see code # # below). In the present case, F(s) = 1/(s-1), whose # # analytical inversion is f(t) = exp(t), and t=1.0 thus # # yielding f(1) = 2.71828180837 what is "e" with an error # # of -7.4e-007%. Oscillatory systems may become innacurate # # after the second cycle, though. # # # # Fs(s) has to be changed accordingly everytime the user # # wants to invert a different function. # # # # Creator: Fausto Arinos de Almeida Barbuto (Calgary, Canada) # # Date: May 16, 2002 # # E-mail: fausto_barbuto@yahoo.ca # # # # Reference: # # Huddleston, T. and Byrne, P: "Numerical Inversion of # # Laplace Transforms", University of South Alabama, April # # 1999 (found at http://www.eng.usouthal.edu/huddleston/ # # SoftwareSupport/Download/Inversion99.doc) # # # # Usage: just invoke Zakian(t), where t>0. # # # ################################################################# # We need cmath because F(s) is a function operating on the # complex argument s = a + bj from cmath import * def Zakian(t): # Tupple "a", of five complex members. a = 12.83767675+1.666063445j, 12.22613209+5.012718792j,\ 10.93430308+8.409673116j, 8.776434715+11.92185389j,\ 5.225453361+15.72952905j # Tupple "K", of five complex members. K = -36902.08210+196990.4257j, 61277.02524-95408.62551j,\ -28916.56288+18169.18531j, +4655.361138-1.901528642j,\ -118.7414011-141.3036911j sum = 0.0 # Zakian's method does not work for t=0. Check that out. if t == 0: print "\n" print "ERROR: Inverse transform can not be calculated for t=0" print "WARNING: Routine zakian() exiting. \n" return ("Error"); # The for-loop below is the heart of Zakian's Inversion Algorithm. for j in range(0,5): sum = sum + (K[j]*Fs(a[j]/t)).real return (2.0*sum/t) # The Laplace function F(s) is defined here. def Fs(s): return 1.0/(s-1.0) # Function Zakian(t) is invoked. print Zakian(1.0)