0

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_] = raCot[a] + Ira; dsda[r_, a_] = (1 + I(a + Cot[a](aCot[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!

user516076
  • 373
  • 1
  • 8
  • 2
    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: https://mathematica.stackexchange.com/a/6966/1871 – xzczd Feb 24 '21 at 11:20
  • @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. – user516076 Feb 24 '21 at 11:44
  • What you need is not For or While. As a start, have a look at Sum. – xzczd Feb 24 '21 at 12:31
  • @xzczd Hi. I've edited my question with an attempt. But i still failed. Hope you can help me. – user516076 Feb 24 '21 at 12:37
  • @xzczd Gauss-Legendre is good. But at this time, i just want to use simpson even if it might not that good. Cz i've seen some people use that and i want something new even if it's not effective. – user516076 Feb 24 '21 at 12:39
  • You haven't typed simpson's rule correctly, please double-check it. 2. \[Theta] and f should be defined as argument of function, too.
  • – xzczd Feb 24 '21 at 12:49
  • @xzczd what about now? See my edit. Sorry i'm too dumb. Still giving no result. – user516076 Feb 24 '21 at 13:13
  • It's even worse. Please calm down, and think about how to code simpson's rule in a general way. If you're still having difficulty, observe the behavior of the following sample, and think carefully about what you can learn from it: trig30degree[f_] := f[30 Degree];{trig30degree[Sin], trig30degree[Cos], trig30degree[Tan]} – xzczd Feb 24 '21 at 13:21
  • @xzczd is it something like composite function? – user516076 Feb 24 '21 at 13:30
  • The key point is, a function relationship can be the argument of a function. Further more: trig30degree[Function[x, Sin[x]+Cos[x]]] – xzczd Feb 24 '21 at 14:20
  • Still, you haven't typed the formula of simpson's rule correctly, please double-check it, pressing Ctrl+Shift+N to transform your code to StandardForm may help you to see the mistake. – xzczd Feb 25 '21 at 03:49
  • @xzczd can you tell me which part is false. I'm really confuse and have no idea. Thanks. – user516076 Feb 25 '21 at 03:54
  • No. Please compare what you have typed with the formula in e.g. wikipedia carefully. – xzczd Feb 25 '21 at 04:01
  • @xzczd I'm using composite simpson's rule 1/3 sir. I've read the wikipedia 3 times but where's my mistake? I'm 100% sure it looks exactly the same with wikipedia. About $f(b)+f(a)$ i write the function with the similar form in paper that is $re^{rt}\hat f(r)$. In that case i'm using $r$ maybe this is why my formula little bit different? – user516076 Feb 25 '21 at 04:23
  • You should first make sure you can code the classic Simpson's rule before making modification. 2. If you're not coding the classic Simpson's rule, proper reference should be given, or others won't even be able to check.
  • – xzczd Feb 25 '21 at 05:03
  • @xzczd Thanks for the advice. Somehow if i put $f(a)$ and $f(b)$ as in usual composite simpson it makes the integration diverges since $\cot(\pi)=-\infty$. That's my confusion. Using the formula above with the trapezoidal's rule gives me the correct answer but i still don't understand why i can't apply to simpson. Maybe i should give up. – user516076 Feb 25 '21 at 09:02
  • 1
    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. – xzczd Feb 25 '21 at 09:05