In this post, I will present an extension of the method previously posted in Quick-and-dirty numerical inversion of the Laplace transform.

The method, due to Piessens, is effectively a Gauss-Kronrod extension of the complex Gaussian quadrature method by Salzer. Recall that the Gauss-Kronrod idea is to use two different quadrature abscissa/weight pairs, where one’s abscissas is a subset of the other, to be able to estimate the accuracy of the higher-order method.

Here is the implementation of the extended version of `NInverseLaplaceTransform[]`

:

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

piessensQ[n_Integer, x_] := Module[{cof},
cof = LinearSolve[SparseArray[{{j_, k_} /; j <= k :> 1/((k - j)!
Binomial[2 n + k - j, n + k - j])},
{n + 1, n + 1}],
-Table[1/(j! Binomial[2 n + j, n + j]),
{j, n + 1, 1, -1}]];
Expand[FromDigits[Prepend[Reverse[cof], 1], x]]]

SalzerPiessensRuleData[n_Integer, prec_: MachinePrecision] :=
Module[{c, pk, pl, sg, sk, sl},
sl = \[FormalX] /.
NSolve[HypergeometricPFQ[{-n, n}, {}, \[FormalX]], \[FormalX], prec];
pl = \[FormalX] /.
NSolve[piessensQ[n, \[FormalX]], \[FormalX], prec];
sk = HypergeometricPFQ[{1 - n, n - 1}, {}, sl];
sg = (-1)^(n + 1) sl ((2 n - 1)/sk)^2/n; c = n! Binomial[2 n, n];
sk = sg - (2 n - 1) sl/(c n sk piessensQ[n, sl]);
pk = 1/(c pl HypergeometricPFQ[{-n, n}, {}, pl]
(D[piessensQ[n, \[FormalX]], \[FormalX]] /.
\[FormalX] -> pl));
{Join[sl, pl], Join[sk, pk]}]

Options[NInverseLaplaceTransform] = {Method -> Automatic,
"Points" -> 6, WorkingPrecision -> MachinePrecision};
NInverseLaplaceTransform[f_, s_, t_, opts : OptionsPattern[]] :=
Module[{met, xa, wt},
met = Switch[OptionValue[Method],
"SalzerPiessens" | Automatic, SalzerPiessensRuleData,
"Salzer", SalzerRuleData, _, Return[$Failed]];
{xa, wt} = met[OptionValue["Points"],
OptionValue[WorkingPrecision]];
Chop[(wt.Function[s, f][1/xa/t])/t]]

For error estimation purposes, one could then generate the approximant corresponding to both Salzer and Salzer-Piessens for the same setting of `"Points"`

and look at their difference.

Using the Bessel function example from the previous post,

f1[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t,
"Points" -> 9, WorkingPrecision -> 30];
f2[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t,
Method -> "Salzer", "Points" -> 9, WorkingPrecision -> 30];

GraphicsGrid[{{Plot[{BesselJ[0, t], Re[f1[t]]}, {t, 0, 19},
Frame -> True, ImageSize -> Medium,
PlotLegends -> {{TraditionalForm[BesselJ[0, t]],
"Salzer-Piessens"}},
PlotStyle -> {AbsoluteThickness[4], Automatic}],
SpanFromLeft},
{Plot[Re[f1[t] - f2[t]], {t, 0, 19}, Frame -> True,
PlotLabel -> Style["Difference between Salzer-Piessens
and Salzer methods", Small],
PlotRange -> All],
Plot[Re[f1[t] - BesselJ[0, t]], {t, 0, 19}, Frame -> True,
PlotLabel -> Style["Difference between Salzer-Piessens and
true transform", Small],
PlotRange -> All]}}]