Welcome, guest | Sign In | My Account | Store | Cart
```#################################################################
#   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/      #
#                                                               #
#   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)
```