3
$\begingroup$

I try to plot the following function:

Plot3D[Gamma[1+0.5*(n+m)]/Sqrt[Gamma[1+n]*Gamma[1+m]],{n,0,1000},{m,0,1000}]

I expect that for m=n and near to it the value should be equal to 1. This is true but not for a particular range of m and n. See WA response for query

plot (Gamma(1+(x+T)/2)/sqrt(Gamma(1+x)Gamma(1+T)))  x=0 to 1000 T=0 to 1000

Is it possible that in this range the function is equal to one ONLY for m=n and so the plot3d can not be able to resolve the maximum ?

Thanks for the help

|improve this question
$\endgroup$
  • $\begingroup$ I am guessing, comparing your MMA code with your WA code that you meant Sqrt[Gamma[1+n]*Gamma[1+m]] instead of Sqrt[Gamma[1+n]*Gamma[1+n]] $\endgroup$ – Bill Mar 26 '15 at 16:45
  • $\begingroup$ @Bill you was right I made a type error and I correct it. Thanks, but the problem persist . $\endgroup$ – Panichi Pattumeros PapaCastoro Mar 26 '15 at 16:50
  • $\begingroup$ Another possible workaround would be to reformulate your function entirely in terms of LogGamma[]. $\endgroup$ – J. M. will be back soon♦ Jul 30 '15 at 17:32
6
$\begingroup$

Gamma is an extremely quickly increasing function, so you're dealing with the ratio of huge numbers here. Something similar to catastrophic cancellation can happen.

Fortunately, Mathematica is very good at dealing with this situation if you let it use arbitrary precision instead of machine precision. Change 0.5 to 1/2 and add something like WorkingPrecision -> 30.

Plot3D[Gamma[1 + (n + m)/2]/Sqrt[Gamma[1 + n]*Gamma[1 + m]], {n, 0, 
  1000}, {m, 0, 1000}, PlotRange -> All, WorkingPrecision -> 30]

Mathematica graphics

|improve this answer
$\endgroup$
9
$\begingroup$

While @Szabolcs has provided the correct work-around, the actual issue has to do with automatic use of Compile for function evaluation. The affected range is such that the value of Gamma[1+n] is still a machine real, but their product is out-of-bound.

The issue comes, because the default setting of Compile's RuntimeOptions is to tolerate machine arithmetic underflows. 

(* In[48]:= *)Compile[{{a, _Real}, {n, _Real}, {m, _Real}}, {a, 
   1/Sqrt[n m]}][Gamma[1 + 116.], Gamma[1 + 115.], Gamma[1 + 117.]]

(* Out[48]= {3.39311*10^190, 0.} *)

Product of such numbers gives zero. The work-around is to either use extended precision, as suggested by @Szabolcs, or to avoid using Compile:

Plot3D[Divide[Gamma[1 + (m + n)/2], 
  Sqrt[Gamma[m + 1] Gamma[n + 1]]], {m, 0, 1000}, {n, 0, 1000}, 
 Compiled -> False, PlotRange -> All]
|improve this answer
$\endgroup$
  • $\begingroup$ I didn't know about the Compiled option, thanks for mentioning it! $\endgroup$ – Szabolcs Mar 26 '15 at 18:25
  • $\begingroup$ Fantastic. Thanks a lot :) $\endgroup$ – Panichi Pattumeros PapaCastoro Mar 26 '15 at 18:33
3
$\begingroup$

As I've previously noted, when computing functions that involve ratios of gamma functions, it is manifestly better to re-express in terms of LogGamma[] and then do a final exponentiation afterwards. This deftly sidesteps the issue of huge intermediate values being used to compute a modestly-sized result.

Plot3D[Exp[LogGamma[1 + (n + m)/2] - (LogGamma[1 + n] + LogGamma[1 + m])/2],
       {n, 0, 1000}, {m, 0, 1000}, PlotRange -> All]

plot of a ratio of gamma functions

|improve this answer
$\endgroup$

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.