4

Happy new year! :) Recently, I am practicing laplace transform technique for solving pdes.

In this extemely helpful post, @xzczd mentioned that "The last step is to transform the solution back, but sadly InverseLaplaceTransform can't handle tsol.".

So, I tried to numerically do inervese laplace transfrom using both built in command InversLaplaceTransform and external package GWR command, but the code is so slow. Also, the result seems not correct when comparing the answer therein.

So, where is the problem. Am i using command wrong? Below is my code. Thanks.

Remove["Global`*"] // Quiet;
<< NumericalLaplaceInversion.m;

(LaplaceTransform[x'[t],t,s])
(s LaplaceTransform[x[t],t,s]-x[0])


f[x_] = x (1 - x);


pde = D[u[t, x], {t, 2}] + D[u[t, x], {x, 4}] == 0;
ic = {u[0, x] == f@x, Derivative[1, 0][u][0, x] == 0};
bc = {u[t, 0] == 0, u[t, 1] == 0, Derivative[0, 2][u][t, 0] == 0, 
   Derivative[0, 2][u][t, 1] == 0};
odeUgly = LaplaceTransform[{pde, bc}, t, s] /. Rule @@@ ic;
ode = odeUgly /. 
    HoldPattern@LaplaceTransform[a_[t, x_], __] :> a[s, x] // 
   Flatten;
(s is a parameter, not a variable, no its derivatives
u[s,x]===U[x]
)
odesol = u[s, x] /. First@DSolve[ode, u[s, x], x] // Simplify;
(odesol=DSolveValue[ode,u[s,x],x]//Simplify)


uatx[x_] = Function[s, odesol // Evaluate];
(uatx[z])


(usolGWR[ti_,xi_]:=GWR[uatx[xi//Rationalize],ti]//Chop;
datGWR=Table[{x,usolGWR[0.2,x]},{x,$MachineEpsilon,1,0.1}]
ListLinePlot[datGWR,PlotRange[Rule]{{0,1},{-0.3,0.3}},Mesh->All])


usol[ti_, xi_] := 
  InverseLaplaceTransform[uatx[Rationalize@xi][s], s, N@ti] // Chop;
dat = Table[{x, usol[0.2, x]}, {x, 0, 1, 0.1}] // Quiet;
ListLinePlot[dat, PlotRange -> {{0 - 0.01, 1 + 0.01}, {-0.3, 0.3}}, 
 Mesh -> All]
xzczd
  • 65,995
  • 9
  • 163
  • 468
xinxin guo
  • 1,323
  • 9
  • 11
  • The DSolveValue[……] line is wrong. After transforming teqn is actually an BVP of ODE, so it should be DSolveValue[……, x]. But after the correction DSolveValue returns unevaluated at least in v12.3.1 (undoubtedly a bug), so, use DSolve instead. – xzczd Jan 06 '22 at 02:29
  • Oh, thanks. Sorry for the silly mistake. I have learned so much from you in this forum. I corrected the code and everything seems to be right now. However, the speed is still low, is it possilbe to accelerate it? How about the compile trick? – xinxin guo Jan 06 '22 at 04:11

1 Answers1

4

Seems that the automatic method choosing of InverseLaplaceTransform isn't working well, at least in v12.3.1. Let's test the speed of the available methods:

AbsoluteTiming@
   InverseLaplaceTransform[uatx[0.2][s], s, 0.1, Method -> #] & /@ {"Stehfest", 
  "Piessens", "Talbot", "Durbin", "Crump", "Papoulis", "Weeks"}
(*
{{8.25883,   0.0760733 + 5.4345*10^-17 I}, 
 {0.501488,  0.0758119 - 1.91165*10^-15 I}, 
 {0.0237084, 0.075812}, 
 {0.47879,   0.0743344}, 
 {0.134007,  0.0846437}, 
 {0.132482,  0.100495 - 0.000824042 I}, 
 {5.25566,   0.0760956}}
 *)

Hmm, "Talbot" method looks like a much better choice, let's turn to it:

usol[ti_, xi_] := InverseLaplaceTransform[uatx[xi][s], s, N[ti], Method -> "Talbot"]

ListLinePlot[Table[usol[0.2, x], {x, 1/10^8, 1, 0.01}], 
  DataRange -> {10^-3, 1}] // AbsoluteTiming

Mathematica graphics

xzczd
  • 65,995
  • 9
  • 163
  • 468