2
$\begingroup$

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$:

$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$enter image description here

My questions are as follows:

  • What is going on here?
  • Is there any settings that I can turn on to fix this?
$\endgroup$
  • $\begingroup$ Psst, Beta[] is built-in... $\endgroup$ – J. M. is away Jun 3 '16 at 17:12
  • $\begingroup$ @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. $\endgroup$ – Winther Jun 3 '16 at 17:15
  • 5
    $\begingroup$ 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[]... $\endgroup$ – J. M. is away Jun 3 '16 at 17:17
  • $\begingroup$ @J.M. Thanks I'll keep that in mind. $\endgroup$ – Winther Jun 3 '16 at 17:19
3
$\begingroup$

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]]

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ It's easy to check: use Beta[x, n] instead. $\endgroup$ – J. M. is away Jun 3 '16 at 17:46
  • $\begingroup$ 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! $\endgroup$ – Winther Jun 4 '16 at 1:03
  • $\begingroup$ @Winther Great then! $\endgroup$ – Anton Antonov Jun 4 '16 at 1:06
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Winther Jun 4 '16 at 1:08

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.