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]}}]