Kris Kuhlman

2008-02-04 17:32:00 UTC

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

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/

Kristopher L. Kuhlman

PhD candidate Department of Hydrology & Water Resources

University of Arizona

(520) 621-1380

http://www.u.arizona.edu/~kuhlman/