id,summary,reporter,owner,description,type,status,milestone,component,version,severity,resolution,keywords,cc
2733,Cylindrical bessel fuction gives wrong result if n=-1 and z>1,Ivan Lisenkov <ivan@…>,John Maddock,"If the order of Bessel function is -1 and argument is more than 1.0 the algorithm goes on wrong branch, giving always 0 as result. Here a program to illustrate the problem:


{{{
#include <iostream>
#include <boost/math/special_functions/bessel.hpp>

int main()
{
	std::cout.setf (std::ios::scientific | std::ios::left);
	std::cout.precision(5);
		
	double z=0.99999;	
	std::cout<<""z=""<<z<<"" J_{-1}(z)=""<<boost::math::cyl_bessel_j(-1, z)<<'\n';
	z=1.00007e+00;
	std::cout<<""z=""<<z<<"" J_{-1}(z)=""<<boost::math::cyl_bessel_j(-1, z)<<'\n';

	return 0;
}
}}}

The output of the program is:
{{{
z=9.99990e-01 J_{-1}(z)=-4.40047e-01
z=1.00007e+00 J_{-1}(z)=-0.00000e+00
}}}

I have created a patch against the latest subverversion:
{{{
Index: boost/math/special_functions/detail/bessel_jn.hpp
===================================================================
--- boost/math/special_functions/detail/bessel_jn.hpp   (revision 51047)
+++ boost/math/special_functions/detail/bessel_jn.hpp   (working copy)
@@ -32,9 +32,12 @@
     {
         return bessel_j0(x);
     }
-    if (n == 1)
+    if ((n == 1) || (n == -1))
     {
-        return bessel_j1(x);
+        factor = n;  // J_{-1}(z) = -1*J_1(z)
+        value = bessel_j1(x);
+        value *=factor;
+        return value;
     }
     if (n < 0)
     {
}}}

The output after the patch is:

{{{
z=9.99990e-01 J_{-1}(z)=-4.40047e-01
z=1.00007e+00 J_{-1}(z)=-4.40073e-01
}}}

This is correspond to GNU Octave:
{{{
octave:32> besselj(-1,0.99999)
ans = -0.44005 - 0.00000i

octave:26> besselj(-1,1.00007)
ans = -0.44007 - 0.00000i
}}}",Patches,closed,Boost 1.38.0,math,Boost 1.37.0,Problem,fixed,bessel,