From bug-octave-request at bevo dot che dot wisc dot edu Thu Aug 28 14:24:31 2003 Subject: small argument besselj fails From: "John W. Eaton" To: Fabian at isas-berlin dot de Cc: bug-octave at bevo dot che dot wisc dot edu Date: Thu, 28 Aug 2003 14:24:13 -0500 On 4-Jul-2003, Fabian at isas-berlin dot de wrote: | For small arguments x, Octave's besselj ( nu, x ) | drops suddenly to (complex) zero. | | This contradicts theory. | Zero should be approached as | | ( 1/2*x ).^nu / gamma( 1 + nu ) | | See Abramowitz & Stegun ( Theorem 3.1.7 p 360 ). | | You may run accomp. script for a demo. of the problem. | | Rolf | | | | | | | | % Octave's result for besselj( nu, x ) , x->0 | % doesn't follow the theoretical small angle limit. | % | % Checked with Octave V2.1.49 ( Win NT 5.0 ) | % Rolf Fabian 030704 | | | 1; | | nu = 30; % order of Besselj | x = linspace(0,5,500).'; | | t = (1/2*x).^nu/gamma(nu+1); % Theory for small x->0 | % see Abramowitz & Stegun | | y = besselj(nu,x); | | xlabel('x'); ylabel('y'); | semilogy(x,real(y),';Octave;',x,t,';theor. limit x -> 0;') | xlabel(''); ylabel(''); The following patch should fix the complex problem. The jump to zero appears to be a bug in the zbesj function from the amos library. I don't understand that code well enough to fix it easily. A solution would be greatly appreciated. The following fortran program might make it easier to debug than it would be using Octave. You can get the code for zbesj from www.netlib.org in the amos directory. Octave's version of that code is currently not modified from the original. jwe program foo double precision zr, zi, fnu, cyr, cyi integer i, kode, n, nz, ierr zi = 0.0d0 fnu = 30.0d0 kode = 1 n = 1 do 10 i = 0, 500 zr = i * .01d0 call zbesj (zr, zi, fnu, kode, n, cyr, cyi, nz, ierr) if (zr .ge. 0.0d0 .and. zi .eq. 0.0d0) cyi = 0.0d0 print *, zr, cyr, cyi, nz, ierr 10 continue end 2003-08-28 John W. Eaton * lo-specfun.cc (zbesj, zbesy, zbesi, zbesk, airy, biry): Also zero imaginary part of result if real part of input value is zero. Index: lo-specfun.cc =================================================================== RCS file: /usr/local/cvsroot/octave/liboctave/lo-specfun.cc,v retrieving revision 1.13 retrieving revision 1.14 diff -u -r1.13 -r1.14 --- lo-specfun.cc 15 Nov 2002 04:47:02 -0000 1.13 +++ lo-specfun.cc 28 Aug 2003 19:03:06 -0000 1.14 at @ -219,7 +219,7 @@ F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr); - if (zi == 0.0 && zr > 0.0) + if (zi == 0.0 && zr >= 0.0) yi = 0.0; retval = bessel_return_value (Complex (yr, yi), ierr); at @ -272,7 +272,7 @@ F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz, &wr, &wi, ierr); - if (zi == 0.0 && zr > 0.0) + if (zi == 0.0 && zr >= 0.0) yi = 0.0; } at @ -314,7 +314,7 @@ F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr); - if (zi == 0.0 && zr > 0.0) + if (zi == 0.0 && zr >= 0.0) yi = 0.0; retval = bessel_return_value (Complex (yr, yi), ierr); at @ -365,7 +365,7 @@ { F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr); - if (zi == 0.0 && zr > 0.0) + if (zi == 0.0 && zr >= 0.0) yi = 0.0; } at @ -616,7 +616,7 @@ F77_FUNC (zairy, ZAIRY) (zr, zi, id, kode, ar, ai, nz, ierr); - if (zi == 0.0 && (! scaled || zr > 0.0)) + if (zi == 0.0 && (! scaled || zr >= 0.0)) ai = 0.0; return bessel_return_value (Complex (ar, ai), ierr); at @ -637,7 +637,7 @@ F77_FUNC (zbiry, ZBIRY) (zr, zi, id, kode, ar, ai, ierr); - if (zi == 0.0 && (! scaled || zr > 0.0)) + if (zi == 0.0 && (! scaled || zr >= 0.0)) ai = 0.0; return bessel_return_value (Complex (ar, ai), ierr); ------------------------------------------------------------- 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 -------------------------------------------------------------