NIntegrate fails in a very strange way

Question

I found a strange behavior in Mathematica when trying to evaluate the integral

$$f(n) = \int_{1}^2 \frac{\Gamma(n)\Gamma(x)}{\Gamma(x+n)}{\rm d}x$$

I evaluate this using F[n_] := NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}] For large values of $n$ the two terms $\Gamma(n)$ and $\Gamma(x+n)$ will be huge, but the integrand itself is monotonely decreasing and contained in $[0,\frac{1}{n}]$.

Mathematica (10.3.1.0) is able to evaluate this integral for all values $n\in[1,1000]$ however for $n=171$ the integration breaks down with NIntegrate::izero and gives $0$ as the result:

F[170]  -> 0.00103244
F[171]  -> 0.0 (NIntegrate::izero)
F[172]  -> 0.00105424

If I try non-integer $n$ then I find that it breaks down for all $n \in [170.6,171.6]$ Increasing MinRecursion or AccuracyGoal does not help. The strange part is that Mathematica is able to evaluate the function perfectly well for $x\in[1,2]$ and the plot of the function over this interval looks just like it should. The plot below is $n\cdot\frac{\Gamma(n)\Gamma(x)}{\Gamma(x+n)}$ for $n=171$:

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$

My questions are as follows:

• What is going on here?
• Is there any settings that I can turn on to fix this?

Answer 1

A bit speculative, but I think we can see why 171 is a magic number, if we factor the Gamma[n] from the integral:

 Table[ {n, NIntegrate[Gamma[x] /Gamma[x + n], {x, 1, 2}]} , {n, 168, 173}]

{{168, 7.20924*10^-304}, {169, 4.26119*10^-306}, {170, 2.41873*10^-308}, {171, 0.}, {172, 0.}, {173, 0.}}

171 is where the integral fails as you can see the result is smaller than $MinMachineNumber. I'd guess NIntegrate factors the constant, but stops factoring when the factored value is too large.

note if we hide the structure the problem ( at this n ) goes away:

f[n_?NumericQ, x_?NumericQ] := Gamma[x] Gamma[n]/Gamma[x + n]
Table[{n, NIntegrate[f[n, x], {x, 1, 2}]}, {n, 168, 175}]

{{168, 0.00108399}, {169, 0.00107641}, {170, 0.00106892}, {171, 0.00106153}, {172, 0.00105424}, {173, 0.00104704}, {174, 0.00103993}, {175, 0.00103291}}

This actually (surprisingly) works even for much larger n.

Answer 2

Fixing it ...

Is there any settings that I can turn on to fix this?

Does using larger working precision produce results you expect?

F[n_] := 
 NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, 
  WorkingPrecision -> 60, MinRecursion -> 6, MaxRecursion -> 20]

In[10]:= F[170]

Out[10]= 0.00106892313424594316124137429342057472552302576793438553998338

In[11]:= F[171]

Out[11]= 0.00106153418172547954471366095844979404173963460222797906864894

In[12]:= F[172]

Out[12]= 0.00105424004229628467250748022777662650942896871502716627259853

UPDATE:

Following J.M.'s advice of using Beta we can see that using larger working precision fixes the issue when using Gamma:

FB[n_] := 
 NIntegrate[Beta[x, n], {x, 1, 2}, WorkingPrecision -> 60, 
  MinRecursion -> 6, MaxRecursion -> 20]

Grid[{#, FB[#]} & /@ Range[170, 172]]

What is going on ...

Following george2079's answer and using the option "IntegrationMonitor" we can see what is going on.

With[{n = 171}, 
 Reap@NIntegrate[(Gamma[x] Gamma[n])/Gamma[x + n], {x, 1, 2}, 
   IntegrationMonitor -> (Sow[
       Map[{#1@"Integrand", #1@"Boundaries", #1@"Integral", #1@
           "Error", #1@"GetValues"} &, #1]] &)]
 ]

f[n_?NumericQ, x_?NumericQ] := Gamma[x] Gamma[n]/Gamma[x + n];
With[{n = 171},
 Reap@NIntegrate[f[n, x], {x, 1, 2}, 
   IntegrationMonitor -> (Sow[
       Map[{#1@"Integrand", #1@"Boundaries", #1@"Integral", #1@
           "Error", #1@"GetValues"} &, #1]] &)]
 ]