Numerical inverse Laplace-Hankel transform

Question

When trying to reproduce the result of this paper about numerical solution of Lamb's problem, I encountered the following double integral (to be more precise, the 0-order inverse Hankel-Laplace transform):

$$\frac{1}{2 \pi i} \int _{\gamma -i \infty }^{\gamma +i \infty }e^{pt}\int _0^{\infty }f(p,\xi )\xi J_0(\xi r)d\xi dp$$

where

r = 1;
f[p_, ξ_] = -(5 p Sqrt[(5 p^2)/6 + ξ^2] )/(
  4 (-4 ξ^2 Sqrt[(5 p^2)/6 + ξ^2] Sqrt[(5 p^2)/2 + ξ^2] + ((5 p^2)/2 + 2 ξ^2)^2));

f[p, ξ] is actually derived from equation (16b) in the paper.

Note: There's a typo in the paper, the relationship between $u_z^{norm}$ and $u_z$ should be $u_z^{norm}=\pi \mu (r^2+z^2)^{1/2}u_z/F$

To calculate the transform, the paper, published in 1999, utilized a quite complicated procedure. Given that the numerical calculation capability of Mathematica has improved a lot since then, I wonder if the transform can be calculated in a more straightforward way nowadays?

For comparison, the following is the desired graph of transformed result:

analyticSol = 
  "1:eJzNVktOwzAQTfMpLUKINTtW7LgEQrDoIqp6AQOuOlJKSpzSHo8jcB4EBDtRPrYnrpuGz6Kjun6eeX5+munFfTydu47jsDEPIdAHugFGi598HibAUmXliVXAwwyWlBWbIx6mJIX4iUQw4As446HYE9Aw3tBk7pd5wmjNYCgh8mTwlmWZcoqJeilZg4vgV3qVCu8j+FdD/pqPdBt+ImvULukrMuQwJXPjlK6PqzNXyMr1lAJFUV9PIsvlSSk+kRQID/CseWNEjxHd3xF17Ku0+WAPdrXIJkk6XX/3K+zMVZnKTF3WOUCoeDoVy+vs4dDDLNnXy3W/ak32w6rd7Pdy3ZggPvxpsSs5J5TZ8q+aYVVmIEmRDwgfecUjHq7j5SqiWyStuH3jMtKeGEwwQrwo2UGrAltxbtPmqibyUiDP29zeRF4J5NggrGdT70Qgg0PruWU9bPRYNH5sYP1C57Mbp7oUuQ2G+jkN6MlAQ8pKwQMn0/8Vu84VyLlae3tHt7d05T7+YvzZ6BPf2CkPd5S8AINHOlvQlBj7vlGNPsdxcecvbqvSA8We6C23CSUpTW6e1ySS+jTwj/MNm2tw+Q==" // 
   Uncompress;
Plot[analyticSol, {tau, 0, 1.2}, PlotRange -> {-0.8, 1.6}, Exclusions -> None, 
 AxesLabel -> {Style["τ", 20], None}, AxesOrigin -> {0, -0.1}]

analyticSol i.e. analytic solution of the transform can be found in p743 of A.Cemal Eringen, Erdogan S. Suhubi, Elastodynamics. Linear theory, Academic Press Inc, 1975.

Notice here tau = c2 t/r; c2 = Sqrt[2/5]; r = 1. analyticSol == -3/8 ≈ 0.375 when tau > Sqrt[3 + Sqrt[3]]/2.

I myself have made several trial, but none of them works well.

Trial 1, direct integration:

With[{t = 2}, 
 AbsoluteTiming[
  1/(2 Pi I) NIntegrate[ξ BesselJ[0, ξ] f[p, ξ] Exp[p t], {p, -I ∞, I ∞}, 
   {ξ, 0, ∞}, WorkingPrecision -> 16, MaxRecursion -> 40, 
   Method -> "ExtrapolatingOscillatory", Exclusions -> Denominator@f[p, ξ] == 0]]]
(* Returned unevaluated after some warning *)

Trial 2, GWR:

hankel[p_, prec_: 16] := 
 NIntegrate[ξ BesselJ[0, ξ] f[p, ξ], {ξ, 0, ∞}, 
  WorkingPrecision -> prec, MaxRecursion -> 40, Method -> "ExtrapolatingOscillatory", 
  Exclusions -> Denominator@f[p, ξ] == 0]

