4
$\begingroup$

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]
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Jan 6 at 2:29
  • $\begingroup$ 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? $\endgroup$
    – xinxin guo
    Jan 6 at 4:11
4
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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