From bug-octave-request at bevo dot che dot wisc dot edu Fri Jun 8 02:50:16 2001 Subject: Re: besseli output range apparently limited to about 10^33 From: Daniel Alt To: "John W. Eaton" CC: bug-octave at bevo dot che dot wisc dot edu Date: Fri, 8 Jun 2001 02:49:30 -0500 John, Thanks for the pointer! I hadn't noticed the scaling feature of Octave's Bessel functions (I don't think Matlab has that). Besides allowing a workaround for my problem, it will let me use Bessel functions with values exceeding the floating point range. The workaround I have decided to use may work as a real fix, though it's not super-efficient. First I calculate my result using the unscaled Bessel function. If the result is not finite, I try again with the scaled Bessel function (OPT=1) and use that answer multiplied by exp(X). In my application the second argument (X) to the Bessel function is always real and positive, so I don't bother checking the sign of X. (I was planning to check if the result is zero, and conditionally redo the calculation using the Bessel function's other scaling (OPT=2), but the value of OPT doesn't seem to make any difference: see the next paragraph). After playing around a little with the Bessel function in Octave it appears that the correct thing to do for besseli with purely real arguments is multiply besseli(ALPHA,X,1) by exp(abs(X)) to get besseli(ALPHA,X). Since besseli(ALPHA,X) is, I think, defined to be i^(-ALPHA)*besselj(ALPHA,i*X), I initially only looked at the besselj function for complex arguments. It appears that, contrary to the information provided by "help", besselj(ALPHA,X,OPT) is scaled by exp(-abs(imag(X))) regardless of the value of OPT. For besseli the scaling uses real(X) instead of imag(X). Here's a transcript showing this: octave:1> besseli(19,20) ans = 7659.0 octave:2> besseli(19,20,1) ans = 1.5786e-05 octave:3> besseli(19,20,2) ans = 1.5786e-05 octave:4> besseli(19,20,-1) ans = 1.5786e-05 octave:5> besseli(19,20,0) ans = 1.5786e-05 octave:6> besseli(19,20,Inf) ans = 1.5786e-05 octave:7> besseli(19,20,-Inf) ans = 1.5786e-05 octave:8> besseli(19,20,i) ans = 1.5786e-05 octave:9> besseli(19,20,-i) ans = 1.5786e-05 octave:10> besseli(19,20,-i)*exp(20) ans = 7659.0 octave:11> besselj(19,20+5i) ans = 0.29279 + 0.97593i octave:12> besselj(19,20+5i,1) ans = 0.0019728 + 0.0065758i octave:13> besselj(19,20+5i,2) ans = 0.0019728 + 0.0065758i octave:14> besselj(19,20+5i,0) ans = 0.0019728 + 0.0065758i octave:15> besselj(19,20+5i,Inf) ans = 0.0019728 + 0.0065758i octave:16> besselj(19,20+5i,-Inf) ans = 0.0019728 + 0.0065758i octave:17> besselj(19,20+5i,-Inf)*exp(5) ans = 0.29279 + 0.97593i octave:18> besselj(19,20-5i,-Inf)*exp(5) ans = 0.29279 - 0.97593i octave:19> besselj(19,20-5i) ans = 0.29279 - 0.97593i octave:20> besselj(20,19-6i) ans = -0.69458 - 0.27057i octave:21> besselj(20,19-6i,1) ans = -1.7217e-03 - 6.7067e-04i octave:22> besselj(20,19-6i,1)*exp(6) ans = -0.69458 - 0.27057i octave:23> besseli(5,7+18i) ans = 53.231 - 59.319i octave:24> besseli(5,7+18i,1)*exp(7) ans = 53.231 - 59.319i octave:25> besseli(5,-13+18i) ans = -1.7694e+04 - 2.0504e+04i octave:26> besseli(5,-13+18i,1)*exp(13) ans = -1.7694e+04 - 2.0504e+04i Cheers, Dan Alt "John W. Eaton" wrote: > > On 4-Jun-2001, Daniel Alt wrote: > > | "Bug" report for Octave 2.1.31 configured for MS Windows > | (octave-windows-2000oct25d.exe downloaded 1 June 2001 from > | sourceforge.net) > | > | (I attempted, unsuccessfully, to use Octave's bug_report command to > | generate this email. My PC doesn't have Unix shell command emulators > | installed, so bug_report fails with "not found" errors for grep, cat, > | chmod, etc.) > | > | It might be stretching the definition to call this a bug, but besseli's > | output goes to infinity at around 10^33 (instead of at the maximum > | floating point value of about 10^308). Presumably this is because some > | intermediate result exceeds the floating point range, but both Matlab > | and Microsoft Excel can calculate besseli values of at least 10^300. > | > | Here is a bit of an Octave transcript showing the transition: > | > | octave:87> besseli(19,81.3565) > | ans = 1.0327e+33 > | octave:88> besseli(19,81.357) > | ans = Inf + Infi > | > | As one might expect, besselj has similar behaviour for equivalent > | inputs: > | > | octave:92> besselj(19,81.357i) > | ans = Inf + Infi > | octave:93> besselj(19,81.3565i) > | ans = -6.3238e+16 - 1.0327e+33i > | > | (One question on this besselj result: Shouldn't ans be purely imaginary > | in the second case?) > > Octave uses the subroutines provided by the "Amos" library from Netlib > to compute Bessel functions. > > You can ask for the results to be scaled to avoid overflow. For > example, I noticed that > > besseli (19, 700) > > fails, but > > exp (700) * besseli (19, 700, 1) > > produces the same result as Matlab does for besseli (19, 700). > > I'm not sure why the Fortran code doesn't get this right. > > I'm not sure that it would be correct for Octave to simply call the > Fortran code and request the scaled result, then multiply that by the > inverse scale factor before returning. It seems that it might not > always be as accurate. > > I'd appreciate it if anyone could shed some light on this and maybe > provide a fix. > > Thanks, > > jwe ------------------------------------------------------------- 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 -------------------------------------------------------------