Compute inverse Laplace transform with Integrate

Question

Since inverse Laplace transform is just an integral, while we already have InverseLaplaceTransform built in, can we compute it with Integrate? The direct implementation of the definition of the inverse Laplace transform seems not to work:

In the following piece of code, il is the direct implementation of the Bromwich integral, elements of ex are all taken from the document of InverseLaplaceTransform.

il[f_, s_, t_] := Module[{r}, 1/(2 π I) Integrate[f Exp[s t], {s, r - I ∞, r + I ∞}]]

ex = {1/(1 + s), Log[s]^2/s, Erf[5 Sqrt[s]], s/(s^2 + 1), 1/(s^2 - 4), 
      1/(1 + s + 5 s^2), 1/(1 + Sqrt[s]), E^(-s)/s, Log[s]/s^4,
      ArcTan[2/s], Gamma[s, 0, 1], PolyGamma[2, s], 
      Exp[s + 1] ExpIntegralEi[-s], 1/(s^2 + s + 1) 1/s, 1/(s^2 + 1), 
      MeijerG[{{1/3}, {}}, {{0, 1/3, 1/3, 2/3}, {}}, s^3/27]/(2 Sqrt[3] π)};

Grid[ParallelMap[{InverseLaplaceTransform[#, s, t], il[#, s, t]}&, ex], Frame -> True]

As shown above, il fails in all cases, but is there really nothing more that Integrate can do?

Answer 1

Direct Integration

Possible issues in performing the integrations include choice of Assumptions, branch cuts in the integrands, and how limits are taken. Addressing the first of these gives five solutions.

il[f_, s_, t_] := Module[{r}, 1/(2 π I) Integrate[f Exp[s t], {s, r - I ∞, r + I ∞}, 
    Assumptions -> r > 2 && t > 0]]

Grid[ParallelMap[{InverseLaplaceTransform[#, s, t], il[#, s, t]} &, 
    Extract[ex, {{1}, {2}, {5}, {6}, {10}}], Method -> "FinestGrained"], Frame -> True]

The last result does not agree with the corresponding result from InverseLaplaceTransform, perhaps due to the location of the branch cuts of ArcTan[2/s].

Indefinite Integrals

Although Integrate returns unevaluated for all definite integrals, il, except those in the preceding table, Integrate returns indefinite integrals for many of the expressions in ex. For instance,

int = Grid[ParallelMap[{InverseLaplaceTransform[#, s, t], 
    Integrate[# Exp[s t]/(2 π I), s, Assumptions -> t > 0]} &, 
    Extract[ex, {{1}, {4}, {5}, {6}, {14}, {15}}], Frame -> True];

returns results in which s appears only in terms of the form, ExpIntegralEi[s t + c[t]], where c[t] is a polynomial in t. For instance, ex[[4]] integrates to

(* -((I E^(I t) ExpIntegralEi[(-I + s) t])/(4 π)) 
   - (I E^(-I t) ExpIntegralEi[(I + s) t])/(4 π) *)

For st a large positive or negative imaginary number, c[t] can be ignored, and

ExpIntegralEi[s t] - ExpIntegralEi[- s t] == 2 Pi I

in this limit, as can be seen from Eq (7) the Wolfram MathWorld description of this function or simply by plotting the left side of this last expression.

Plot[Evaluate[ReIm[ExpIntegralEi[I/x] - ExpIntegralEi[-I/x]]], {x, 0, 0.2}]

Hence, int can be simplified in this limit to

(int /. ExpIntegralEi[_] -> 2 Pi I) // Expand // ExpToTrig // FullSimplify

Inverse Laplace Transforms obtained by integration have, of course, been obtained in the first table for elements {1, 5, 6} of ex, and those transforms agree with those in the second table. Incidentally, it may not be obvious that the two results for ex[[15]] in the second table are equal, They are.

FullSimplify[1 - (E^(-t/2) (Sqrt[3] Cos[(Sqrt[3] t)/2] + Sin[(Sqrt[3] t)/2]))/Sqrt[3]]
    - (E^(-(-1)^(1/3) t) (-1 - I Sqrt[3] + 2 (1 + (-1)^(1/3)) E^((-1)^(1/3) t) 
        - 2 E^(I Sqrt[3] t)))/(2 (1 + (-1)^(1/3)))
(* 0 *)

Thus, we now have results for {1, 2, 4, 5, 6, 10, 14, 15} of ex, although the result for ex[[10]] does not agree with the InverseLaplaceTransform result, as explained above.

At first glance, the indefinite integral of ex[[13]] should be evaluated in a similar fashion.

(* -((I E^(1 + s (1 + t)) ExpIntegralEi[-s])/(2 π (1 + t))) 
   + (I E ExpIntegralEi[s t])/(2 π (1 + t)) *)

The second term in this expression reduces in the limit of s t a large imaginary number to -(E/(1 + t)), the result returned InverseLaplaceTransform. However, the first term does not vanish in this limit but instead oscillates at an every faster rate. Probably, this is so, because ex[[13]] has a branch cut along the entire positive real axis, so the path integral in the question cannot be carried out legitimately.

Indefinite integrals also can be computed for ex[[3]] in terms of Erf and Erfi functions, and for ex[[9]] (and also ex[[2]]) in terms of HypergeometricPFQ functions. I have not attempted to obtain limits of these expressions for s t a large imaginary number. Integrate again returns unevaluated for elements {7, 11, 12, 16} of ex.

Answer 2

Still, this isn't a complete answer, but a stepfoward.

There're at least 3 issues:

1. The capability of Integrate seems to be limited, in other words, InverseLaplaceTransform seems to be targetedly optimized.

2. As noticed by bbgodfrey, the assumption t > 0 is important to make Integrate work.

3. The r in my il function actually doesn't have any effect.

3rd item is worth elaborating. We know (OK, perhaps many of us don't know) Integrate doesn't have Hold* attribute:

Attributes@Integrate
(* {Protected, ReadProtected} *)

This means the arguments of Integrate will be evaluated before going into Integrate. Then guess what will happen on the r - I ∞ and r + I ∞?:

{s, r - I ∞, r + I ∞}

Aha, the r simply disappears! This feature isn't surprising, but easy to overlook.

To make the r play a role in integration, we modify il to the following:

ilfixed[f_, s_, t_, r_ : 10] := 
 1/(2 π I) Integrate[
   f Exp[s t] /. s -> s + r, {s, -I ∞, +I ∞}, 
   Assumptions -> t > 0]

Since v13.3 ContourIntegrate is introduced, we'll test one more approach:

ilnew[f_, s_, t_, r_ : 10] := 
 1/(2 π I) ContourIntegrate[f Exp[s t], 
   s ∈ InfiniteLine[{r, 0}, {0, 1}], Assumptions -> t > 0]

Given that $\gamma$ (r in code) is an arbitrary positive constant chosen so that the contour of integration lies to the right of all singularities for inverse Laplace transform, let's first check the singularities of ex. Luckily FunctionSingularities has been introduced in v12.2:

Grid[MapIndexed[{#2[[1]], N@Reduce[FunctionSingularities[#, s], s]} &, ex], 
 Frame -> All]

11th result looks a bit suspicious, let's check the ComplexPlot for safety:

pic = 
   ParallelTable[
    ComplexPlot[e, {s, -10 - 10 I, 10 + 10 I}, PlotLabel -> e], {e, 
     ex}]; // AbsoluteTiming
(* {545.658, Null} on 2-core laptop. *)

Grid[MapIndexed[{#2[[1]], #} &, pic]\[Transpose] // Map[Partition[#, 8] &] // 
   Transpose // Flatten[#, 1] &, Frame -> All]

OK, FindSingularities indeed mis-calculates ex[[11]]. From the analysis above, we know that, except for the special ex[[16]], a naive positive $\gamma$ (r in code) should work in principle. Let's try it:

freeze = # -> Inactive@# & /@ {InverseLaplaceTransform, Integrate, 
    ContourIntegrate};

lst = 
   ParallelTable[{InverseLaplaceTransform[e, s, t], ilfixed[e, s, t, 0], 
      ilfixed[e, s, t], ilnew[e, s, t]} /. freeze, {e, ex}]; // AbsoluteTiming
(* {1245.389, Null} on 2-core laptop. *)

For comparison I've also calculated the integration on $\gamma=0$. The results are as follows:

typeset = Function[expr, 
   If[! FreeQ[expr, Inactive], "×", expr] /. _Piecewise :> 
     "Piecewise[…]", Listable];

Grid[Prepend[
  MapIndexed[Prepend[#, #2[[1]]] &, typeset@lst] // FullSimplify[#, t > 0] &, {"", 
   "InverseLaplaceTransform", "γ = 0, Integrate", "γ = 10, Integrate", 
   "γ = 10, ContourIntegrate"}], Frame -> All]

As we can see, by moving the contour, Integrate finds the correct answer for ex[[10]], and finds solutions for ex[[{4, 8, 14, 15}]], but it's not immediately clear if results for ex[[{8, 14}]] are correct or not. Let's plot them as a quick test:

Prepend[MapThread[
    Plot[#, {t, 0, 2}, WorkingPrecision -> 16, 
      PlotLabel -> #2] &, {lst[[#, {1, 3, 4}]], {InverseLaplaceTransform, Integrate, 
      ContourIntegrate}}], #] & /@ {8, 14}

As shown above, Integrate doesn't handle ex[[8]] properly, but the inverse transform of ex[[14]] is correct.

To sum up, by setting $\gamma$ (r in code) in correct way, we now manage to calculate the inverse transform with Integrate for ex[[{1, 2, 4, 5, 6, 10, 14, 15}]]. If we turn to ContourIntegrate, we can further calculate inverse transform of ex[[{7, 8, 9}]].

Remark

1. It's a bit funny that for ex[[3]], even InverseLaplaceTransform itself fails nowadays. By checking the legacy document, we can find this example has been removed in v12.2:

   v12.1 document

   v12.2 document

2. The assumption t > 0 isn't actually necessary for ContourIntegrate.

3. Though owning the attribute HoldAll, NIntegrate suffers from the I ∞ issue, too. See this answer for more info.

4. The following should be a correct setup for calculating Laplace inversion of ex[[16]], but sadly it doesn't work at least for now:

    ilnew2[f_, s_, t_, line_] := 
     1/(2 π I) ContourIntegrate[f Exp[s t], s ∈ line, Assumptions -> t > 0]
   
    ilnew2[ex[[16]], s, t, HalfLine[{{0, 0}, {1, 1}} + Threaded@{1, 0}]] - 
      ilnew2[ex[[16]], s, t, 
       HalfLine[{{0, 0}, {1, -1}} + Threaded@{1, 0}]] // AbsoluteTiming
    (* ContourIntegrate returns unevaluated. *)