A Bessel & Struve functions related integral

Question

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

Answer 1

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:

FindSequenceFunction[
  Rest @ CoefficientList[
    Normal @ Series[
      1/2 (-2 - π BesselI[1, x] StruveL[0, x] + BesselI[0, x] (2 + π StruveL[1, x])),
      {x, 0, 30}],
  x]][n]
(* (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]
 ]

3.989672183*10^-7

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.