1
$\begingroup$

I am trying to do some numerical integration with the following codes:

EI[m_, x_] := x^m/m! (0.5772156649 + Log[x] - HarmonicNumber[m]) + Sum[If[k != m, x^k/((k - m)*k!), 0], {k, 0, Infinity}];
phi[m_, n_, x_] := (-1)^(n + 1)*m^n*Sum[(-Pi^2)^k*Binomial[n, 2*k + 1]*(EI[m, x])^(n - 2*k - 1)*(x^m/m!)^(2*k + 1), {k, 0, Floor[(n - 1)/2]}];
F[n_, m_, t_] := NIntegrate[((1 - Exp[-t*x])*Exp[-n*x]*phi[m, n, x])/x, {x, 0, Infinity}, WorkingPrecision -> 30];

They worked fine when I calculated F[10, 24, 0.4] // Timing. But when $m>24$, for example, F[10, 25, 0.4]// Timing, I got an error saying

NSum::nsnum: Summand (or its derivative) If[k!=25,x^k/((k-25) k!),0] is not numerical at point k = 15.

and forF[10, 26, 0.4]// Timing, an error saying

NSum::nsnum: Summand (or its derivative) If[k!=26,x^k/((k-26) k!),0] is not numerical at point k = 15.

Can anyone help? Thanks.

$\endgroup$
1
$\begingroup$

When defining piecewise functions in Mathematica, If etc. is often less robust compared to Piecewise, and your problem is exactly the case.

With PiecewiseExpand, one can easily translate If into Piecewise:

With[{mid = PiecewiseExpand@If[k != m, x^k/((k - m) k!), 0], 
      mid2 = Rationalize[0.5772156649, 0]}, 
      EI[m_, x_] := x^m/m! (mid2 + Log[x] - HarmonicNumber[m]) + Sum[mid, {k, 0, Infinity}]]

phi[m_, n_, x_] := (-1)^(n + 1) m^n 
   Sum[(-Pi^2)^k Binomial[n, 2 k + 1] (EI[m, x])^(n - 2 k - 1) (x^m/m!)^(2 k + 1), 
       {k, 0, Floor[(n - 1)/2]}]

F[n_, m_, t_] := 
  NIntegrate[((1 - Exp[-t x]) Exp[-n x] phi[m, n, x])/x, {x, 0, Infinity}, 
             WorkingPrecision -> 30, Method -> {Automatic, "SymbolicProcessing" -> False}]

F[10, 26, 4/10] //Quiet // AbsoluteTiming
{2.7180000, 49219.9010173053372129784300811}
$\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.