Hello everyone,

I try to plot the digamma function of (1/2 + 1/x) but I’m not sure that I’m plotting the good one.

I’ve tried:

special.polygamma(0, (1/2 + 1/x))

and

special.polygamma(1, (1/2 + 1/x))

but I don’t have the same result as with mathcad.

I’ve tried to code it like that:

def F(x): return mpmath.diff(lambda x: gamma(1/2 + 1/x),1)/gamma(1/2 + 1/x)

But It returns zero division error even when x is in ]0,1]

Any idea?

Hello everyone,

I try to plot the digamma function of (1/2 + 1/x) but I'm not sure that I'm

plotting the good one.

I've tried:

special.polygamma(0, (1/2 + 1/x))

and

special.polygamma(1, (1/2 + 1/x))

You want special.polygamma(0, (1/2 + 1/x)). See

http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.polygamma.html

The number specifies which derivative of the digamma function you want.

Surely you want the 0th derivative?

But It returns zero division error even when x is in ]0,1]

I think it blows up at x = 0. What is the type of x in your usecase? Is

it an array? If x contains the element 0, you will get a zero

division error. You could try plotting the points explicitly:

from numpy import linspace

from pylab import *

x = linspace(0.5, 2, num=100, endpoint=True)

y = special.polygamma(0, (1/2 + 1/x))

plot(x, y)

show()

You can compare output against this:

http://www.wolframalpha.com/input/?i=digamma(1%2F2+%2B+1%2Fx)+between+0.5+and+2

Hope this helps.

## ···

On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:

--

Damon McDougall

http://damon-is-a-geek.com

B2.39

Mathematics Institute

University of Warwick

Coventry

West Midlands

CV4 7AL

United Kingdom

Another problem might be the “1/2” part, which in python2.x would yield 0 unless one does “from **future** import division”.

Ben Root

## ···

On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall <damon.mcdougall@…287…> wrote:

On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:

Hello everyone,

I try to plot the digamma function of (1/2 + 1/x) but I’m not sure that I’m

plotting the good one.

I’ve tried:

special.polygamma(0, (1/2 + 1/x))

and

special.polygamma(1, (1/2 + 1/x))

You want special.polygamma(0, (1/2 + 1/x)). See

http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.polygamma.html

The number specifies which derivative of the digamma function you want.

Surely you want the 0th derivative?

But It returns zero division error even when x is in ]0,1]

I think it blows up at x = 0. What is the type of x in your usecase? Is

it an array? If x contains the element 0, you will get a zero

division error. You could try plotting the points explicitly:

from numpy import linspace

from pylab import *

x = linspace(0.5, 2, num=100, endpoint=True)

y = special.polygamma(0, (1/2 + 1/x))

plot(x, y)

show()

You can compare output against this:

http://www.wolframalpha.com/input/?i=digamma%281%2F2+%2B+1%2Fx%29+between+0.5+and+2

Hope this helps.

Wow, I can't believe I didn't spot that. Nice one.

I will update my answer according to Ben's astute observation:

from scipy import special

from pylab import *

x = linspace(0.5, 2.0, num=100, endpoint=True)

y = special.polygamma(0, 0.5 + 1.0/x)

plot(x, y)

show()

Thanks Ben.

## ···

On Tue, Jul 10, 2012 at 08:57:24AM -0400, Benjamin Root wrote:

On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall > <damon.mcdougall@...287...>wrote:

> On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:

>

> > But It returns zero division error even when x is in ]0,1]

>

> I think it blows up at x = 0. What is the type of x in your usecase? Is

> it an array? If x contains the element 0, you will get a zero

> division error. You could try plotting the points explicitly:

>

Another problem might be the "1/2" part, which in python2.x would yield 0

unless one does "from __future__ import division".

Ben Root

--

Damon McDougall

http://damon-is-a-geek.com

B2.39

Mathematics Institute

University of Warwick

Coventry

West Midlands

CV4 7AL

United Kingdom

Thanks! The problem came from the 1/2 ! For convenience I’ve found the function digamma on numpy

**http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.psi.html**

But it’s quite hard to find it! Maybe we can ask to add digamma in the title between parenthesis?

Fabien

2012/7/10 Damon McDougall <damon.mcdougall@…287…>

## ···

On Tue, Jul 10, 2012 at 08:57:24AM -0400, Benjamin Root wrote:

On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall > > > <damon.mcdougall@…1003…7…>wrote:

On Tue, Jul 10, 2012 at 12:27:59PM +0200, Fabien Lafont wrote:

But It returns zero division error even when x is in ]0,1]

I think it blows up at x = 0. What is the type of x in your usecase? Is

it an array? If x contains the element 0, you will get a zero

division error. You could try plotting the points explicitly:

Another problem might be the “1/2” part, which in python2.x would yield 0

unless one does “from **future** import division”.

Ben Root

Wow, I can’t believe I didn’t spot that. Nice one.

I will update my answer according to Ben’s astute observation:

from scipy import special

from pylab import *

x = linspace(0.5, 2.0, num=100, endpoint=True)

y = special.polygamma(0, 0.5 + 1.0/x)

plot(x, y)

show()

Thanks Ben.

–

Damon McDougall

http://damon-is-a-geek.com

B2.39

Mathematics Institute

University of Warwick

Coventry

West Midlands

CV4 7AL

United Kingdom