Example notebook for the package  'NumericalLaplaceInversion'

Authors

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

Summary

The package provides only one function:GWR.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 other words, approximate numbers (with decimal point) or Mathematica functions starting with the letter 'N' are not allowed.

Needs

The "NumericalLaplaceInversion.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 NumericalLaplaceInversion.m file) into the same folder where the calling notebook (in this case NumericalLaplaceInversionExample.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:

<< NumericalLaplaceInversion.m

The only function available is GWR.

Help

Help is available:

?GWR

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 = GWR[fun, t0];

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

Note 1: There is no garantee, that the result is as accurate as the precision returned by the Kernel indicates. 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 = 4;  ft = 0; 
    Print["    M ", "  GWR"];
While[
 ftold = ft; i++;  M = 2^i; 
        ft = GWR[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)] ];
GWR[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[
      GWR[1/(1 + #^2) &, 50, 128], {29, 29}]
N[Sin[50], 29]

References

The method is described in:
1) Valkó, P.P.and Abate J.: Comparison of Sequence Accelerators for the Gaver Method of Numerical Laplace Transform Inversion, Computers & Mathematics with Applications, (2002), accepted for publication (CAM  5307)

More information on multi-precision inversion can be found in:
2) Valkó, P.P.and Vajda, S : Inversion of noise-free Laplace transforms: Towards a standardized set of test problems,  Inverse Problems in Engineering, (2002) vol .10.No.5, pp 467-483.

or 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