A simple piece of code in Python that inverts function in the Laplace field to the real field. Accurate, fast and easy to use.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 | ```
#################################################################
# 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)
``` |

Many engineering problems, especially the ones expressed by linear differential equations (either ordinary or partial) at initial values can be solved via the Laplace transforms. However, it is sometimes complicate, tedious or even impossible to invert some transforms to the realm of the real numbers. In these cases, a numerical inversion can be a life saver. Petroleum reservoir engineering and mechanics of rocks usually make use of numerical inversion methods like the one presented. Zakian's algorithm is easy to implement and use -- and produces amazingly good results.