I wanted to share a couple more pieces of info. 

1.  I tried your pseudo code below

Y[8.5](4*Pi)

using my value of pi (I believe that it was in an earlier post) and did not get a bad answer.  so if this is really a cancellation error could it be related to the precision of pi (how many decimals did you use)?

2.  Since we are calling bessel of a half integer order I would suspect that the Gamma( ) function is called.  Could a bad value from that function caused the error?  If so I might expect the debug statement to trace back that far but they didn't.

3.  Matlab gave a good answer for the same (similar) input.

I have GSL on a linux machine at work and plan to try that too. 

Thanks for looking into this.

David






John,

Thanks, I am glad that you can reproduce the problem.  Based on the debug messages I would guess that your sph_neumann wraps your cyl_neumann, this is typically how its done.  I tried calling cyl_neumann with nu = 8.5 and the same value of x (ref Arfkin).

double ytest = boost::math::cyl_neumann(double(n)+0.5,a*kw);

This works too.  I'm not sure if my type casting matches what boost does.
So, I can wrap my own and they should work.
Based on what I expect of the bessel function I do not believe that there should be any ill-behaved points for non-zero values of x.
When I got the overflow error I initially thought it was do to a recursive assignment as it was wrapped in a template function of my own and called from within a loop (not numerical instability).  But when I could reproduce the problem outside the loop and the wrapper I figured I'd email the user group.

Please let me know what you find.

David


From: John Maddock <boost.regex@virgin.net>
To: boost-users@lists.boost.org
Sent: Friday, March 23, 2012 6:10 AM
Subject: Re: [Boost-users] spherical bessel functions

> It doesn't matter if I define double x = a*kw and pass this into the interface, what seems to matter is the value of the index,
> when n = 8 it doesn't work. It works for 9 and 10 and all integers before 8.
> I can replace x with a, lambda or pi and it works for n = 8, so it seems to be this combination of numbers or something I'm doing
> wrong w/r to type casting variables etc.

OK I can reproduce here, and it's very specific to the calculation of Y[8.5](4*Pi), if you change nu at all, or add a tiny increment to the 4PI arg then it all works OK.  The problem BTW is exquisite cancellation inside the algorithm, which leads me to suspect there's something "special" about this value.  Unfortunately I don't have a fix at present (or understand the root cause yet), more investigation needed....

Thanks for the report, John.
_______________________________________________
Boost-users mailing list
Boost-users@lists.boost.org
http://lists.boost.org/mailman/listinfo.cgi/boost-users