I am using boost spherical bessel functions and get an overflow error for the following call. double ytest; double x = 12.6; int n = 8; ytest = boost::math::sph_neumann(n, x); the value of x in not relevant as long as you're away from 0. it's the index that causes the problems, for n = 8, the call to cyl_neumann(n,x) causes the overflow. I tested this with other index values less than 8 and they all work fine. I have not tried greater than 8. MATLAB is well defined for these values and I can see no reason why this case would be special. Any ideas? Thanks, David
John,
Thanks for the reply. I did not give you enough details and apologize.
Here is a more detailed code snippet.
double a = 0.1;
double lambda = 0.05;
double pi = 3.1415926535897932384626433832795;
double kw = 2.0*pi/lambda;
int n = 8; // try 7, 9, etc.
cout << kw*a << endl;
double xtest = boost::math::sph_neumann(n,a*kw);
//double ytest = boost::math::cyl_neumann(n,a*kw);
cout << xtest << endl; //<< " " << ytest << endl;
Here is some more info to diagnose the problem.
The problem seems to be in the spherical neumann function sph_neumann, and it seems to be for n = 8.
When I posted the problem I quoted cyl_neumann( ) because ultimately that's the function that appeared in the debugger.
I inappropriately assumed that that function would cause problems, but the problem seems to be in sph_neumann.
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.
I am using VS2005 and boost version 1.0, (104400).
I hope this helps.
David
________________________________
From: John Maddock
Any ideas?
Which Boost version? This works for me as: std::cout << boost::math::cyl_neumann(8, 12.6) << std::endl; prints: 0.244094 John. _______________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
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.
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
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
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
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
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)?
I used: boost::math::cyl_neumann(8.5, 4 * pi); using your value of pi to reproduce the issue - I did say "exquisite" cancellation though - the values have to be just right - down to the last bit - to cause the issue. In fact the cancellation problem occurs after quite a bit of calculation so you have to be very unlucky indeed for things to work out "just so" :-(
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.
No it's an error in the logic inside bessel_jy in bessel_jy.hpp, I'm testing a fix now.
I have GSL on a linux machine at work and plan to try that too.
GSL doesn't support bessel J or Y with non-integer orders (unless they're changed since I last looked), they do support the spherical bessel functions but use a different calculation method so it probably works there. Still testing yours, John.
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)?
I used: boost::math::cyl_neumann(8.5, 4 * pi); using your value of pi to reproduce the issue - I did say "exquisite" cancellation though - the values have to be just right - down to the last bit - to cause the issue. In fact the cancellation problem occurs after quite a bit of calculation so you have to be very unlucky indeed for things to work out "just so" :-(
Ok, I could swear it worked for me but like you said, if the input is nudged a bit the issue goes away (or returns).
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.
No it's an error in the logic inside bessel_jy in bessel_jy.hpp, I'm testing a fix now.
I have GSL on a linux machine at work and plan to try that too.
GSL doesn't support bessel J or Y with non-integer orders (unless they're changed since I last looked), they do support the spherical bessel functions but use a different calculation method so it probably works there. Still testing yours, John. It sounds like you've tracked it down. Thanks. David. _____________________________________________ Boost-users mailing list Boost-users@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users
It sounds like you've tracked it down. Thanks. David.
Here's the patch to the headers: Index: bessel_jy.hpp =================================================================== --- bessel_jy.hpp (revision 77307) +++ bessel_jy.hpp (working copy) @@ -491,6 +485,16 @@ CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy T t = u / x - fu; // t = J'/J gamma = (p - t) / q; + // + // We can't allow gamma to cancel out to zero competely as it messes up + // the subsequent logic. So pretend that one bit didn't cancel out + // and set to a suitably small value. The only test case we've been able to + // find for this, is when v = 8.5 and x = 4*PI. + // + if(gamma == 0) + { + gamma = u * tools::epsilon<T>() / x; + } Ju = sign(current) * sqrt(W / (q + gamma * (p - t))); Jv = Ju * ratio; // normalization I'm re-running the full tests now before committing (there are changes to the tests not shown above), but the Bessel function tests all look good so far, John.
John Thanks, I added your patch and these functions no longer produce the errors I was seeing before. Please update me on your other testing. David
I'm re-running the full tests now before committing (there are changes to the tests not shown above), but the Bessel function tests all look good >so far,
John.
participants (2)
-
David Bergman -
John Maddock