2
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ This integral does not converge for x->0. Due to Series it goes with x^-3/2 $\endgroup$
    – Akku14
    Nov 13 '21 at 11:53
  • 1
    $\begingroup$ @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}$ $\endgroup$
    – yarchik
    Nov 13 '21 at 11:54
  • 1
    $\begingroup$ @yarchik, o.k., if you use EulerGamma. But OP never says, 0.57721566 is EulerGamma. The rationalized form does not converge. $\endgroup$
    – Akku14
    Nov 13 '21 at 12:01
  • $\begingroup$ Using LogGamma[x + 1.] keeps thing MachinePrecision and fast. (Doesn't fix the convergence problem, but shows why Matlab seems fast.) $\endgroup$
    – Michael E2
    Nov 13 '21 at 16:03
4
$\begingroup$

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    *)
$\endgroup$
1
  • 1
    $\begingroup$ @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. $\endgroup$
    – Roman
    Nov 13 '21 at 20:28
4
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ 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 -> 10`1^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10`6^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10`12^-9}, (LogGamma[x + 1] + EulerGamma*x)/x^(5/2) /. {x -> 10`18^-9}. (+1) $\endgroup$
    – Michael E2
    Nov 13 '21 at 20:19
3
$\begingroup$

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   *)
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – yarchik
    Nov 13 '21 at 18:09
  • $\begingroup$ @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. $\endgroup$
    – Akku14
    Nov 13 '21 at 18:28

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.