[SciPy-user] Possible bug in bessel special function h2vp

[Robert Kern]

Julian Swartz wrote:
> Hi,
> 
> I've encountered what appears to be a bug in the function h2vp, which
> calculates the nth order derivative of a hankel function of the second kind.
> Here is an example with a second derivative being requested. The
> function returns two different values on successive calls. Any
> subsequent calls will yield the same result as the second call shown
> below. This problem is also evident in function h1vp, calculating the
> derivatives of a hankel function of the first kind.
> 
> 
>>>> from scipy import *
>>>> from scipy.special import *
>>>> n = 1.0
>>>> x = 2*0.3*pi
>>>> h2vp(n, x, 2)
> (2.9279764804171628e+165+0.22981683442664608j)
>>>> h2vp(n, x, 2)
> (-0.23515420288453626+0.54623524288750835j)
> 
> Even more interesting is that given that the hankel function of the
> second kind is defined as:
> 
>   (2)
> H    (x) = J  (x)  - jY  (x)
>   n            n          n
> 
> Where J and Y are the standard cylindrical bessel functions of the first
> and second kind. One would therefore expect that the second derivative
> of the Hankel function should be the sum of the second derivatives of
> the ordinary bessel functions. Using the scipy functions jvp and yvp,
> 
>>>> jvp(n, x, 2)
> -0.40831349982163045
>>>> yvp(n, x, 2)
> -0.18651604740953789
> 
> Which is quite different to either of the results above.
> 
> Any ideas on what the problem could be. (or am I just missing something?)

Not a clue. On OS X with the last CVS of scipy before the migration to
SVN, I don't see differences when repeating the call, but I do see a
difference from the stated definition for derivatives of order > 1:

In [1]: n = 1.0

In [2]: x = 2*0.3*pi

In [3]: special.h2vp(n,x,2)
Out[3]: (-0.26294530063194937+0.22981683442664608j)

In [4]: special.h2vp(n,x,2)
Out[4]: (-0.26294530063194937+0.22981683442664608j)

In [8]: special.jvp(n,x,2) - 1j*special.yvp(n,x,2)
Out[8]: (-0.40831349982163045+0.18651604740953792j)

In [23]: special.h2vp(n,x,1)
Out[23]: (-0.017916685502942564-0.58616758520739864j)

In [24]: special.jvp(n,x,1) - 1j*special.yvp(n,x,1)
Out[24]: (-0.017916685502942648-0.58616758520739876j)

In [21]: special.h2vp(n,x,3)
Out[21]: (-0.2382065153405013+0.02078698962754455j)

In [22]: special.jvp(n,x,3) - 1j*special.yvp(n,x,3)
Out[22]: (0.050806004596450884+0.10554382826770611j)

In [25]: special.hankel2(n,x)
Out[25]: (0.58147279675872476+0.1732031480684324j)

In [26]: special.jv(n,x) - 1j*special.yv(n,x)
Out[26]: (0.58147279675872487+0.17320314806843251j)

-- 
Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter