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
> !------------------------------------------------------------------------------------------------

--

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)!

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)!

> 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

 Page 1 of 1 [ 5 post ]

Relevant Pages