TransWikia.com

How can I invert a Laplace transform numerically?

Mathematica Asked on January 9, 2021

I have a very complicated expression, which I want to transform using the inverse Laplace transform. The built-in function InverseLaplaceTransform doesn’t work. So, I would like to do it numerically. What are the best ways/packages to do that?

Function

fnm[s_, n_] := a1[s] y[s]^n + a2[s] y[s]^-n;

where a1, a2, y and c2 are complicated functions of s. For example,

a1[s_] :=
  (kon (y[s]^m + y[s]^(2 L - m + 1))) /
    (c2[s] (y[s]^m + y[s]^(1 - m)) (y[s]^m  + y[s]^(2 L - m + 1)) + 
      u (1 - y[s]^2) (1 - y[s]^(2 L)));

y[s_] := 
  (s + 2 u + kon - (kon koff)/(s + koff) - 
    Sqrt[(s + 2 u + kon - (kon koff)/(s + koff))^2 - 4 u^2])/(2 u);

c2[s_] := (kon koff)/(s + koff);

a2[s_] := 
  (kon y[s] (y[s]^m + y[s]^(2 L - m + 1))) /
    (c2[s] (y[s]^m + y[s]^(1 - m)) (y[s]^m  + y[s]^(2 L - m + 1)) + 
      u (1 - y[s]^2) (1 - y[s]^(2 L)));

One Answer

As of v12.2, InverseLaplaceTransform supports numeric Laplace inversion. In addition, there exist at least 6 Mathematica packages for numeric inverse Laplace transform. They're:

Currently I only have experience with the first 2 packages so I'll solve OP's problem only with GWR and FT in this answer. (I might add illustrations for other packages later. )

Since OP doesn't show the possible value of parameters, they'll be casually chosen as:

test[s_] = Block[{kon = 1, koff = 2, u = 3, L = 4, m = 5, n = 6}, fnm[s, n]];

InverseLaplaceTransform solution

As mentioned above, this feature is added in v12.2:

InverseLaplaceTransform[test[s], s, 10^-3.]
(* -0.983223 *)

With[{eps = 10^-3}, 
     Plot[InverseLaplaceTransform[test[s], s, N@t],
          {t, eps, 3 + eps}, MaxRecursion -> 0, PlotRange -> All]] // AbsoluteTiming 

enter image description here

Notice 3rd argument of InverseLaplaceTransform cannot be an exact number, or the numeric methods won't be invoked.

GWR solution

The usage of GWR function has been explained clearly in NumericalLaplaceInversionExample.nb so I'd like not to repeat it here.

With[{eps = 10^-3.}, 
   rstGWR = Table[{t, GWR[test, t]}, {t, eps, 3 + eps, 1/25}]]; // AbsoluteTiming
(* {2.546700, Null} *)

ListLinePlot[rstGWR, PlotRange -> All]

Mathematica graphics

FT solution

The usage of FT function has been documented in FixedTalbotNumericalLaplaceInversionExample.nb so, once again, I'd like not to repeat.

With[{eps = 10^-3.}, 
   rstFT = Table[{t, FT[test, t]}, {t, eps, 3 + eps, 1/25}]]; // AbsoluteTiming
(* {1.765505, Null} *)

ListLinePlot[rstFT, PlotRange -> All]

The result just looks the same as that of GWR so I'll omit it here. The only thing worth adding is, FT can give symbolic expression as output, and the output can even be Compiled in some cases, but sadly:

symbolicinversion = FT[test, t];

cf = Compile[t, #] &@symbolicinversion;

With[{eps = 10^-3.}, 
   rstFTcompiled = Table[{t, cf@t}, {t, eps, 3 + eps, 1/25}]]; // AbsoluteTiming
(* {0.015624, Null} *)

ListLinePlot[rstFTcompiled, PlotRange -> {-1, 1/10}]

Mathematica graphics

test isn't the case.

inverseLaplaceHyperbolic trial

A quick test shows that, inverseLaplaceHyperbolic works as well on test except for very small t. (Yeah, similar to rstFTcompiled. ) The following is my trial. Feel free to point out if I've used the function in improper way.

With[{eps = 10^-3.}, 
   rstHyper = Table[{t, 
      inverseLaplaceHyperbolic[test[s], {s, t}, WorkingPrecision -> 20]}, {t, eps, 
      3 + eps, 1/25}]]; // AbsoluteTiming
(* {0.557155, Null} *)
ListLinePlot[rstHyper, PlotRange -> {-1, 1/10}, Mesh -> All]

enter image description here

ILT trial

A quick test shows that, ILT doesn't seem to work well on test. The following is my trial. Feel free to point out if I've used the function in improper way.

fneval = 250;
With[{eps = 10^-3}, 
   rstILT = {#, ILT[test, #, fneval]}[Transpose] &@
     Range[eps, 3 + eps, 1/25]]; // AbsoluteTiming
ListLinePlot[rstILT, PlotRange -> {-1, 1/10}]

Mathematica graphics

Correct answer by xzczd on January 9, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP