Re: strange besselj behavior

[David Bateman]

Ahn Kyung wrote:

> Folks,
>
> Test this.
>
> x=linspace(0,100,10000)';
> plot(x, besselj(0,x));
> plot(x, besselj(2,x));
>
> When I tried this, besselj(0,x) did not have any problem, but
> plot(x,besselj(2,x)) spit out
>
> octave:4> plot(x,besselj(2,x))
>
> gnuplot> pl '/tmp/oct-PJiXlG' t "line 1"
> ^
> line 0: Bad data on line 5
>
> So I checked what besselj(2,x) looked like. Here's the comparison with
> besselj(0,x):
>
> octave:5> besselj(0,x)
> ans =
>
> 1.0000e+00
> 9.9997e-01
> 9.9990e-01
> 9.9977e-01
> 9.9960e-01
> 9.9937e-01
> 9.9910e-01
> ...
>
> octave:6> besselj(2,x)
> ans =
>
> 0.00000 - 0.00000i
> 0.00001 + 0.00000i
> 0.00005 + 0.00000i
> 0.00011 + 0.00000i
> 0.00020 + 0.00000i
> 0.00031 + 0.00000i
> 0.00045 + 0.00000i
>
> ---------------
>
> So, besselj in some cases(besselj(3,x)->bad; besselj(4,x)->OK, etc.)
> thinks that it can have imaginary values like Hankel functions?? The
> plot looked OK when I tried plot(x,real(besselj(2,x))), so at least I
> may trust the result with real(..).
>
> I tried this just to see how well besselj would work. I'm using
> octave-2.1.50, installed via rpm on redhat core 2 platform. Well, I
> also have a Mandrake 10 box with octave installed with mandrake
> specific rpm and guess what. The same problem.
>
> If any of you experience no such problems in other versions or other
> platforms, it may be just some bad compilation. Otherwise, I'd go
> ahead and report it as a bug.
>
> So I'd appreciate if you could share your experience on this.
>
>
>

Checking against 2.1.64, I don't have this problem.

octave:1> x=linspace(0,100,10000)';
octave:2> besselj(0,x(1:6))
ans =

1.00000
0.99997
0.99990
0.99977
0.99960
0.99937

octave:3> besselj(2,x(1:6))
ans =

0.00000000
0.00001250
0.00005001
0.00011251
0.00020001
0.00031250


Then checking the liboctave changelog I find this

* lo-specfun.cc (is_integer_value): New function.
(zbesj, zbesi, zbesy): Special case negative integer or half
integer orders that cause overflow for small arguments.


dated 2004-07-22.I can't check as I no longer have a version of 2.1.57
or earlier that should be the last version to demonstrate this.However,
reading the changelog message, if x started at a larger value than zero
the problem probably would not have manifested... In any case I believe
the problem is fixed.

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------