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.
N[Integrate[...]]
. UseNIntegrate[...]
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 theWorkingPrecision
, I get a result with the default/automatic method, but also a warning that it might not be correct. $\endgroup$