Bessel function scaling + limited range?

[John W. Eaton]

On  4-Feb-2008, Kris Kuhlman wrote:

| I am trying to understand the scaling that Octave performs with besseli and
| besselk for real arguments.  I have worked with Amos' Bessel function
| routines in fortran, I understand that they return exp(-x)*besseli(n,x) and
| exp(x)*besselk(n,x) when you request a scaled Bessel function.  I have seen
| a couple previous bug reports in the Octave list regarding this topic, but
| they did not address my problem.
| 
| Why does Octave only return ~ 1.0E+33 as a maximum range?  When the argument
| is near 1 the exponential scaling has very little effect, so I can't compute
| very high orders of Bessel functions.  An example of my usage is below.
| Basically, for an argument I need a vector of Bessel functions of order 0 to
| ~35, Octave can only muster 0:23.
| 
| for e=1:ne
|   v(1,e) = sqrt(q)*exp(-eta(e));
|   v(2,e) = sqrt(q)*exp( eta(e));
|   besi(1:R+1,e) = besseli(0:R, v(1,e)).';
|   besk(1:R+1,e) = besselk(0:R, v(2,e)).';
| 
|   for n=1:N
|      Ke(n,e) = sum(AEv(1:R,n).* besi(1:R,e).*besk(1:R,e));
|      Ko(n,e) = sum(AOd(1:R,n).* ...
|            (besi(1:R,e).*besk(2:R+1,e) + besi(2:R+1,e).*besk(1:R,e)));
|   end
| end

This code snippet doesn't help much since you ahven't sent any of the
data to go along with it (waht are the values for ne, e, eta, q, R,
AEv, AOd, ...)?

| In the case I have been looking at just now, v(2,e) is ~0.6, so it goes like
| this:  The scaling option (the third argument) doesn't help me here.
| 
| GNU Octave, version 3.0.0
| -- snip --
| Octave was configured for "x86_64-unknown-linux-gnu".
| -- snip --
| octave:1> besselk(23,0.6)
| ans =  5.9453e+32
| octave:2> besselk(24,0.6)
| ans = Inf + Infi
| octave:3> besselk(23,0.6,1)
| ans =  1.0833e+33
| octave:4>besselk(24,0.6,1)
| ans = Inf + Infi
| 
| I guess I will not be able
| to have neutral code, I will have to write something Octave-specific into my
| code?

Or, you could help us fix the problem in Octave.  You say you have
experience using the Amos functions, so perhaps you can tell us how to
use them properly to solve your problem?

I'm attaching a test program in Fortran that calls ZBESK from Amos
directly.  The test program includes ZBESK and all of its
dependencies.  It calls ZBESK with KODE = 1 and then 2:

  C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
  C                    KODE= 1  RETURNS
  C                             CY(I)=K(FNU+I-1,Z), I=1,...,N
  C                        = 2  RETURNS
  C                             CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
  C

In both cases, for Z = 0.6 and FNU = 24, ZBESK returns IERR = 2, which
indicates that an overflow has occurred.  I don't know how to avoid
that problem.

Is there a newer version of the Amos library around somewhere that
fixes this problem?

jwe


-------------- next part --------------
A non-text attachment was scrubbed...
Name: bessel-test.f.gz
Type: application/octet-stream
Size: 54223 bytes
Desc: not available
Url : https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20080204/1d25fc48/attachment-0001.obj