I try to numerically compute this integral and I don't figure out why on earth Mathematica is not able
to do it. Is my input correct? Does it possibly have a closed form?

N[Integrate[1/2 (-2 - π BesselI[1, x] StruveL[0, x] + 
  BesselI[0, x] (2 + π StruveL[1, x]))/(E^x - 1), {x, 0, Infinity}]]

EDIT: the correct answer is 0.6076895380471122

OK, here is the infinite series representation of the above integral

$$\sum_{k=1}^{\infty} \left(\arcsin\left(\frac{1}{k}\right)-\frac{1}{k}\right)$$

  • $\begingroup$ I can't answer the question, but a few observations might help: 1. if you need to integrate numerically, don't use N[Integrate[...]]. Use NIntegrate[...] and check the available methods. A very long section of the docs is dedicated to this: reference.wolfram.com/mathematica/tutorial/… 2. if I increase the WorkingPrecision, I get a result with the default/automatic method, but also a warning that it might not be correct. $\endgroup$ – Szabolcs Jan 25 '14 at 22:29
  • $\begingroup$ OK. Thanks for your kindness. $\endgroup$ – user 1357113 Jan 25 '14 at 22:37
  • $\begingroup$ Trying to plot the integrand for a wide interval, e.g. $(0,100)$, reveals that in its current form it's impossible to compute the integrand using machine precision arithmetic only (because of the ratio of huge numbers that appears). So using a higher WorkingPrecision is going to be essential (or rewriting the integrand to avoid this problem would be even better, if possible). $\endgroup$ – Szabolcs Jan 25 '14 at 22:38
  • $\begingroup$ Your best bet would be to go back to where you got the integrand from in the first place and try to rewrite this integral in a form that avoids dividing or subtracting huge numbers to get a result of order 1 magnitude. $\endgroup$ – Szabolcs Jan 28 '14 at 16:12
  • $\begingroup$ @Szabolcs the questions you see here come from my mind. Do you want me go back into my mind? $\endgroup$ – user 1357113 Jan 28 '14 at 16:36

I think the reason NIntegrate is hard to use is that the two terms that are products of a Bessel and a Struve are very large numbers nearly equal numbers being subtracted and perhaps that it the integral converges slowly. Since the OP already knows the answer, it makes me wonder about the purpose of the question, whether it is about how to find the answer or why NIntegrate fails to work. If the former, perhaps the following will be useful. Since Mathematica can integrate x^n/(E^x - 1), we can look for a power series solution:

Assuming[{n > 0 && n \[Element] Integers},
 Integrate[x^n/(E^x - 1), {x, 0, Infinity}]
(* n! Zeta[1 + n] *)

The power series coefficients may be determined from FindSequenceFunction:

  Rest @ CoefficientList[
    Normal @ Series[
      1/2 (-2 - π BesselI[1, x] StruveL[0, x] + BesselI[0, x] (2 + π StruveL[1, x])),
      {x, 0, 30}],
(* (2^(-1 - n) (1 + (-1)^n))/((1 + n) Gamma[1 + n/2]^2) *)

The odd-power terms are zero, so we can replace n by 2k. Putting it all together, we get that the integral is equal to

NSum[(4^-k (2 k)! Zeta[1 + 2 k])/((1 + 2 k) Gamma[1 + k]^2), {k, Infinity}]
(* 0.607689461628809` *)

It agrees with the OP's answer to within 10^-7.

The integrand

Just to give an idea of the problematic nature of the integrand: The values get so large that numerical approaches will probably run into serious obstacles. For example, the two large terms of the integrand agree to over 4000 digits at x -> 10^4.

N[{π BesselI[1, x] StruveL[0, x], 
   BesselI[0, x] (2 + π StruveL[1, x])} /. x -> 10^4, 10]
{3.877905409*10^8681, 3.877905409*10^8681}
Block[{$MaxExtraPrecision = 4300},
  N[Differences @
     {π BesselI[1, x] StruveL[0, x], 
      BesselI[0, x] (2 + π StruveL[1, x])} /. x -> 10^4, 10]
  ] // AbsoluteTiming

N::meprec: Internal precision limit $MaxExtraPrecision = 4300.` reached while evaluating {-[Pi] BesselI[1,10000] StruveL[0,10000]+BesselI[0,10000] (2+[Pi] StruveL[1,10000])}. >>

{36.431719, {0.*10^4375}} 
Block[{$MaxExtraPrecision = 4400},
  N[Differences @
     {π BesselI[1, x] StruveL[0, x], 
      BesselI[0, x] (2 + π StruveL[1, x])} /. x -> 10^4, 10]
  ] // AbsoluteTiming
{34.786213, {7.027263539*10^4336}}

The difference is still fairly significant:

Block[{$MaxExtraPrecision = 4400},
 N[1/2 (-2 - π BesselI[1, x] StruveL[0, x] + 
       BesselI[0, x] (2 + π StruveL[1, x]))/(E^x - 1) /. 
   x -> 10^4, 10]

The function appears to be decreasing and convex up and the derivative at this point is 5.984757694*10^-11. It follows from these suppositions that the integral from 10^4 to Infinity would be at least 0.0013 just from computing the area below by the tangent line.

  • $\begingroup$ Thank you for confirming that!!! Many of the mathematicians I usually talk to said I'm far away from the truth (maybe also because I have no background in mathematics they don't give me much credit). $\endgroup$ – user 1357113 Feb 2 '14 at 7:56
  • $\begingroup$ Do you think the numerical answer can be expressed in terms of the known constants or something similar? $\endgroup$ – user 1357113 Feb 2 '14 at 8:00
  • 2
    $\begingroup$ @Chris'ssis Off the top of my head, I don't have an idea how to get to an exact answer. Out of curiosity, how did you get the answer you give in the question? $\endgroup$ – Michael E2 Feb 2 '14 at 16:39

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.