9
$\begingroup$

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)));
$\endgroup$
  • $\begingroup$ What happens if you simply write out the InverseLaplace integral with your integrand (function)? $\endgroup$ – David G. Stork Nov 5 '15 at 21:47
  • $\begingroup$ I did some work on this; I really should find the time to implement the extended version of the algorithm there... $\endgroup$ – J. M. is in limbo Nov 5 '15 at 21:49
  • $\begingroup$ @DavidG.Stork there is no output. Mathematica is running and never stops. $\endgroup$ – Maria Nov 5 '15 at 22:19
  • 1
    $\begingroup$ Many (many many) years ago I used a numerical package that is available in the Wolfram library: library.wolfram.com/infocenter/MathSource/2691 $\endgroup$ – Mike Honeychurch Nov 5 '15 at 22:50
  • 1
    $\begingroup$ Also, you can try this package (BTW I've used the package in this answer.) $\endgroup$ – xzczd Nov 6 '15 at 1:57
7
$\begingroup$

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

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

$\endgroup$
  • $\begingroup$ I'd appreciate if someone tell me how to use the GWR method to solve the above mentioned problem with the parameter kon varying between 1 and 10. So, i need to find the solution as a function of time and kon. Thank you in advance for your reply. $\endgroup$ – JoseOM Apr 21 '19 at 5:11
  • $\begingroup$ @JoseOM Just check the documents of Table and Interpolation. $\endgroup$ – xzczd Apr 21 '19 at 13:15
  • 1
    $\begingroup$ "inverseLaplaceHyperbolic works as well on test except for very small t. (Yeah, similar to rstFTcompiled.)" - not too surprised, both are very similar contour integral methods, except that the former uses a contour that is cheaper to evaluate. $\endgroup$ – J. M. is in limbo Nov 17 '19 at 0:20

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