0
$\begingroup$

I have a mathematica code for solving numerical inverse Laplace transform (Credit to Mr. Patrick O. Kano), but sadly the code is using Nintegrate, and i want to use Simpson's rule with adjustable stepsize, and function. By the way my numerical ILT is based on Talbot's method. For more details about Talbot's method, parametrization and its contour, i've mentioned it on my previous question here: Need help to plot Talbot's Contour (Newbie question)

Hope this doesn't sound complicated. So here is my Mathematica code for numerical ILT.

Timeval = 2;
Rval = 2;
Flap[s_] = 1/(s - 1);
Tfunexact[t_] = Exp[1*t];
Valexact = N[Tfunexact[Timeval], 10]

STalbot[r_, a_] = r*a*Cot[a] + I*r*a;
dsda[r_, a_] = (1 + I*(a + Cot[a]*(a*Cot[a] - 1)));
TimeDfun[r_, t_] = 
  Rval/(Pi)*
   NIntegrate[
    Exp[STalbot[r, a]*t]*Flap[STalbot[r, a]]*dsda[r, a], {a, 0, Pi}, 
    WorkingPrecision -> 10];

{Timeval, Approxval} = Timing[TimeDfun[Rval, Timeval]]
RelError = Abs[Approxval - Valexact]/Valexact

Notice the Nintegrate there, i want to replace it with Simpson's rule, and here is my attempt [New Edit]:

simp[r_, t_, 
  M_] := (r/(3 M)) (r*Exp[r*t]*Flap[r] + 
    2*Sum[Exp[
       t*STalbot[r, 2 i*Pi/M]*Flap[STalbot[r, 2 i*Pi/M]]*
        dsda[r, 2 i*Pi/M]], {i, 1, M/2 - 1}] + 
    4*Sum[Exp[
       t*STalbot[r, (2 i - 1)*Pi/M]*Flap[STalbot[r, (2 i - 1)*Pi/M]]*
        dsda[r, (2 i - 1)*Pi/M]], {i, 1, M/2 }])

Flap[s_] := 1/(s - 1);
TableForm[Table[{M, N[simp[0, Pi, M]]}, {M, 10, 100, 10}], 
 TableHeadings -> {{}, {"M", "s_n"}}]

{Timeval, Approxval} = 
 Timing[Re[simp[Rval, Timeval, 50, WorkingPrecision -> 10]]]
RelError = Abs[Approxval - Valexact]/Valexact

Can you please help me to accomplish this. Now, my code shows the result, but still the result was incorrect. Please comment if there's something unclear about my post. Thanks in advance!

$\endgroup$
17
  • 2
    $\begingroup$ I'm sorry, but this site isn't a free coding service, and this question is out of scope of this site. The translation should be straightforward, please have a try yourself. Also, consider using Gauss-Legendre quadrature: mathematica.stackexchange.com/a/6966/1871 $\endgroup$
    – xzczd
    Commented Feb 24, 2021 at 11:20
  • $\begingroup$ @xzczd ok, i'll try. But please don't vote to close this question. I'll be right back to improve my question. I'm not sure and still need some times to understand translating "for while" in Mathematica. And if there's a case who want to help me. Thanks for the advice btw. $\endgroup$
    – user516076
    Commented Feb 24, 2021 at 11:44
  • $\begingroup$ What you need is not For or While. As a start, have a look at Sum. $\endgroup$
    – xzczd
    Commented Feb 24, 2021 at 12:31
  • $\begingroup$ @xzczd Hi. I've edited my question with an attempt. But i still failed. Hope you can help me. $\endgroup$
    – user516076
    Commented Feb 24, 2021 at 12:37
  • 1
    $\begingroup$ Just move the bounds a little to e.g. $[10^{-6},\pi+10^{-6}]$ when using the usual Simpson. If the $re^{rt}\hat f(r)$ in the original paper is designed for trapezoid rule, I won't be surprised if it cannot be naively applied to Simpson's rule. $\endgroup$
    – xzczd
    Commented Feb 25, 2021 at 9:05

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.