This is a fast and highly accurate numerical method for the inversion of the Laplace transform
Python, 65 lines
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 | # -*- coding: utf-8 -*-
# Talbot suggested that the Bromwich line be deformed into a contour that begins
# and ends in the left half plane, i.e., z → −∞ at both ends.
# Due to the exponential factor the integrand decays rapidly
# on such a contour. In such situations the trapezoidal rule converge
# extraordinarily rapidly.
# For example here we compute the inverse transform of F(s) = 1/(s+1) at t = 1
#
# >>> error = Talbot(1,24)-exp(-1)
# >>> error
# (3.3306690738754696e-015+0j)
#
# Talbot method is very powerful here we see an error of 3.3e-015
# with only 24 function evaluations
#
# Created by Fernando Damian Nieuwveldt
# email:fdnieuwveldt@gmail.com
# Date : 25 October 2009
#
# Reference
# L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer. Talbot quadratures
# and rational approximations. BIT. Numerical Mathematics,
# 46(3):653 670, 2006.
from cmath import *
def Talbot(t,N):
# Initiate the stepsize
h = 2*pi/N;
# Shift contour to the right in case there is a pole on the positive real axis : Note the contour will
# not be optimal since it was originally devoloped for function with
# singularities on the negative real axis
# For example take F(s) = 1/(s-1), it has a pole at s = 1, the contour needs to be shifted with one
# unit, i.e shift = 1. But in the test example no shifting is necessary
shift = 0.0;
ans = 0.0;
if t == 0:
print "ERROR: Inverse transform can not be calculated for t=0"
return ("Error");
# The for loop is evaluating the Laplace inversion at each point theta which is based on the trapezoidal rule
for k in range(0,N):
theta = -pi + (k+1./2)*h;
z = shift + N/t*(0.5017*theta*cot(0.6407*theta) - 0.6122 + 0.2645j*theta);
dz = N/t*(-0.5017*0.6407*theta*(csc(0.6407*theta)**2)+0.5017*cot(0.6407*theta)+0.2645j);
ans = ans + exp(z*t)*F(z)*dz;
return ((h/(2j*pi))*ans).real
# Here is the Laplace function to be inverted, should be changed manually
def F(s):
return 1.0/(s+1.0)
# Define the trig functions cot(phi) and csc(phi)
def cot(phi):
return 1.0/tan(phi)
def csc(phi):
return 1.0/sin(phi)
|
Many problems in applied mathematics can be solved using the Laplace transform such as PDEs for example. The method describe here is fast and accurate. It is based on a deformation of the Bromwich line to a contour that ends in the left half plane, i.e both end points tend to infinity. Here the test function F(s) = 1/(s+1) is used. With only 24 function evaluations, the Talbot method agrees to 15 significant digits. Since it is based on the midpoint rule the number of function evaluations should be manually selected.
Minor point - common programming practice would iterate on range(0,N) and let theta = -pi+(k+1./2)/h
Ahh - semicolons at the end of lines. Old C habit.
BTW I once met the primary author of the book you referenced, Dr. Lloyd Trefethan, back when he & I were still in graduate school. I was impressed and made it a point to remember his name. Super, super smart.
A beautiful programm. The algorithm seems to be important and was unkonown to me. In this version there are, I think, rounding problems with large N. These problems can be overcome by using mpmath.
A variant of this recipe is in Recipe 576938