Wolfram Research

Function Repository Resource:

NInverseLaplaceTransform

Source Notebook

Find the numerical approximation for the inverse Laplace transform

Contributed by: Peter Valko (p-valko@tamu.edu)
professor emeritus of Petroleum Engineering, Texas A&M

ResourceFunction["NInverseLaplaceTransform"][expr,s,t]

gives a numerical approximation to the inverse Laplace transform to expr evaluated at the numerical value t, where expr is a function of the symbol s.

Details and Options

Given a positive t value, the numerical approximation to the inverse Laplace transform of expr can be⁠—in principle⁠—calculated from numerically integrating , where F(s)=expr and γ is an arbitrary positive constant chosen so that the contour of integration lies to the right of all singularities in F(s). However, specialized methods are more reliable. This resource function provides an effective combination of two methods described in the references.
The expr must be numeric whenever symbol s is set to a numeric value.
Approximate numbers (numbers with decimal points) are allowed in the expression.
The t value must be positive, such as 2.9×103 or 4×π2.
ResourceFunction["NInverseLaplaceTransform"] returns $Failed when it cannot determine the numerical result.
ResourceFunction["NInverseLaplaceTransform"] provides the result as a machine number (about 16 significant digits).
There are two options: "Startm" and "Method". (Please notice the quotation marks.)
"Startm" is an integer between 5 and 15. The internal numerical inversion procedures ("FT" and "GWR") are called with 2m terms so "Startm"5 implies that initially 25=32 terms will be used.
"Method" can be "FT", "GWR" or "Automatic". (Please notice the quotation marks.) "Method""Automatic" will use the internal procedures "FT" and "GWR" in an alternating manner, while "FT" or "GWR" will use only the corresponding internal procedure.
The default options are {"Startm"5, "Method""Automatic"}, so ResourceFunction["NInverseLaplaceTransform"][extpr,s,t] is equivalent to ResourceFunction["NInverseLaplaceTransform"][expr,s,t,{"Startm"5, "Method""Automatic"} ]
The curly braces {} can be omitted around the options.

Examples

Basic Examples

The expression of this example has a known symbolic Laplace inverse:

In[1]:=
ResourceFunction["NInverseLaplaceTransform"][
 Exp[1 - Sqrt[s]] s^(-3/4), s, 3.7]
Out[1]=

We can compare the result with the answer from the symbolic evaluation:

In[2]:=
InverseLaplaceTransform[Exp[1 - s^(1/2)] s^(-3/4), s, t]
Out[2]=
In[3]:=
% /. {t -> 3.7}
Out[3]=

This expression cannot be inverted symbolically, only numerically:

In[4]:=
expr2 = ((-1. + 100.*s)*
     Sinh[0.5*Sqrt[s]])/(s*(Sqrt[s]*Cosh[Sqrt[s]] + s*Sinh[Sqrt[s]]));
ResourceFunction["NInverseLaplaceTransform"][expr2, s, 0.2727]
Out[5]=

Scope

The t can be small or large and can contain constants (like π):

In[6]:=
expr3 = Log[p^(1/4)] p^(-1/4);
{ResourceFunction["NInverseLaplaceTransform"][expr3, p, 6.3 10^-3],
 ResourceFunction["NInverseLaplaceTransform"][expr3, p, 6.3 10^3],
 ResourceFunction["NInverseLaplaceTransform"][expr3, p, 2 Pi]}
Out[7]=

Compare with the answer from symbolic evaluation:

In[8]:=
InverseLaplaceTransform[expr3, p, t]
Out[8]=
In[9]:=
% /. t -> {6.3 10^-3, 6.3 10^3, 2. Pi}
Out[9]=

The expression Log[p]p is not a Laplace transform of anything because a Laplace transform must tend to zero at p. Nevertheless, numerical inversion returns a result that makes sense:

In[10]:=
expr4 = Log[p] p;
ResourceFunction["NInverseLaplaceTransform"][expr4, p, 10.]
Out[11]=

One way to look at expr4 is . Note that is the Laplace transform of f[t]=(-EulerGamma-Log[t]), and 0.01 is the second derivative of f[t] at t=10:

In[12]:=
D[InverseLaplaceTransform[Log[p]/p, p, t], {t, 2}]
Out[12]=
In[13]:=
% /. t -> 10.
Out[13]=

In other words, numerical inversion works on a larger class of functions than inversion, but the extension is coherent with the operational rules.

Options

The two options "Startm" and "Method" are introduced here. Considering the following Laplace transform pair:

