An algorithm to numerically invert functions in the Laplace field is presented. It is based on the Fast Fourier Transform (FFT) technique and yields a numerical solution for t=a ("a" is a real number) for a Laplace function F(s) = L(f(t)), where "L" represents the Laplace transformation.
Download
Copy to clipboard1 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 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 | #################################################################
# 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)
|
Linear differential equations at initial values, either ordinary or partial, are often found in engineering's classical problems. For example, flow of fluids through porous media, electronic circuits, heat conduction in solids, etc, are phenomema that are described by differential equations. The Laplace transformation is a technique that can be utilised to solve these equations by transforming them into equations in the Laplace domain, where they can be more easily manipulated and eventually inverted to yield the solution in the original domain.
The analytical inversion may sometimes prove to be either nasty, tedious or impossible, though. When all classic methods fail, a numerical inversion can be the last resort available to render a solution.
There are several methods available. Here the one based on the Fast Fourier Transform method is presented. It is reasonably fast, accurate and easy to use.
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)
Hi,
I have used the code with success. Thanks.
Best Regards, Juan