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

February 1, 2017

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

invlap

Advertisements