In[14]:=
expr5 = Log[s]/(1 + s^2);  
f5[t_] = Cos[t]*SinIntegral[t] - Sin[t]*CosIntegral[t];

The inverse f5[t] is periodic-like, but not exactly periodic. At reasonably small t values, there is no problem:

In[15]:=
{f5[20.], ResourceFunction["NInverseLaplaceTransform"][expr5, s, 20.]}
Out[15]=

But already at the moderate value t=100. we get a surprise:

In[16]:=
{f5[100.], ResourceFunction["NInverseLaplaceTransform"][expr5, s, 100.]}
Out[16]=

The result from the NInverseLaplaceTransform function is wrong, and⁠—what is even worse⁠—we do not get an error message. This is what is happening the background: the internal numerical inversion procedures "GWR" (Gaver-Wynn-rho) and "FT" (Fixed-Talbot) provide the same wrong answer initially, so NInverseLaplaceTransform accepts it.

Tackling the problem with option "Startm":

The default value of the option "Startm" is 5. In the "FT" and "GWR" internal inversion procedures 2m terms are used, so the default means that initially the number of terms is set to 25=32. For the given problem, this does not provide enough resolution. We can try a larger "Startm":

In[17]:=
ResourceFunction[
 "NInverseLaplaceTransform"][expr5, s, 100., {"Startm" -> 6}]
Out[17]=

The result is still wrong, but increasing the starting m again provides the correct answer:

In[18]:=
ResourceFunction[
 "NInverseLaplaceTransform"][expr5, s, 100., {"Startm" -> 7}]
Out[18]=

Any larger "Startm" would do the same. Changing the default to "Startm"7⁠—or a larger number⁠—would be, however, a waste of computational resources in most of the cases, so we stick with the default, "Startm"5 .

Tackling the problem with option “Method”:

The other option is "Method". It can be "FT", "GWR" or Automatic. The default is Automatic, meaning that both internal procedures "FT" and "GWR" will be used – in an alternating manner.

If we explicitly set the method to "FT", only that internal procedure will be called. The above problem can be solved with the default "Startm", if we specify "Method""GWR" :

In[19]:=
ResourceFunction[
 "NInverseLaplaceTransform"][expr5, s, 100., {"Method" -> "GWR"}]
Out[19]=

(However, for significantly larger t values an elevated "Startm" is needed, whatever "Method" we specify. See additional remarks in the Author’s Notes.)

Applications

Find the maximum value and location of the inverse Laplace transform of example 2, first graphically and then using FindMaximum:

In[20]:=
expr6 = ((-1. + 100.*s)*
     Sinh[0.5*Sqrt[s]])/(s*(Sqrt[s]*Cosh[Sqrt[s]] + s*Sinh[Sqrt[s]]));
LogLinearPlot[
 ResourceFunction["NInverseLaplaceTransform"][expr6, s, t], {t, 0.001,
   100.}, PlotRange -> All, PlotTheme -> "Detailed"]
Out[21]=
In[22]:=
FindMaximum[
 ResourceFunction["NInverseLaplaceTransform"][expr6, s, t], {t, 0.001,
   100.}]
Out[22]=

In the next application, we need to calculate average temperature. Assume that solving a differential equation we’ve obtained the following analytical result for the temperature as a function of time:

In[23]:=
f7[t_] = (EllipticTheta[3, Pi/2, E^(-(Pi^2*t))] + EllipticTheta[3, Pi, E^(-(Pi^2*t))])/2;

Symbolic integration is unable to calculate the average temperature in the time-interval between zero and one:

In[24]:=
Integrate[f7[t], {t, 0, 1}]
Out[24]=

Now we try numerical integration with various options:

