Message Boards Message Boards

2
|
3119 Views
|
15 Replies
|
10 Total Likes
View groups...
Share
Share this post:

Question on the Laplace Inversion of a function using Mathematica

Posted 10 months ago

I need to find the following inverse Laplace transform:

$$q(t)=\mathcal L^{-1} \left [\frac{12 e^{24 s} (1 - e^{168 s} + 168 e^{180 s} s)}{ e^{12 s} + e^{168 s} - e^{180 s} + 2016 e^{204 s} s^2-1} \right ](t), $$

12 e^(24 s)(1 - e^(168 s) + 168 e^(180 s) s))/(-1 + e^(12 s) + e^(168 s) - e^(180 s) + 2016 e^(204 s) s^2)

obtained as the solution of a convolution Volterra integral equation; see here for more details.

What I know about the inverse function is that $$q(t)=1, \, $0<t<12$$ and it tends to $16/19$ as $t$ tends to infinity (we always have $ 0\le q(t) \le1$ for $t \ge0$).

WolframAlpha cannot obtain the solution within a (non-pro) computational time, so I guess there is no closed-form solution for this inverse.

I know that there is a package that finds inverse Laplace transforms numerically, but I don't know how to use it in that environment.

Could you derive the final solution and show me how it is obtained?

POSTED BY: Amir A.
15 Replies

I am no expert on Laplace transforms, I just tried out the built-in commands of Mathematica:

f[t_] = InverseLaplaceTransform[(12  (1 - E^(168  s) + 
       168  E^(180  s)  s))/(-1 + E^(12  s) + E^(168  s) - 
     E^(180  s) + 2016  E^(204  s)  s^2), s, t]

gives no closed-form output. I tried numerical values:

In[1986]:= {f[1.], f[10.], f[1000.], 16/19.}

Out[1986]= {1.27295*10^-111, -1.82483*10^-10, 0.8421, 0.842105}

which partially agree with what you know.

POSTED BY: Gianluca Gorni

If I compute first term of asymptotic I have: $\theta (t-24)$

and looks like for: 0<t<24 is zero and for t > 24 is 1.

 func = (12 (1 - E^(168 s) + 168 E^(180 s) s))/(-1 + E^(12 s) + E^(
     168 s) - E^(180 s) + 2016 E^(204 s) s^2);
 g[t_] := InverseLaplaceTransform[func, s, N@t, Method -> "Durbin"];
 f = InverseLaplaceTransform[Asymptotic[func, {s, Infinity, 1}], s, t]
  Show[{Plot[f, {t, 0, 50}, PlotStyle -> Red], 
    ListLinePlot[Table[{t, g[t]}, {t, 1/10, 100, 1/3}]] // Quiet}, 
   PlotRange -> All]

enter image description here

Plot looks correct to me.

Regards M.I

POSTED BY: Mariusz Iwaniuk

I had copied and pasted from your second line. Now it's better:

f[t_] = InverseLaplaceTransform[(12  E^(24 s) 
      (1 - E^(168   s) + 168   E^(180   s)   s))/
    (-1 + E^(12   s) + E^(168   s) -
      E^(180   s) + 2016   E^(204   s)   s^2), s, t];
{f[1.], f[10.], f[1000.], 16/19.}

{1., 0.999996, 0.842102, 0.842105}
POSTED BY: Gianluca Gorni

Which method is the best (you already used Durbin)?

This is a difficult question and I have no competence and I am not an expert in this field.In my humble opinion, it depends on the given function what method is the best.You just have to experiment yourself.

See attached file.

Attachments:
POSTED BY: Mariusz Iwaniuk
Posted 10 months ago

I have used GWR for many years, with very good results. This is the illusttration how to use it:

(*http://library.wolfram.com/infocenter/MathSource/4738/*)GWR[F_,t_,M_: 32,precin_: 0]:=Module[{M1,G0,Gm,Gp,best,expr,τ=Log[2]/t,Fi,broken,prec},If[precin<=0,prec=21 M/10,prec=precin];
If[prec<=$MachinePrecision,prec=$MachinePrecision];
broken=False;
If[Precision[τ]<prec,τ=SetPrecision[τ,prec]];
Do[Fi[i]=N[F[i τ],prec],{i,1,2 M}];
M1=M;
Do[G0[n-1]=τ (2 n)!/(n! (n-1)!) Sum[Binomial[n,i] (-1)^i Fi[n+i],{i,0,n}];
If[Not[NumberQ[G0[n-1]]],M1=n-1;G0[n-1]=.;Break[]];,{n,1,M}];
Do[Gm[n]=0,{n,0,M1}];
best=G0[M1-1];
Do[Do[expr=G0[n+1]-G0[n];
If[Or[Not[NumberQ[expr]],expr==0],broken=True;Break[]];
expr=Gm[n+1]+(k+1)/expr;
Gp[n]=expr;
If[OddQ[k],If[n==M1-2-k,best=expr]];,{n,M1-2-k,0,-1}];
If[broken,Break[]];
Do[Gm[n]=G0[n];G0[n]=Gp[n],{n,0,M1-k}];,{k,0,M1-2}];
best]
SetAttributes[GWR,{Listable,NHoldAll}]

F[s_]:=12 (1-Exp[168 s]+168 Exp[180 s] s)/(-1+Exp[12 s]+Exp[168 s]-Exp[180 s]+2016 Exp[204 s] s^2);

t={10^-1,1,10,100,1000,10000,100000};

N[GWR[F,t,60],20]
N[GWR[F,t,80],20]
N[GWR[F,t,100],20]
N[GWR[F,t,120],20]
N[GWR[F,t,140],20]

You should increase the number of digits (the third calling parameter of GWR) until you get a satisfactory accuracy of the results. I notice that your original formula has one extra parenthesis, so that I am not sure if after removing it the formula for F is correct. But if it is, then I see that your transform may require a rather high value of the third parameter. Of course, instead of using a list of time values (t) you can just put a single value. Leslaw

POSTED BY: Leslaw Bieniasz

See attached file.

Attachments:
POSTED BY: Mariusz Iwaniuk

See attached file.

Attachments:
POSTED BY: Mariusz Iwaniuk
Posted 10 months ago

Gianluca thank you so much for your reply! The error for $t \le 25$ is very high, but for t>25, the inverse seems reasonable. For $0<t<12$, the function needs to be 1, while it is 0 and in some cases is negative! Hope some Mathematica developers let us know how this can de fixed.

POSTED BY: Amir A.
Posted 10 months ago

Thank you so much Mariuz! I think you and Gianluca ignored $e^{24 s}$ in the numerator. Could you please make sure you have used the function given in OP, and see what will happen next? It is also given here.

Which method is the best (you already used Durbin)?

POSTED BY: Amir A.
Posted 10 months ago

Thank you so much Leslwa for the input! You may see that the formula is

12 e^(24 s)(1 - e^(168 s) + 168 e^(180 s) s))/(-1 + e^(12 s) + e^(168 s) - e^(180 s) + 2016 e^(204 s) s^2)

