There is a convergence problem that arises with the use of specific values for x and a in the incomplete gamma function on p.218-219 of the book. I am currently working with small values for x and a, e.g. x=0.03122 and a=0.1468558, and the algorithm does not converge for these or for values of the same magnitude. (more specifically, the gcf function).


I have tried increasing the number of maximum iterations, but I find that even with 1000000 iterations, the algorithm will not converge for some values. Furthermore, 10000000 iterations would take a considerable amount of time on an average machine, since I am using the function in a loop.

Has anybody encountered the same problem, and if so, how we could go about it?

Thanks

Note : I am using a Visual Basic version of this function

If you are using gcf, the code will work only for the cases
where x < a+1.0. Otherwise use gser.
I use the C++ code but I would assume that VB includes the
series representation code for the incomplete Gamma function.
I hope this is useful.

Thanks for your message Clarence.

I am using both functions (gcf and gser) according to the values of x and a. I don't know why the algorithm fails to converge for certain sets of values.

I have now found an alternative algorithm in Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996) which works ok.

However, I am still curious to find out what causes the problem...therefore any further comments welcome!

I am no expert but I think that NR provides 2 functions
for the incomplete gamma function depending on whether
you want what they call P(a,x) or Q(a,x)=1.0-P(a,x).
These are gammp and gammq. In each of these the use the
continued fraction or the series representation depending
on the values of x and a. I ASSUME that the regions of
convergence are different for the two methods and that
the authors have checked that all of the cases are covered
but the code. Obviously, you have found out differently.
I would be interested inhaving some of your data for the
cases where these fail since they claim that all is OK
for all positive values of x and a.

I know this isn't much help but I am curious as well as to what is
happening.

Clarence

Hi,

Clarence gave you the correct answer. You can only use
gcf for x > 1+a, gser otherwise. They have different
regimes of convergence. For the values of x and a you
gave, you must use gser.

Saul Teukolsky