11
$\begingroup$

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]

enter image description here

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

$\endgroup$
6
  • $\begingroup$ I guess it has to do with the Bromwich integral. $\endgroup$
    – Dimitris
    Apr 29, 2015 at 12:23
  • $\begingroup$ As it stands, Mathematica does not know how to do contour integrations of this sort. You can, however, supply an appropriate Hankel contour for computations like these. $\endgroup$ May 4, 2015 at 10:39
  • $\begingroup$ Long time no see, @J.M. :) How about giving an answer? I'm not that familiar with Hankel contour. $\endgroup$
    – xzczd
    May 4, 2015 at 13:07
  • $\begingroup$ Not having Mathematica on hand, I won't be able to write an answer for a while. However, the last time I looked at this, I managed to use Hankel contours for numerically evaluating inverse Laplace transforms; whether they'll also work symbolically, I am not as certain. $\endgroup$ May 4, 2015 at 13:36
  • 1
    $\begingroup$ To recap: you require a symbolic solution resulting from Integrate[], correct? $\endgroup$ Apr 9, 2016 at 8:53

2 Answers 2

7
+100
$\begingroup$

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]

enter image description here

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

enter image description here

Hence, int can be simplified in this limit to

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

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ Though not a complete answer, thanks for your effort :) $\endgroup$
    – xzczd
    Apr 16, 2016 at 16:42
  • $\begingroup$ @xzczd Thank you. I believe that the integration contour would need to be changed to obtain most of the other answers. In a few cases, the branch cuts might also need to be moved. $\endgroup$
    – bbgodfrey
    Apr 16, 2016 at 16:47
1
$\begingroup$

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

enter image description here

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]

enter image description here

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]

enter image description here

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]

enter image description here

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}

enter image description here

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. *)
    
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.