Fwd: Bessel function scaling + limited range?

[Kris Kuhlman]

I have used a modified fortran90 version of these libraries, as distributed
by Allan Miller;  TOMS 644.  The fortran code was downloaded from the
address below.

http://users.bigpond.net.au/amiller/toms.html

Kris


On Feb 4, 2008 2:08 PM, John W. Eaton <jwe at bevo.che.wisc.edu> wrote:

> 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
>
>
>


-- 

Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/ <http://www.u.arizona.edu/%7Ekuhlman/>



-- 

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/ee91d0d4/attachment.html