search

rtoy / maxima

A Clone of Maxima's repo
Other
0 stars 0 forks source link

gamma limit error #3287

Closed rtoy closed 3 months ago

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:10 Created by kcrisman on 2013-08-07 13:34:36 Original: https://sourceforge.net/p/maxima/bugs/2621

Maxima 5.29.1 in ECL. Apparently the limit is actually 1, see https://groups.google.com/forum/#!topic/sage-support/y9UnyWbUagY

(%i2) y:gamma(x+1/2)/(sqrt(x)*gamma(x));
                                           1
                                 gamma(x + -)
                                           2
(%o2)                          ----------------
                               sqrt(x) gamma(x)
(%i3) limit(y,x,0);
(%o3)                                  0
rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:11 Created by pbruin on 2014-05-18 00:19:25 Original: https://sourceforge.net/p/maxima/bugs/2621/#e442

It seems to me that the limit as x -> 0 is correctly computed as 0. The question originally asked in the discussion linked to, however, was about the limit as x -> infinity.

For any real number a, Stirling's approximation shows that the function

f: gamma(x + a)/(x^a*gamma(x));

tends to 1 as x -> infinity. This is also what

limit(x, a, inf);

gives you (after asking a few questions). However, for specific rational but non-integral values of a the result returned by Maxima seems to be either 0 or inf. Besides the above example (using limit(y, x, inf)) one also gets, for example,

(%i23) f: gamma(x - 2/5)/(x^(-2/5)*gamma(x));
                                2/5           2
                               x    gamma(x - -)
                                              5
(%o23)                         -----------------
                                   gamma(x)
(%i24) limit(f,x,inf);
(%o24)                                 0

I tried to do a bit of debugging. It seems that the limit (in the case a = 1/2, say) is computed via limit(g, x, inf), where

g: 1/sqrt(x)*exp(x*log(2*x-1)-(x-1/2)*log(x-1)-1/2)/2^x;

as obtained by Stirling's approximation. Indeed, trying to compute this limit also yields infinity instead of the correct value 1.

When we instead simplify the expression for g to

g1: 1/sqrt(x)*exp(x*log(x-1/2)-(x-1/2)*log(x-1)-1/2);

the limit is still computed as infinity, but this time it takes several minutes. I don't know if the slowness is related to the incorrect answer. Finally, further simplifying g1 to

g2: exp(x*log(x-1/2)-(x-1/2)*log(x-1)-1/2-log(x)/2);

and computing limit(g2, x, inf) does instantaneously return the correct answer 1.

(All of the above is in Maxima 5.33.0 on GCL.)

Thanks,

Peter Bruin

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:15 Created by pbruin on 2014-05-20 14:03:36 Original: https://sourceforge.net/p/maxima/bugs/2621/#2766

I finally discovered the reason for this bug: the series expansion step in Gruntz's algorithm (see his thesis [1]) is not implemented correctly. On page 49, it is stated that when expanding a logarithm, log w must be replaced by h. In Gruntz's Maple implementation (appendix of [1]), he does this by means of a magical declaration ln(W) := e3; that is somehow picked up by Series().

In the Maxima implementation, the substitution log w -> h is only done after expanding the whole expression. The problem is that by time, the log w may have been transformed, e.g. because it appears inside an exponential.

The patch limit-rewrite-logs.patch [see message below] fixes this, at least for "pure exp-log functions" (the ones for which Gruntz's algorithm was designed) by doing a suitable substitution on relevant subexpressions of the form log(f(x)) before the series expansion.

The patch currently keeps the existing substition after the series expansion because (1) if the expression involves other transcendental functions whose Taylor expansion involves logarithms, we still want to transform those if we see them, and (2) not doing this substitution led to errors even in cases where it is apparently unnecessary, e.g. in the test

integrate(x^6/(1 + x + x^2)^(15/2), x, 0, inf);

I hope this is an acceptable fix for this bug.

[1] D. Gruntz, On computing limits in a symbolic manipulation system. Ph.D. thesis, ETH Zürich, 1996, http://e-collection.library.ethz.ch/view/eth:40284

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:18 Created by pbruin on 2014-05-20 14:16:04 Original: https://sourceforge.net/p/maxima/bugs/2621/#e7c0

Here is a slightly simplified version of the patch. [Edit: simplified once more]

Attachments:

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:22 Created by robert_dodier on 2014-05-27 17:42:52 Original: https://sourceforge.net/p/maxima/bugs/2621/#9f63

rtoy commented 3 months ago

Imported from SourceForge on 2024-07-07 19:53:25 Created by robert_dodier on 2014-05-27 17:42:52 Original: https://sourceforge.net/p/maxima/bugs/2621/#af40

Fixed by commit [b0579c08a] (applied patch).