# NIntegrate fails in a very strange way

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?
• Psst, Beta[] is built-in... – J. M. is away Jun 3 '16 at 17:12
• @J.M. Yes that does works as a solution to this particular problem. I was wondering more in general how to solve such problems) and I still want to know what is going on here. It's really strange. – Winther Jun 3 '16 at 17:15
• It's usually a bad idea to compute ratios of gamma functions (or factorials, for that matter) in this way for large arguments. Either you re-express everything in terms of LogGamma[] and exponentiate at the very end, or figure out an alternative expression using any of Pochhammer[], Binomial[], FactorialPower[], Beta[]... – J. M. is away Jun 3 '16 at 17:17
• @J.M. Thanks I'll keep that in mind. – Winther Jun 3 '16 at 17:19

## 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]] &)]
]


• It's easy to check: use Beta[x, n] instead. – J. M. is away Jun 3 '16 at 17:46
• Thanks for your answer and the interesting features you mention that I was not aware of like IntegrationMonitor. WorkingPrecision does fix it, should have though of that! – Winther Jun 4 '16 at 1:03
• @Winther Great then! – Anton Antonov Jun 4 '16 at 1:06

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.

• Thank you for your answer. That sounds like a plausible reason. Also nice to learn about the use of NumericQ in the variable definition, was unaware of that you could do that. – Winther Jun 4 '16 at 1:08