Sum over incomplete Gamma Functions 
Author Message
 Sum over incomplete Gamma Functions

Not really a fortran specific question, but perhaps you have some
clever hints for me anyway:

The result of an integration leads to an infinite sum over incomplete
Gamma functions:

  inf         ( Gamma[ 2n, k*z] - Gamma[ 2n, z] )
Sum     ---------------------------------------------------------
 n=0              n**(2n+1) * n! * (2n+1)!

where 0 < k <= 1 and z > 0, Gamma[.,.] denotes the incomplete Gamma
function in Mathematica notation - the sum is convergent.

However, already for moderate values of k and z, the number of
iterations exceeds the maximum number possible before overflows occurs
and convergence is still far off.

For example z=0.5 and k=1000: convergence in the 11th digit is at
iteration 560 as calculated with Mathematica.

Context: this sum is used in an numerical integration scheme as part
of the integral kernel and called many times.

I consulted some of the usual literature (i.e. Abramowitz&Stegun,
etc.) but there seems to be no specific (semi)closed form for this
sum, or is it?
Has somebody stumbled over this before?

All hints most welcome - Cheers!

The code calls an implementation of the incomplete Gamma function
(Num.Recipes) and for smaller values of kappa and z is equivalent to
Mathematica's output. The factorial is evaluated in double precision.

!------------------------------------------------------------------------------------------------
function GammaSum(kappa,z)
use typedef, only: ft ! double precision
use specfunc, only: Fac, IncGamma ! Factorial and Incomplete Gamma
implicit none
integer   , parameter  :: GAMMAITR = 80
real( ft ), intent(in) :: kappa,z
real( ft )             :: GammaSum
real( ft )             :: gsm,prv,ord,num,den
integer                :: n

gsm  = 0
Do n=0,GAMMAITR

   prv = gsm
   ord = real(2*n,ft)

   num = IncGamma(ord,kappa*z) - IncGamma(ord,z)
   den = 2.0_ft**(2*n+1) * Fac(n)*Fac(n+1)

   gsm = gsm + num/den

  ! Exit if convergence .... (not shown)
EndDo
GammaSum = gsm

end function GammaSum
!------------------------------------------------------------------------------------------------



Wed, 14 Sep 2011 12:15:06 GMT  
 Sum over incomplete Gamma Functions

Quote:

> Not really a Fortran specific question, but perhaps you have some
> clever hints for me anyway:

> The result of an integration leads to an infinite sum over incomplete
> Gamma functions:

>   inf         ( Gamma[ 2n, k*z] - Gamma[ 2n, z] )
> Sum     ---------------------------------------------------------
>  n=0              n**(2n+1) * n! * (2n+1)!

Your code does 2**(2n+1), not n**(2n+1). Which is correct? If the
former, then it's easy to generate each denominator from the previous
one, and avoid computing the factorials and power.

Quote:

> where 0 < k <= 1 and z > 0, Gamma[.,.] denotes the incomplete Gamma
> function in Mathematica notation - the sum is convergent.

> However, already for moderate values of k and z, the number of
> iterations exceeds the maximum number possible before overflows occurs
> and convergence is still far off.

> For example z=0.5 and k=1000: convergence in the 11th digit is at
> iteration 560 as calculated with Mathematica.

> Context: this sum is used in an numerical integration scheme as part
> of the integral kernel and called many times.

> I consulted some of the usual literature (i.e. Abramowitz&Stegun,
> etc.) but there seems to be no specific (semi)closed form for this
> sum, or is it?
> Has somebody stumbled over this before?

> All hints most welcome - Cheers!

> The code calls an implementation of the incomplete Gamma function
> (Num.Recipes) and for smaller values of kappa and z is equivalent to
> Mathematica's output. The factorial is evaluated in double precision.

> !------------------------------------------------------------------------------------------------
> function GammaSum(kappa,z)
> use typedef, only: ft ! double precision
> use specfunc, only: Fac, IncGamma ! Factorial and Incomplete Gamma
> implicit none
> integer   , parameter  :: GAMMAITR = 80
> real( ft ), intent(in) :: kappa,z
> real( ft )             :: GammaSum
> real( ft )             :: gsm,prv,ord,num,den
> integer                :: n

> gsm  = 0
> Do n=0,GAMMAITR

>    prv = gsm
>    ord = real(2*n,ft)

>    num = IncGamma(ord,kappa*z) - IncGamma(ord,z)
>    den = 2.0_ft**(2*n+1) * Fac(n)*Fac(n+1)

>    gsm = gsm + num/den

>   ! Exit if convergence .... (not shown)
> EndDo
> GammaSum = gsm

> end function GammaSum
> !------------------------------------------------------------------------------------------------

--
*********** To reply by e-mail, make w single in address **************


Thu, 15 Sep 2011 01:56:03 GMT  
 Sum over incomplete Gamma Functions

[...]

Quote:
> Your code does 2**(2n+1), not n**(2n+1). Which is correct? If the
> former, then it's easy to generate each denominator from the previous
> one, and avoid computing the factorials and power.

thanks Ian, that's a typo -

Ok, I see your idea and it might be worthwhile to follow up -

Cheers - Ralf



Thu, 15 Sep 2011 06:02:43 GMT  
 Sum over incomplete Gamma Functions

Quote:
> The result of an integration leads to an infinite sum over incomplete
> Gamma functions:
>  inf         ( Gamma[ 2n, k*z] - Gamma[ 2n, z] )
> Sum     ---------------------------------------------------------
> n=0              n**(2n+1) * n! * (2n+1)!

If we take your denominator to be 2**(2*n+1)*n!*(n+1)!, your original
integral would have been:

Integral_k*z^z{I_1(t)*exp(-t)/t**2}*dt,

right?  Why should it be easier to evaluate the series rather
than approximate the integral directly?

--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end



Thu, 15 Sep 2011 07:35:56 GMT  
 Sum over incomplete Gamma Functions


Quote:


> > The result of an integration leads to an infinite sum over incomplete
> > Gamma functions:
> > ?inf ? ? ? ? ( Gamma[ 2n, k*z] - Gamma[ 2n, z] )
> > Sum ? ? ---------------------------------------------------------
> > n=0 ? ? ? ? ? ? ?n**(2n+1) * n! * (2n+1)!

> If we take your denominator to be 2**(2*n+1)*n!*(n+1)!, your original
> integral would have been:

> Integral_k*z^z{I_1(t)*exp(-t)/t**2}*dt,

> right? ?Why should it be easier to evaluate the series rather
> than approximate the integral directly?

> --
> write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
> 6.0134700243160014d-154/),(/'x'/)); end

Yes - that looks about right, haven't double checked ...

In fact, that is what I am evaluating right now...I just thougth the
Gammafunction solution is more 'beautiful'...

Cheers - Ralf



Thu, 15 Sep 2011 07:40:24 GMT  
 
 [ 5 post ] 

 Relevant Pages 

1. complementary incomplete gamma function of complex argument

2. Incomplete Gamma function

3. Gamma and Bessel Functions in APL

4. The Gamma Function

5. Gamma function in PBCC

6. Gamma function in PBCC

7. Gamma function in formula node

8. gamma function, corrected...

9. gamma function...

10. Reciprocal gamma function of complex argument

11. Fastest code for Gamma function

12. Thanks to everybody for Gamma function algorithm

 

 
Powered by phpBB® Forum Software