Bessel function scaling + limited range?

Kris Kuhlman kuhlman at hwr.arizona.edu
Mon Feb 4 11:32:00 CST 2008


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

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 would appreciate any pointers or advice as to how I can move my matlab
code over to octave (this works fine in Matlab), I guess I will not be able
to have neutral code, I will have to write something Octave-specific into my
code?

Thanks,

Kris

-- 

Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20080204/228f0e20/attachment-0001.html 


More information about the Bug-octave mailing list