# -*- 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)