On 1/29/06, yigal <yigal.weinstein at gmail.com> wrote:
> I want to compute the incomplete gamma function, gamma(x,z) via a
> continued fraction. Is there an easy or at least straightforward way in
> Maxima to do this? The sequences for numerators and denominators are
> simple:
>
> num = (1,(a-1)*1,(a-2)*2,...)
> den = (x+1-a, x+3-a, x+5-a,...)
Well, here is a simple-minded approach. Maybe this is enough.
gcf (cf_num, cf_denom,n) := block ([s : 0],
for k from n - 1 step -1 thru 0 do
s : ev(cf_num) / (ev(cf_denom) + s), s);
cf_num is an expression to give the numerators.
cf_denom is an expression to give the denominators.
n is the number of terms you want, counting from 0 to n - 1.
gcf just builds up a nested quotient of the specified depth.
gcf ('(if k=0 then z^a*exp(-z) else k*(a - k)), z - a + 2*k + 1, 5);
=> (nested quotients expression)
gammainc_approx (a, z) := '' %;
gammainc_approx (3/4, 1/2); => (some stuff)
'' %, numer; => 0.5753834162127
Another program (Octave) seems to think that's correct to
about 2 decimal places. I didn't try any other examples.
Note that '' in the stuff above is two single quotes, not a double quote.
'' % means re-evaluate the preceding output.
I looked to see if the maxima/share directory has any CF
stuff, but I didn't find any. I also looked for a built-in incomplete
gamma function, but didn't find one.
Hope this helps,
Robert Dodier