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!
For
orWhile
. As a start, have a look atSum
. $\endgroup$