In[25]:=
Do[
 {meth, NIntegrate[f7[t], {t, 0, 1}, Method -> meth]} // Print, {meth, {"GlobalAdaptive", "LocalAdaptive", "DoubleExponential", "MonteCarlo", "AdaptiveMonteCarlo", "QuasiMonteCarlo", "AdaptiveQuasiMonteCarlo"}}]

We see that numerical integration of this function is troublesome, so we recall that the Laplace transform of f7[t] is:

In[26]:=
expr7 = (0.5*(1 + Cosh[Sqrt[s]])*Csch[Sqrt[s]])/Sqrt[s];

Using the operational rule for integration, the average temperature is:

In[27]:=
ResourceFunction["NInverseLaplaceTransform"][expr7/s, s, 1]
Out[27]=

Possible Issues

In this example, the expr is not a Laplace transform, and the inversion attempt fails:

In[28]:=
ResourceFunction["NInverseLaplaceTransform"][Exp[s^2]/s, s, 10]
Out[28]=

For some Laplace transforms, a too large (or too small) t can be detrimental and the inversion fails. Consider the transform of Sin[t]:

In[29]:=
ResourceFunction["NInverseLaplaceTransform"][1/(1 + s^2), s, 1200]
Out[29]=

A larger value of "Startm" will help but it is better to avoid large t-s if the inverse is periodic:

In[30]:=
ResourceFunction["NInverseLaplaceTransform"][1/(1 + s^2), s, 1200, "Startm" -> 6]
Out[30]=

The method fails in the vicinity points where the Laplace inverse, or its derivatives, have a discontinuity; in the following example at t=2, the third derivative of the inverse (the SawtoothWave function) has a discontinuity:

In[31]:=
expr10=(1-s/(Exp[s]-1))/s^5;
Do[{t, NInverseLaplaceTransform[expr10, s, t]} // Print, {t, 1.5, 2.5, 0.1}]

The following inverse has a discontinuity at t=25 and it is zero for t<25 :

In[32]:=
InverseLaplaceTransform[Erf[5*Sqrt[s]], s, t]
Out[32]=

Leaving out the Dirac delta component, the result of numerical inversion should agree with the simple expression:

In[33]:=
f11[t_] = -(5*HeavisideTheta[-25 + t])/(Pi*Sqrt[-25 + t]*t);

Approaching from above with the numerical inversion, we can get quite near to the discontinuity (though with increasing computation time):

In[34]:=
Do[{t, ResourceFunction["NInverseLaplaceTransform"][Erf[5*Sqrt[s]], s,
     t], f11[t], "Current time: " <> DateString["Time"]} // Print, {t,
   25.1, 25., -0.02}]

Numerical inversion works the best if the inverse transform is an analytic function.

Neat Examples

Here the result is a familiar number because the inverse transform is known to be Log[1+Sqrt[t]]:

In[35]:=
expr12 = 0.5 (Pi*Erfi[Sqrt[s]] - ExpIntegralEi[s])/(Exp[s]*s);
ResourceFunction["NInverseLaplaceTransform"][expr12, s, 81.]
Out[36]=
In[37]:=
Log[1 + Sqrt[t]] /. t -> 81.
Out[37]=

In this case the Laplace transform involves the MeijerG function (a rather complicated mathematical function that is nevertheless suitable for numerical evaluation in the Wolfram Language):

In[38]:=
expr13 = MeijerG[{{3/5}, {}}, {{0, 1/5, 2/5, 3/5, 3/5, 4/5}, {}}, s^5/3125]/(4*Sqrt[5]*Pi^2);
ResourceFunction["NInverseLaplaceTransform"][expr13, s, 1.]
Out[39]=

We can check symbolically that the result of the numerical inversion is correct:

In[40]:=
InverseLaplaceTransform[expr13, s, t]
Out[40]=
In[41]:=
% /. t -> 1.
Out[41]=

Here the inverse Laplace transform happens to be the two-parameter Mittag–Leffler function:

In[42]:=
expr14 = (Sqrt[1 + p] - ArcCoth[Sqrt[1 + p]])/(Sqrt[Pi] (1 + p)^(3/2));
ResourceFunction["NInverseLaplaceTransform"][expr14, p, 5.]
Out[43]=
In[44]:=
MittagLefflerE[1, 1/2, -5.]
Out[44]=

Here the Laplace transform expression cannot be cast easily into a “real-only” form (though it is real for positive s values):

In[45]:=
expr15 = 0.5  s^(-1 - I) (s^(2 I) Gamma[-I] + Gamma[I]);

We prepare a log-linear plot:

In[46]:=
LogLinearPlot[
 ResourceFunction["NInverseLaplaceTransform"][expr15, s, t], {t, 0.001, 1000.}, PlotRange -> All, PlotTheme -> "Detailed"]
Out[46]=

The log-linear plot suggests that the inverse is simply Sin[Log[t]]. A quick check proves the assertion:

In[47]:=
SeedRandom[1234]
Do[
 tran = 10^RandomReal[{-2, 8}];
 {tran, ResourceFunction["NInverseLaplaceTransform"][expr15, s, tran] == Sin[Log[tran]]} // Print
 , 5]

Resource History

Source Metadata

See Also

License Information