minushankel = Function[{p}, #] &@(DownValues@hankel /. {f@a__ :> -f@a, 
       NIntegrate -> nIntegrate})[[1, -1]] /. nIntegrate -> NIntegrate

-GWR[minushankel[#, 64] &, 2] // AbsoluteTiming
(* {165.339792, -0.31950524676} <- large error *)    

Here I introduced a minushankel function because of another issue, I've started a separate post for it.

Trial 3, FT:

FT[hankel, 2] // AbsoluteTiming
(* {300.864984, -0.4128181} <- Large error *)

Trial 4, customized Crump method (coding according to the paper):

crumpsp[f_, t_, prec_: MachinePrecision] := 
 Module[{c2 = Sqrt[2/5], r = 1, z = 0, taumax = 2, tmax, T, w, q = 0, ee = 10^-8, 
   lstk = SetPrecision[Range@50, prec], mid}, tmax = taumax Sqrt[r^2 + z^2]/c2; 
  T = 8/10 tmax; w = q - Log[ee]/(2 T); 
  mid =(*Parallel*)
   Map[f, (w + (I lstk π)/T)(*,DistributedContexts->Automatic*)]; 
  Exp[w t]/T (1/2 f@w + Total[Re@mid Cos[lstk Pi t/T] - Im@mid Sin[lstk Pi t/T]])]

crumpsp[hankel[#, 16] &, 2] // AbsoluteTiming
(* {304.024717, -0.362223}  <- Not that small error *)

BTW, it's frustrating yet worth to mention, the author of the paper claimed that the computation time on a PC Pentium/120 MHz for obtaining each graph was typically 5 min.

How to obtain an accurate enough result of the transform (in a not that slow way with built-in functions/existing packages, if possible)?

Answer 1

The following is not particularly fast and could be more accurate but does make progress toward the goals set in the Question. To begin, consider the analytical properties of f, as defined in the Question. By observation, it has branch points at

{I Sqrt[6/5] ξ, I Sqrt[2/5] ξ}

and their conjugates. Poles are obtained by

p /. Simplify[Solve[Denominator[f[p, ξ]] == 0, p], ξ > 0] // Flatten
(* {0, -2 I Sqrt[1/15 (3 - Sqrt[3])] ξ, 2 I Sqrt[1/15 (3 - Sqrt[3])] ξ} *)

with corresponding residues,

Simplify[Residue[f[p, ξ], {p, #}] & /@ %, ξ > 0]
(* {-(3/(8 ξ)), Sqrt[1 + Sqrt[3]]/(8 (Sqrt[33 - 19 Sqrt[3]] - 2 Sqrt[3 - Sqrt[3]] + 
    Sqrt[9 + 5 Sqrt[3]]) ξ), ...} *)

Remark

If you're in v9.0.1, the Residues of $\pm 2 i \sqrt{\frac{1}{15} \left(3-\sqrt{3}\right)} \xi$ won't be calculated correctly due to a bug, please use

Simplify[(SeriesCoefficient[Exp[p t] f[p, ξ], {p, #, -1}] &) /@ %, ξ > 0]

instead to calculate the residues.

(The third Residue is the same as the second.) The poles and branch cuts can be seen in the plots, here for ξ == 1,

Row[Plot3D[#, {pr, -1, 1}, {pi, -1, 1}, ImageSize -> Medium, 
PlotRange -> {-2, 2}, AxesLabel -> {pr, pi, f}] & /@ ReIm[f[pr + I pi, 1]]]

all lying along the imaginary axis, and the branch cuts running from the corresponding branch points to I Infinity and its conjugate.

From these results it is evident that the Bromwich contour of the first expression in the Question must lie to the right of the imaginary axis and should be somewhat to the right to minimize integration issues. The integral over p as a function of ξ for t == 2 is

Plot[(ξ NIntegrate[f[1 + I pi, ξ] Exp[2 (1 + I pi)], {pi, -Infinity, Infinity}] // Chop)
    /(2 Pi), {ξ, 0, 50}, PlotRange -> All]

This computation, which takes a few minutes on my PC, seems quite accurate, as can be seen by varying the real part of p for the contour integration. The integration over ξ is neither fast nor accurate, unfortunately.

pinv[ξ_?NumericQ] := NIntegrate[f[1 + I pi, ξ] Exp[2 (1 + I pi)]/(2 Pi), 
    {pi, -Infinity, Infinity}] // Chop
NIntegrate[ξ BesselJ[0, ξ] pinv[ξ], {ξ, 0, Infinity}, 
    MaxRecursion -> 40, Method -> "ExtrapolatingOscillatory"]
(* SequenceLimit::seqlim: The general form of the sequence could not be determined, 
   and the result may be incorrect. *)
(* -0.423105 *)

rather than the desired -3/8.

However, further progress can be made be rewriting f as

f1[p_,ξ_] = -5 p Sqrt[Sqrt[5/6] p + I ξ] Sqrt[Sqrt[5/6] p - I ξ]/
    (4 (-4 ξ^2 Sqrt[Sqrt[5/6] p + I ξ] Sqrt[Sqrt[5/6] p - I ξ] Sqrt[Sqrt[5/2] p + I ξ] 
    Sqrt[Sqrt[5/2] p - I ξ] + ((5 p^2)/2 + 2 ξ^2)^2))

which causes the branch cuts to extend along lines of constant Im[p] from - Infinity to the branch points. The Bromwich contour then can be distorted to yield a solution consisting of the three Residues plus integrals of the discontinuities along the four branch cuts. The Residues alone give a qualitatively reasonable rendition of the plot in the Question.

Plot[Integrate[BesselJ[0, ξ] Cos[ξ t 2 Sqrt[1/15 (3 - Sqrt[3])]], {ξ, 0, Infinity}] 
    2 Sqrt[1 + Sqrt[3]]/(8 (Sqrt[33 - 19 Sqrt[3]] - 2 Sqrt[3 - Sqrt[3]] + 
    Sqrt[9 + 5 Sqrt[3]]) ) - 3/8, {t, 0, 2}, PlotRange -> {-.5, 1.5}, AxesLabel -> {t, F}]

Although the value of the integral is -3/8, as desired, for large t, the curve is not accurate at small t and, of course, does not include the integrals along the branch cuts. These are given for t == 2 by

inv65[ξ_?NumericQ] := ξ Im[NIntegrate[(f1[ξ Sqrt[6/5] I + pr - .00001 I, ξ] - 
     f1[ξ Sqrt[6/5] I + pr + .00001 I, ξ]) Exp[2 pr], {pr, -Infinity, 0}] 
     Exp[2 ξ Sqrt[6/5] I]]/Pi
inv25[ξ_?NumericQ] := ξ Im[NIntegrate[(f1[ξ Sqrt[2/5] I + pr - .00001 I, ξ] - 
     f1[ξ Sqrt[2/5] I + pr + .00001 I, ξ]) Exp[2 pr], {pr, -Infinity, 0}] 
     Exp[2 ξ Sqrt[2/5] I]]/Pi
Plot[{inv65[ξ], inv25[ξ]}, {ξ, 0, 50}, PlotRange -> {-.04, .04}]

These branch cut integrals offer two advantages over the Bromwich integration discussed above. First, the integrals themselves converge more rapidly due to the factors Exp[2 pr]. Second, the results decrease rapidly with ξ. Hence, the contribution to the t == 2 total result,

NIntegrate[BesselJ[0, ξ] (inv65[ξ] + inv25[ξ]), {ξ, 0, Infinity}, 
    MaxRecursion -> 40, Method -> "ExtrapolatingOscillatory"]
(* 0.000191503 *)

is acceptably close to the desired value of 0 despite error messages. Presumably, a smaller value could be obtained using larger values of WorkingPrecision, although the already slow computation would be far slower still. The curve in the Question can be reproduced reasonably well by

inv[ξ_?NumericQ, t_?NumericQ] := 
    ξ Im[NIntegrate[(f1[ξ Sqrt[6/5] I + pr - .00001 I, ξ] - 
      f1[ξ Sqrt[6/5] I + pr + .00001 I, ξ]) Exp[t pr], 
      {pr, -Infinity, 0}] Exp[t ξ Sqrt[6/5] I]]/Pi + 
    ξ Im[NIntegrate[(f1[ξ Sqrt[2/5] I + pr - .00001 I, ξ] - 
      f1[ξ Sqrt[2/5] I + pr + .00001 I, ξ]) Exp[t pr], 
      {pr, -Infinity, 0}] Exp[t ξ Sqrt[2/5] I]]/Pi
tinv = Table[NIntegrate[BesselJ[0, ξ] inv[ξ, t], {ξ, 0, Infinity}, 
    MaxRecursion -> 40, Method -> "ExtrapolatingOscillatory"] + 
    Integrate[BesselJ[0, ξ] Cos[ξ t 2 Sqrt[1/15 (3 - Sqrt[3])]], {ξ, 0, Infinity}] 
    2 Sqrt[1 + Sqrt[3]]/(8 (Sqrt[33 - 19 Sqrt[3]] - 2 Sqrt[3 - Sqrt[3]] + 
    Sqrt[9 + 5 Sqrt[3]]) ) - 3/8, {t, .01, 2, .01}];
ListLinePlot[tinv, DataRange -> {0.01 Sqrt[2/5], 2 Sqrt[2/5]}, 
    PlotRange -> {{0, 1.2}, {-0.8, 1.6}}, AxesLabel -> {Style["τ", 20], None}, 
    AxesOrigin -> {0, -.1}]

This curve required a few hours to produce, although significant savings could have been obtained with ParallelTable and some simplification of the integrands. That the computations in the reference given in the Question are so much faster reinforces the observation that superior algorithms typically trump superior coding.

Addendum

In response to a comment below, this additional material elaborates on the inverse Laplace transform of f1 and compares it with the inverse Laplace transform of f. As already discussed, f1 exhibits four branch cuts and three poles, shown in red in the plot immediately below for ξ == 1.

The Bromwich integration contour described in the Question is a vertical line in the p plane, lying to the right of all branch cuts and poles, shown in the plot as a dashed line. Although the integration indicated by the first expression in the Question can be performed along this contour, the discussion in the latter part of this Answer indicated that it can be useful to transform the contour by shifting it far to the left, where it becomes a set of contours around the branch cuts and poles, as illustrated in the plot. (To minimize clutter, the contours are shown for only one each of the cuts and poles.) The portions of the shifted contour that are at p == - Infinity are ignored, because the integral vanishes there due to the integrand factor Exp[p t]. Integrals along the branch cuts in effect integrate the discontinuity of f1 across those cuts, and the integrals around the poles are, of course, their Residues.

Just as the contour segments shown in the plot can be separated into components arising from the poles and from the branch cuts, the same can be done for the result of the original Bromwich contour integration in the early part of this Answer. The next plot compares (1) the sum of the branch cut integrals from the fourth plot above, and (2) the result in the second plot above after subtracting the pole contributions.

Plot[{ξ Im[NIntegrate[(f1[ξ Sqrt[6/5] I + pr - .00001 I,ξ] - 
    f1[ξ Sqrt[6/5] I + pr + .00001 I, ξ]) Exp[2 pr], {pr, -Infinity, 0}] 
    Exp[2 ξ Sqrt[6/5] I]]/Pi 
    + ξ Im[NIntegrate[(f1[ξ Sqrt[2/5] I + pr - .00001 I, ξ] - 
    f1[ξ Sqrt[2/5] I + pr + .00001 I, ξ]) Exp[2 pr], {pr, -Infinity, 0}] 
    Exp[2 ξ Sqrt[2/5] I]]/Pi, 
    (ξ NIntegrate[f[1 + I pi, ξ] Exp[2 (1 + I pi)], {pi, -Infinity, Infinity}] // Chop)
    /(2 Pi) - 
    Cos[ξ 4 Sqrt[1/15 (3 - Sqrt[3])]] 2 Sqrt[1 + Sqrt[3]]/(8 (Sqrt[33 - 19 Sqrt[3]]
    - 2 Sqrt[3 - Sqrt[3]] + Sqrt[9 + 5 Sqrt[3]]) ) + 3/8}, {ξ, 0, 50}, 
    PlotRange -> {-.08, .08}]

Not surprisingly, the two curves are virtually identical. As explained earlier, this decomposition is desirable, because the ξ integration can be performed symbolically for the pole contributions, resulting in the third plot. Indeed, the seemingly simple appearance of the branch cut contributions in the fourth plot suggest that they too could be represented by symbolic expressions, which might lead to further symbolic ξ integrations.