#################################################################
# Function InvLap(t,omega,sigma,nint), numerically inverts a #
# Laplace transform F(s) into f(t) using the Fast Fourier #
# Transform (FFT) algorithm for a specific time "t", an #
# upper frequency limit "omega", a real parameter "sigma" #
# and the number of integration intervals "nint" . #
# #
# Function F(s) is defined in separate as Fs(s) (see code #
# below). Fs(s) has to be changed accordingly everytime the #
# user wants to invert a different function. #
# #
# I suggest to use omega>100 and nint=50*omega. The higher #
# the values for omega, the more accurate the results will be #
# in general, but at the expense of longer processing times. #
# #
# Sigma is a real number which must be a little bigger than #
# the real part of rightmost pole of the function F(s). For #
# example, F(s) = 1/s + 1/(s-2) + 1/(s+1) has poles for s=0, #
# s=2 and s=-1. Hence, sigma must be made equal to, say, #
# 2.05 so as to keep all poles at the left of this value. #
# The analytical inverse for this simple function is #
# f(t) = 1 + exp(-t) + exp(2t). For t=1.25, omega=200, #
# nint=10000 and sigma=2.05, the numerical inversion yields #
# f(1.25) ~= 13.456844516, or -0.09% away from the actual #
# analytical result, 13.468998757 (results truncated to 9 #
# decimal places). This is the example used in this code. #
# #
# Creator: Fausto Arinos de Almeida Barbuto (Calgary, Canada) #
# Date: May 18, 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: invoke InvLap(t,omega,sigma,nint), for t>0. #
# #
#################################################################
# We need cmath because F(s) is a function operating on the
# complex argument s = a + bj
from math import ceil
from cmath import *
# *** Driver InvLap function ***
def InvLap(t,omega,sigma,nint):
# Sanity check on some parameters.
omega = ceil(omega)
nint = ceil(nint)
if omega <= 0:
omega = 200
if nint <= 0:
nint = 10000
return (trapezoid(t,omega,sigma,nint))
# *** Function trapezoid computes the numerical inversion. ***
def trapezoid(t,omega,sigma,nint):
sum = 0.0
delta = float(omega)/nint
wi = 0.0
# The for-loop below computes the FFT Inversion Algorithm.
# It is in fact the trapezoidal rule for numerical integration.
for i in range(1,(nint+1)):
witi = complex(0,wi*t)
wf = wi + delta
wfti = complex(0,wf*t)
fi = (exp(witi)*Fs(complex(sigma,wi))).real
ff = (exp(wfti)*Fs(complex(sigma,wf))).real
sum = sum + 0.5*(wf-wi)*(fi+ff)
wi = wf
return ((sum*exp(sigma*t)/pi).real)
# *** The Laplace function F(s) is defined here. ***
def Fs(s):
return (1.0/s + 1.0/(s+1.0) + 1.0/(s-2.0))
# Function InvLap(t,omega,sigma,nint) is invoked.
print InvLap(1.25,200,2.05,10000)