Could you check the result in this case?

Do you think this package developed by Peter Valkó provides more accurate results compared to those used in other replies. Have you ever checked it?

POSTED BY: Amir A.
Posted 10 months ago

Mariusz Thanks a ton! I am still struggling to run the code in a notebook in ww.wolframcloud.com, but I have been unsuccessful. The following can be compiled (while I could not plot f[t_])

f[t_] = InverseLaplaceTransform[(12 E^(24 s) (1 - E^(168 s) + 168 E^(180 s) s))/ (-1 + E^(12 s) + E^(168 s) - E^(180 s) + 2016 E^(204 s) s^2), s, t].

However, I am getting the error "No Wolfram Language translation found" with your code (I guess I should work in another environment and sorry if I do not have sufficient experience on working with Wolfram environments and Mathematica). To make sure that the output is correct, could you please provide a figure of $p(t)=1-q(t)$ over t= 1 to 1500 (at time 1, 2, 3, ...., 400, and 450, 500, ..., 1500) using "Weeks" method (which seems to be one of most stable methods)? I guess the function have significant fluctuations before 400 ( $p(t)=0$ for $0<t<12$ and tends to 3/19 in at infinity).

POSTED BY: Amir A.
Posted 10 months ago

I appreciate it Gianluca!

POSTED BY: Amir A.
Posted 10 months ago

I have verified the GWR package many times. You should run it yourself, otherwise you will not get any experience with it. One thing that worries me, regarding your transform, is that due to the presence of exponential terms it can be rather sensitive to variations of s. There is also a subtraction of exponential terms, which for large s (corresponding to small t) can possibly cause a considerable loss of significant digits. It may be necessary to do some rearrangements of your transform, if you need accurate results, but I have no concrete suggestions how to rearrange. Leslaw

POSTED BY: Leslaw Bieniasz
Posted 10 months ago

Wow! A nice comparison between Weeks method and GWR package. Thank you! Could you kindly also include the plot of $1-q(t)$ in the file as I already can see the file in Wolfram cloud, but cannot run it there.

POSTED BY: Amir A.
Posted 10 months ago

Thanks! I just updated my answer here. You may check it.

POSTED BY: Amir A.
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract