Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)$$

share|improve this question
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:… 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. – Szabolcs Jan 25 '14 at 22:29
OK. Thanks for your kindness. – A_math_ninja Jan 25 '14 at 22:37
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). – Szabolcs Jan 25 '14 at 22:38
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. – Szabolcs Jan 28 '14 at 16:12
@Szabolcs the questions you see here come from my mind. Do you want me go back into my mind? – A_math_ninja Jan 28 '14 at 16:36
up vote 5 down vote accepted

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.

share|improve this answer
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). – A_math_ninja Feb 2 '14 at 7:56
Do you think the numerical answer can be expressed in terms of the known constants or something similar? – A_math_ninja Feb 2 '14 at 8:00
@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? – Michael E2 Feb 2 '14 at 16:39

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.