From: Fabien Lafont <lafont.fabien@gm...>  20120710 10:28:05
Attachments:
Message as HTML

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? 
From: Damon McDougall <damon.mcdougall@gm...>  20120710 11:05:14

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.  Damon McDougall http://damonisageek.com B2.39 Mathematics Institute University of Warwick Coventry West Midlands CV4 7AL United Kingdom 
From: Benjamin Root <ben.root@ou...>  20120710 12:58:02
Attachments:
Message as HTML

On Tue, Jul 10, 2012 at 7:05 AM, Damon McDougall <damon.mcdougall@...>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. > > Another problem might be the "1/2" part, which in python2.x would yield 0 unless one does "from __future__ import division". Ben Root 
From: Damon McDougall <damon.mcdougall@gm...>  20120710 13:02:40

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@...>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://damonisageek.com B2.39 Mathematics Institute University of Warwick Coventry West Midlands CV4 7AL United Kingdom 
From: Fabien Lafont <lafont.fabien@gm...>  20120710 13:21:49
Attachments:
Message as HTML

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@...> > 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@...>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://damonisageek.com > B2.39 > Mathematics Institute > University of Warwick > Coventry > West Midlands > CV4 7AL > United Kingdom > 