Example notebook for the package  'FixedTalbotNumericalLaplaceInversion'

Authors

Joe Abate and Peter P.Valkó
e-mail: p-valko@tamu.edu
web:  http://pumpjack.tamu.edu/~valko/

Summary

This package provides only one function: FT. The function calculates the 
value of the inverse of a Laplace transform at a specified time point. The 
Laplace transform should be provided as a function ready for 
multiple-precision evaluation in the complex plane. In other words, 
approximate numbers (with decimal point) or Mathematica functions starting 
with the letter 'N' are not allowed. The algorithm used is called the Fixed 
Talbot algorithm.



There are two main groups of methods for the numerical inversion of the 
Laplace transform. One of them uses values of the transform calculated on the 
real line. Such an algorithm (GWR) is available from  
http://library.wolfram.com/infocenter/MathSource/4738/ . This submission 
falls into the other group of methods, requiring the evaluation of the 
transform in the complex plane. The best such method is due to Talbot and 
this realization gives an improved version of his algorithm. Deatils of the 
algorith are discussed in the paper Abate, J. and Valkó, P. P.:  
Multi-Precision Laplace Transform Inversion, International Journal of 
Numerical Methods for Engineering, (2003), accepted for publication (IJNME 
A1032.C)

Needs

The "FixedTalbotNumericalLaplaceInversion.m" file should be placed into  the AddOns folder on the user's computer. It is my experience, however, that  some users prefer to put the package file (in this case the  FixedTalbotNumericalLaplaceInversion.m file) into the same folder where the  calling notebook (in this case  FixedTalbotNumericalLaplaceInversionExample.nb) is located. The following  lines will make sure that Mathematica finds the package in the folder where the calling notebook is located:

thisDir = 
        ToFileName[
("FileName" /. 
                NotebookInformation[
                  EvaluationNotebook[]])
[[1]]];
SetDirectory[thisDir];

The package is loaded by the command:

<< FixedTalbotNumericalLaplaceInversion.m

The only function available is FT.

Help

Help is available:

?FT

Example and Notes

Our example is to invert  F = \(TraditionalForm\`e\^\(-\@\(\(s\^2 + \(37\ s\)\/100 + 1\)\/\(s\^2 + s + π\)\)\)\/s\) at time  t0 = 100.

fun[
          s_] = (1/
              s) Exp[-Sqrt[ (s^2 +  (37/100) s + 
                      1)/(s^2 +  s + Pi)] ];
t0 = 100;
ft = 
        FT[fun, t0];


Print["Laplace transform: ", 
        fun[s] // TraditionalForm];
Print["Inverse at t0: ", t0, " is  ft: ", ft];

Note 1: There is no guarantee that the result is as accurate as the precision 
returned by the Kernel. If you wish, for instance, 17-digit accuracy, 
it is suggested to double the number of terms 'M' and stop when the first 17 
digits are stabilized:

i = 3;  ft = 0; 
    Print["    M ", "  FT"];
While[
 ftold = ft; i++;  M = 2^i; 
        ft = FT[fun, t0, M]; 
 
        Print[PaddedForm[ M, 4], 
          PaddedForm[ ft, {20, 17}] ]; 
 
        Abs[ft - ftold] ≥ Abs[ft] (10^(-16))];

Note 2: Approximate numbers (decimal points) should be avoided. In the next 
example the result is nonsense:

badfun[
          s_] = (1/
              s) Exp[-Sqrt[ (s^2 +  0.37 s + 1)/(s^2 +  s + 
                      Pi)] ];
FT[badfun, t0, 128]

Note 3: The Laplace transform function 'F' should have only one argument. 
Only its name is needed in the calling statement. 

In the above examples we used the symbol 's' for the Laplace variable, but 
any other available symbol could be used. 

TheLaplace transform also can be provided as a 'pure' function:

PaddedForm[
      FT[1/(1 + #^2) &, 50, 128], {29, 29}]
N[Sin[50], 29]

Note 4: While Talbot's method is extremely effective, it needs the  evaluation of the transform at complex arguments. Therefore, the result might  depend on how a certain complex function is defined (the location of the  branch cut.) Mathematica's unique way of handling some of the simple complex functions differs  from the intuitive definition most people are familiar with. The user is  cautioned to check out the relevant comments in Reference [1].

References

The method is described in
1)  Abate, J. and Valkó, P. P.:  Multi-Precision Laplace  Transform Inversion, International Journal of Numerical Methods for  Engineering, (2003), accepted for publication (IJNME A1032.C)

Before using the FT procedure we recommend to try the GWR  algorithm available from  
2) http://library.wolfram.com/infocenter/MathSource/4738/ . 

More information on multi-precision inversion can be found at the  WEB address:  
3) http://pumpjack.tamu.edu/~valko/

For other approaches to the same problem see also:
4) MathSource/Enhancements/Algebraic/0210-968: Numerical Inversion of  Laplace Transform, Arnaud  Mallet
and

5) MathSource/Enhancements/Algebraic/0206-738: Approximate Inversion  of Laplace Transform, Alexander Cheng, Paston Sidauruk