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