Bessel function zeros in Fortran?

[Nasser M. Abbasi]

What do folks use to find Bessel function zeros in fortran?

http://mathworld.wolfram.com/BesselFunctionZeros.html

I see that gfortran has Intrinsic for Bessel functions,
but not for finding the zeros

https://gcc.gnu.org/onlinedocs/gfortran/Intrinsic-Procedures.html#Intrinsic-Procedures

"BESSEL_J0: Bessel function of the first kind of order 0
BESSEL_J1: Bessel function of the first kind of order 1
BESSEL_JN: Bessel function of the first kind
BESSEL_Y0: Bessel function of the second kind of order 0
BESSEL_Y1: Bessel function of the second kind of order 1
BESSEL_YN: Bessel function of the second kind"

The above is really nice that Fortran has these. But what about
for the zeros?

I see some people have written Fortran functions for this, such
as this page

http://jean-pierre.moreau.pagesperso-orange.fr/Fortran/mjyzo_f90.txt

"SUBROUTINE JYZO(N,NT,RJ0,RJ1,RY0,RY1)
! Purpose: Compute the zeros of Bessel functions Jn(x),
! Yn(x), and their derivatives
"

But wanted to ask first, if there is another way to do
this in Fortran, or another standard library to use, other
than having to download code from the internet and use it.
I searched http://www.netlib.org also, and see many
functions there for Bessel but did not see one for its zeros
so far. Looked at Lapack, but saw nothing there either.

I think Bessel functions zeros should be an Intrinsic. This
is used when solving PDE's such a wave and heat PDE,
and easier to use than one having to do root finding
and change the interval each time.

There is also spherical bessel functions zeros which is useful
to have as well.

I think Fortran should spend its efforts in adding more
such useful functions such as these. These functions and
much more, come build-in in Mathematica, Maple and may be
Matlab as well, out of the box, which makes it easy to use.

http://reference.wolfram.com/language/ref/BesselJZero.html

--Nasser

[steve kargl]

Nasser M. Abbasi wrote:

> What do folks use to find Bessel function zeros in fortran?

Bisection with a good guess for the value of zero. The guess
is trivial.

> The above is really nice that Fortran has these. But what about
> for the zeros?

Other than x = 0 for J_n(x) for n > 0, the zeros of the Bessel and
Neumann functions are not exactly representable. At best, you
can try to locate the closest floating point argument that
approximates a zero.


--
steve

[spectrum]

Though not sure if this is a popular approach, I often look at the implementation
in Scipy and find an underlying Fortran code, for example,

1) Find zeros + bessel things.
https://docs.scipy.org/doc/scipy/reference/special.html#zeros-of-bessel-functions

2) Choose the function of closest interest and click the [source] button on the right, e.g.
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jnjnp_zeros.html#scipy.special.jnjnp_zeros

3) Check the code and get the literature or underlying library info.
https://github.com/scipy/scipy/blob/v1.1.0/scipy/special/basic.py#L158-L196

which seems to trace back to the code by "specfun.f".
https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f

A particular routine "jdzo()" is here:
https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f#L380

It looks like f2py is used to wrap such a function:
https://github.com/scipy/scipy/blob/master/scipy/special/specfun.pyf#L36

The original code may be here (scroll down to find Fortran codes):
https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.f90

-- Comparing the code for "jdzo()", the Scipy version still uses ".eq." etc while
the above special_functions.f90 uses "==". On the other hand, the former uses
"double precision" while the latter uses literal "kind = 8". So the former seems
to be modified by Scipy for their use.

