# Integral using LogGamma fails

For the Matlab code:

integral(@(x)(gammaln(x+1)+0.57721566 *x)./x.^2.5,0,1e9)


I get =5.469... (almost instantaneously)

Trying the equivalent in Mathematica 12.3:

NIntegrate[(LogGamma[x+1]+0.57721566*x)/x^2.5,{x,0,1*10^9}]


the calculation would not finish after a long wait.

Will you show me how to obtain the Matlab (correct) answer?

• This integral does not converge for x->0. Due to Series it goes with x^-3/2 Nov 13 '21 at 11:53
• @Akku14 This is not true, Series[(LogGamma[x + 1] + EulerGamma x)/x^(5/2), {x, 0, 3}] yields the leading term of $\frac{1}{360} \pi ^4 x^{3/2}+\frac{\pi ^2}{12 \sqrt{x}}+\frac{\sqrt{x} \psi ^{(2)}(1)}{6}$ Nov 13 '21 at 11:54
• @yarchik, o.k., if you use EulerGamma. But OP never says, 0.57721566 is EulerGamma. The rationalized form does not converge. Nov 13 '21 at 12:01
• Using LogGamma[x + 1.] keeps thing MachinePrecision and fast. (Doesn't fix the convergence problem, but shows why Matlab seems fast.) Nov 13 '21 at 16:03

The integrand, with all numbers exact:

f[x_] = (LogGamma[x + 1] + EulerGamma x)/x^(5/2);


Capture the singular behavior for small $$x$$ in the asymptotic expression for $$x\to0^+$$:

g[x_] = Series[f[x], {x, 0, 1}] // Normal
(*    π^2/(12 Sqrt[x]) + 1/6 Sqrt[x] PolyGamma[2, 1]    *)


Integrate parts separately and add up: the singular integral is treated analytically in the first term,

Integrate[g[x], {x, 0, 1}] +
NIntegrate[f[x] - g[x], {x, 0, 1}] +
NIntegrate[f[x], {x, 1, ∞}]

(*    5.47135    *)


If you need more precision:

With[{wp = 10^2},
Integrate[g[x], {x, 0, 1}] +
NIntegrate[f[x] - g[x], {x, 0, 1}, WorkingPrecision -> wp] +
NIntegrate[f[x], {x, 1, \[Infinity]}, WorkingPrecision -> wp]]

(*    5.471346135899603136108128448540649799821138398820913558897482509129509200122440032933242531582252312    *)

• @yarchik I took the integration limit of $10^9$ to mean "integrate to infinity", in the same way as I took $0.57721566$ to be the Euler–Mascheroni constant. Yes, if the OP is actually interested in integrating only to $10^9$, other tricks may be needed, for example what you propose. Nov 13 '21 at 20:28

You can get an answer almost instantly with

i1 = NIntegrate[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2), {x, 0, Infinity}, WorkingPrecision -> 12, MaxRecursion -> 20];
i2 = NIntegrate[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2), {x, 10^9, Infinity}, WorkingPrecision -> 12, MaxRecursion -> 20];
i1 - i2
(* 5.46993572961 *)


Notice that by explicitly specifying the WorkingPrecision the singularity problem at x==0 is avoided.

Finally, I would like to add why there is a problem at zero with machine precision. Compare

N[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10^-9}]
N[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10^-9}, 1]
N[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10^-9}, 6]
N[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10^-9}, 12]
N[(LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10^-9}, 18]
(* -1.48426*10^6 *)
(* 3.*10^4 *)
(* 26008.7 *)
(* 26008.6912475 *)
(* 26008.6912475495376 *)


As you can see, all arbitrary precision results are consistent, whereas machine precision is not.

• I believe NIntegrate does not use N[] this way, but does something like this (probably while blocking $MinPrecision and $MaxPrecision): (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10.^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 101^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 106^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 1012^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 1018^-9}. (+1) Nov 13 '21 at 20:19

This is not an issue of Mathematica (rather than MatLab)

If you use 0.57721566 as constant, the integrand does not converge for x-> 0 (even not the rationalized form). Only for EulerGamma beeing the constant, it converges as Series shows.

f1 = (LogGamma[x + 1] + 0.57721566 x)/x^2.5;

ser1 = Series[f1, {x, 0, 1}] // Normal // Expand

(*   -(4.90153*10^-9/x^1.5)   *)

Integrate[ser1, {x, 0, 10^9}]

(*   Integrate::idiv:  "Integral of 1/x^1.5 does not converge on {0,1000000000}"  *)


NIntegrate explodes.

f2 = (LogGamma[x + 1] + EulerGamma x)/x^(5/2);

ser2 = Series[f2, {x, 0, 1}] // Normal

(*   \[Pi]^2/(12 Sqrt[x]) + 1/6 Sqrt[x] PolyGamma[2, 1]   *)

Integrate[ser2, {x, 0, 10^9}]

(*   5000/9 Sqrt[10] (3 \[Pi]^2 + 2000000000 PolyGamma[2, 1])   *)


With given workingprecision get the result at once.

NIntegrate[f2, {x, 0, 10^9}, WorkingPrecision -> 15]

(*   5.46993572961007   *)

• Just curious, how long does it take to execute your last command? MA11.1.1 runs too long, however, splitting the integral helps as my answer shows. Nov 13 '21 at 18:09
• @yarchik, last command is done in 0.078 seconds with version 8.0, and takes half the time with WorkingPrecision->10. But any attempts to work with MachinePrecision like MichealE2 propsed, or simpliy 1. *f2 fail to converge. Strange. Nov 13 '21 at 18:28