Alex8
09-01-2006, 01:10 PM
Hello,
I use Fortran 90 IMSL Library routine DBCIS (dbl precision) to compute Modified Bessel Function of the 1st kind for any positive real order and complex argument. I know it is not directly from Num. Recipes, but the methods used in that routine are the same pretty much as those described in Num. Recipes.
I get a discontinuity in the Bessel function. I know that it should be continuous except when you cross branch cuts. But in my case, to my best knowledge, there are no branch cut crossings anywhere in sight.
Here is what I get:
Let
z - complex argument
nu - real order of the Bessel Function
I_nu(z) =is the value of the function
Here are the numbers:
nu=12.2998027613412
if z1=(24.9999744517597,-3.574092057101672E-002)
I_nu(z1) = (21.9616730993373,-3.578340339596549E-002)
if z2=(24.9999744517135,-3.574092057232411E-002)
I_nu(z2) = (19.4517591395845,-3.925690074051879E-002)
The only thing I know is that the second value is correct. I obtained it from integral representation of the Bessel function.
Z is very close to being real and there are no branch cuts in such a tight interval. As you can see, no matter how fine I make the subdivision of the interval [z1,z2], the Bessel function value still jumps.
I could not have made a programming mistake here as it is just a direct call of the IMSL routine. This very thing, though, precludes me from fixing anything as the routine is the black box.
Any ideas as to what is going on here would be greatly appreciated. If this is due to the problem with the DBCIS routine (rather than me being dumb; I hope this is the case), then maybe there are other routines that are free from this problem.
Thanks in advance.
I use Fortran 90 IMSL Library routine DBCIS (dbl precision) to compute Modified Bessel Function of the 1st kind for any positive real order and complex argument. I know it is not directly from Num. Recipes, but the methods used in that routine are the same pretty much as those described in Num. Recipes.
I get a discontinuity in the Bessel function. I know that it should be continuous except when you cross branch cuts. But in my case, to my best knowledge, there are no branch cut crossings anywhere in sight.
Here is what I get:
Let
z - complex argument
nu - real order of the Bessel Function
I_nu(z) =is the value of the function
Here are the numbers:
nu=12.2998027613412
if z1=(24.9999744517597,-3.574092057101672E-002)
I_nu(z1) = (21.9616730993373,-3.578340339596549E-002)
if z2=(24.9999744517135,-3.574092057232411E-002)
I_nu(z2) = (19.4517591395845,-3.925690074051879E-002)
The only thing I know is that the second value is correct. I obtained it from integral representation of the Bessel function.
Z is very close to being real and there are no branch cuts in such a tight interval. As you can see, no matter how fine I make the subdivision of the interval [z1,z2], the Bessel function value still jumps.
I could not have made a programming mistake here as it is just a direct call of the IMSL routine. This very thing, though, precludes me from fixing anything as the routine is the black box.
Any ideas as to what is going on here would be greatly appreciated. If this is due to the problem with the DBCIS routine (rather than me being dumb; I hope this is the case), then maybe there are other routines that are free from this problem.
Thanks in advance.