-- It also seems like some indexing of local arrays are different (e.g. Z0(1400) vs
Z0(0:1400), which might be due to 0- vs 1-based indexing thing (not sure at all).
So for pure Fortran, the original version might be more straightforward to test or use
(not sure again...). Anyway, I think we can validate the code by comparing
with some reference results (available from other sites).

-----
Other things found from the Internet search:

Zeros of Regular Bessel Functions (GSL)
https://www.gnu.org/software/gsl/manual/html_node/Zeros-of-Regular-Bessel-Functions.html

Bessel function zeros (in mpmath for arbitrary precision?)
http://docs.sympy.org/0.7.3/modules/mpmath/functions/index.html
http://docs.sympy.org/0.7.3/modules/mpmath/functions/bessel.html#bessel-function-zeros
http://mpmath.org/doc/current/functions/bessel.html#bessel-function-zeros

Boost: Finding Zeros of Bessel Functions of the First and Second Kinds
https://www.boost.org/doc/libs/1_67_0/libs/math/doc/html/math_toolkit/bessel/bessel_root.html

[spectrum]

Hmmm, the Scipy version uses floating-literal constants without "D0" etc at several
places...

https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f#L380

while the original version attaches "D0" for the same piece of codes (e.g., please
look at jdzo()). So this latter version might be more reliable
(unless the single-precision literals have virtually no effect on the final results,
or some compiler options (??) are used to promote such literals by Scipy).

https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.f90

[Clive Page]

On 19/05/2018 03:40, Nasser M. Abbasi wrote:
> I think Fortran should spend its efforts in adding more
> such useful functions such as these. These functions and
> much more, come build-in in Mathematica, Maple and may be
> Matlab as well, out of the box, which makes it easy to use.

I'm not convinced of this. Bessel functions are already in the Standard, of course.

The difficulty is that Fortran is used in a very wide variety of fields, and each of us has a pet collection of functions that we use regularly and therefore think should be built in, but there is generally a different set for different programmers. But if we were to add all of them: well firstly there is not always a really good universally agreed algorithm for them (take the arguments on how to generate pseudo-random numbers as an example of the difficulties here), and secondly the result would be a language standard hugely expanded and dominated by the definition of all the new intrinsic procedures. I think these are best left in external modules, personally.


--
Clive Page

[Ron Shepard]

The reason that all of these other functions are included in these other
languages is that they are interpreted, so it is difficult, sometimes
almost impossible, to write efficient low-level code in a native way. In
Matlab, for example, the built-in functions are often implemented in
fortran. That feature does not apply to fortran, which is compiled not
interpreted, and which allows efficient low-level code to be written
natively.

Derivatives or zeros of Bessel functions seems like kind of an obscure
request. But if that is what you need, and you needed the code to be
written yesterday, then I can see how digging through the literature,
exploring various algorithms, and writing the code yourself might be an
obstacle. On the other hand, once this is done in fortran, your
implementation will likely be just as efficient as if it were written in
assembler, something that you would not expect if you were writing the
code in an interpreted language. And also, as the OP noted, there is a
vast library of such codes written in fortran over the last 50+ years
that can be used either directly or after modifications.

$.02 -Ron Shepard

[Charlie Roberts]

On Fri, 18 May 2018 21:40:06 -0500, "Nasser M. Abbasi" <n...@12000.org>
wrote:

>What do folks use to find Bessel function zeros in fortran?
>
>http://mathworld.wolfram.com/BesselFunctionZeros.html

It depends on what you want to pay! Over many years,
as a user of the NAG Fortran Library I had ready access
to what you need.

https://www.nag.co.uk/numeric/CL/nagdoc_cl23/pdf/INDEXES/KWIC/bessel.html

S17ALC is what you need .... and I am sure that IMSL
has a routine, too. But I think that you are
looking for something 'free'. Unfortunately, Netlib
does not appear to have what you need.

Not sure if Numerical Recipes has a zero finder .... I
do not have my copy handy. Their code is well thought
out in most cases.

If you want to do some coding, you may want to read

https://www.sciencedirect.com/science/article/pii/037704279190169K

but, having spent a lifetime in numerically intensive work
I would not advice writing such code unless you have experience
in producing stable, numerically accurate, fast code.

[Themos Tsikas]

Just to bring this up to date and to give a Fortran routine

https://www.nag.co.uk/numeric/fl/nagdoc_latest/html/s/s17alf.html


Best Regards
Themos Tsikas

On Monday, 21 May 2018 14:27:10 UTC+1, Charlie Roberts wrote:
> On Fri, 18 May 2018 21:40:06 -0500, "Nasser M. Abbasi"

[shmar...@gmail.com]

Of course it would depend on your specific application, but do you really need to generate the zeros on the fly? Can you not generate a table of zeros for the various orders you'll ever need out to some maximum x; and then just load the table every time you run your program?

[James Van Buskirk]

"Themos Tsikas" wrote in message
news:512ea752-b802-4cdb...@googlegroups.com...

> Just to bring this up to date and to give a Fortran routine

> https://www.nag.co.uk/numeric/fl/nagdoc_latest/html/s/s17alf.html

So this subroutine doesn't handle Robin boundary conditions?

[Themos Tsikas]

Hello,

I am not sure what you mean, in this context. It doesn't handle ANY boundary conditions, it does not solve a general ODE or PDE.

Themos Tsikas

[William Clodius]

Nasser M. Abbasi <n...@12000.org> wrote:

> What do folks use to find Bessel function zeros in fortran?
>
> http://mathworld.wolfram.com/BesselFunctionZeros.html
>
> I see that gfortran has Intrinsic for Bessel functions, but not for
> finding the zeros
>
> https://gcc.gnu.org/onlinedocs/gfortran/Intrinsic-Procedures.html#Intrinsi
> c-Procedures
>
> "BESSEL_J0: Bessel function of the first kind of order 0
> BESSEL_J1: Bessel function of the first kind of order 1
> BESSEL_JN: Bessel function of the first kind
> BESSEL_Y0: Bessel function of the second kind of order 0
> BESSEL_Y1: Bessel function of the second kind of order 1
> BESSEL_YN: Bessel function of the second kind"
>
> The above is really nice that Fortran has these. But what about
> for the zeros?

> <snip>
For Bessel functions of the first kind with integer orders, you can use
a polynomial root finder. Note that you rapidly loose precision as the
order increases and the roots of the numerical function can differ
significantly from the roots of the ideal function.