## Quick-and-dirty numerical inversion of the Laplace transform

I present here a set of Mathematica routines implementing a quick-and-dirty method for numerically evaluating the inverse Laplace transform. As it is well-known, the problem is a numerically ill-conditioned one, and thus checking the sensibility of your results is very much important here. The value of the following method, however, compared to more elaborate transform inversion methods, is that if one wants a quick estimate for arguments not too far from the imaginary axis, the approximation derived from this method does not take too much effort to evaluate.

The method, due to Herbert Salzer, is a complex-valued $n$-point Gaussian quadrature, where the nodes are the (complex!) roots of (a special case of) the Bessel polynomial. Due to this, if you want to use this method, the function you wish to transform should be able to evaluate complex arguments.

The method below does not give an error estimate; what you can try for verifying the accuracy of your results is to evaluate for an initial number of points (the default is ten points in the routine below), and then increase the number of points taken, and then compare the results you have. (You might also wish to use a sufficiently high value of the option WorkingPrecision.)

Here are the the Mathematica routines:

SalzerRuleData[n_Integer, prec_:MachinePrecision] := Block[{x}, Transpose[Map[{#, (-1)^(n + 1) # ((2 n - 1)/HypergeometricPFQ[{1 - n, n - 1}, {}, #])^2/n} &, x /. NSolve[HypergeometricPFQ[{-n, n}, {}, x], x, prec]]]]

Options[NInverseLaplaceTransform] = {Points -> 10, WorkingPrecision -> MachinePrecision}; NInverseLaplaceTransform[f_, s_, t_, opts:OptionsPattern[]] := Module[{xa, wt}, {xa, wt} = SalzerRuleData[OptionValue[Points], OptionValue[WorkingPrecision]]; Chop[(wt.Function[s, f][1/xa/t])/t]]

As an example, we try the inverse Laplace transform of the function $\frac1{(1+s)^2}$, which is $t\exp(-t)$:
 f[t_] = NInverseLaplaceTransform[(1 + s)^(-2), s, t];
 With[{t = N[2]}, {t Exp[-t], f[t], t Exp[-t] - f[t]}] {0.270671, 0.270671 + 5.68434*10^-14 I, 2.91759*10^-8 - 5.68434*10^-14 I}

(The tiny imaginary part is an artifact of the evaluation; you can use Chop[] if you wish.)

As another example, here is a comparison of the approximate inverse transform of $\frac1{\sqrt{1+s^2}}$ with the true inverse transform, $J_0(t)$ (where $J_0(t)$ is the Bessel function of the first kind):

f[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, WorkingPrecision -> 25];
 GraphicsRow[{Plot[{f[t], BesselJ[0, t]}, {t, 0, 15}, Frame -> True], Plot[Re[f[t] - BesselJ[0, t]], {t, 0, 15}, Frame -> True, PlotRange -> All]}]

The plot on the left shows that the true function and the approximation are almost indistinguishable within the range plotted. The plot on the right shows the difference between the approximate and the true function.