`besselj_zero` failing to return some zeros (v1.1.0)

[ChristopheBerthod]

In v1.1.0, the function besselj_zero(n, j) returns 0.0 for some pairs (n, j):

using FunctionZeros

for n in 0:100
    for j in 1:100
        besselj_zero(n, j) == 0.0 && println(n, " ", j)
    end
end

90 6

Other cases can be found by extending the loops to larger values (43 cases for n in 0:1000 and j in 1:1000).

[simonp0420]

This issue is partially fixed in the pending PR #28 , where at least the incorrect (but improved relative to version 1.1) behavior is described in the updated Readme.md . With that PR:

• The first 10 zeros are correctly returned for any choice of order nu.
• For order nu ≤ 150, any choice of zero enumerator j produces a correct result.

Thus, the above code produces no output. When the ranges are extended to 1000 we get:

using Pkg: Pkg
Pkg.activate(temp=true)
Pkg.add(url="https://github.com/simonp0420/FunctionZeros.jl", rev="fix_overflow")
using FunctionZeros: besselj_zero

count::Int = 0
for n in 0:1000
    for j in 1:1000
        if iszero(besselj_zero(n, j))
            println(n, " ", j)
            global count += 1
        end
    end
end
@show count;

Activating new project at `C:\Users\peter\AppData\Local\Temp\jl_ubTerh`
    Updating git-repo `https://github.com/simonp0420/FunctionZeros.jl`
   Resolving package versions...
    Updating `C:\Users\peter\AppData\Local\Temp\jl_ubTerh\Project.toml`
  [b21f74c0] + FunctionZeros v1.1.0 `https://github.com/simonp0420/FunctionZeros.jl#fix_overflow`
    Updating `C:\Users\peter\AppData\Local\Temp\jl_ubTerh\Manifest.toml`
  [7d9f7c33] + Accessors v0.1.43
  [38540f10] + CommonSolve v0.2.6
  [a33af91c] + CompositionsBase v0.1.2
  [187b0558] + ConstructionBase v1.6.0
  [ffbed154] + DocStringExtensions v0.9.5
  [b21f74c0] + FunctionZeros v1.1.0 `https://github.com/simonp0420/FunctionZeros.jl#fix_overflow`
  [3587e190] + InverseFunctions v0.1.17
  [92d709cd] + IrrationalConstants v0.2.6
  [692b3bcd] + JLLWrappers v1.7.1
  [2ab3a3ac] + LogExpFunctions v0.3.29
  [1914dd2f] + MacroTools v0.5.16
  [21216c6a] + Preferences v1.5.1 [loaded: v1.5.0]
  [f2b01f46] + Roots v2.2.10
  [276daf66] + SpecialFunctions v2.6.1
  [efe28fd5] + OpenSpecFun_jll v0.5.6+0
  [56f22d72] + Artifacts v1.11.0
  [ade2ca70] + Dates v1.11.0
  [8f399da3] + Libdl v1.11.0
  [37e2e46d] + LinearAlgebra v1.12.0
  [de0858da] + Printf v1.11.0
  [fa267f1f] + TOML v1.0.3
  [4ec0a83e] + Unicode v1.11.0
  [e66e0078] + CompilerSupportLibraries_jll v1.3.0+1
  [4536629a] + OpenBLAS_jll v0.3.29+0
  [05823500] + OpenLibm_jll v0.8.7+0
  [8e850b90] + libblastrampoline_jll v5.15.0+0
186 13
309 30
347 36
413 47
499 11
502 63
514 20
539 70
557 15
655 93
684 99
703 57
725 60
811 132
844 29
861 16
897 25
912 149
916 35
925 152
938 155
963 95
989 71
997 72
count = 24

[ChristopheBerthod]

Thank you for working on this.

If I may, let me give one reason why these cases are problematic. I suppose that returning 0.0 is a way to indicate a failure of the root-finding algorithm. However, this breaks codes that use the Bessel functions as basis functions, because the normalization uses the zeros. For instance, the set of functions
$$\phi_{ni}(x)=\frac{\sqrt{2}}{RJ_{n+1}(\alpha_{ni})}J_n\left(\alpha_{ni}\frac{x}{R}\right)$$
with $\alpha_{ni}$ the $i$-th zero of $J_n(x)$ form a complete orthonormal basis (relative to index $i$) to represent functions vanishing at the edge of a disk of radius $R$.

Therefore, if the result is anyway wrong in case of failure, why not returning the asymptotic zero rather than 0.0?

[simonp0420]

Thank you for working on this.

You're welcome 😄. My reason for working on this package is because I need the functionality for my own uses: electromagnetic modes of circular metallic waveguides, where the math is very similar to your example. For my application, I only need low-order Bessel functions, where the order nu is always less than 10. The PR I created (#28) ensures correct results for any requested zero whenever nu ≤ 150, which far exceeds my needs.

Therefore, if the result is anyway wrong in case of failure, why not returning the asymptotic zero rather than 0.0?

The value 0.0 that is returned in (only some of!) the incorrect cases is an accidental result of the root-finder routine. Apparently it occurs because the seed value (the asymptotic root) is too far from the true zero sought.

It is difficult to achieve correct results for arbitrary nu and arbitrary enumeration index n because for large nu, the asymptotic formula originally implemented in this package is not accurate enough to use as a starting value until n becomes large enough relative to nu. My PR implements a different asymptotic formula for the complementary (large nu/small n) regime, for which it works well. But to obtain an asymptotic estimate for the n'th Bessel function zero in this regime, it requires as an input the n'th zero of a corresponding Airy function (or its derivative). I've hard-coded the first 10 zeros of the various Airy functions (and derivatives) needed. I did not have the time or energy to create a new package that returns arbitrary zeros of Airy functions and their derivatives. So as it now stands, the original asymptotic formula that assumes large n (relative to nu) must be used as the starting guess for any additional zeros desired beyond the first 10. Here is a scatter plot of the first 15 asymptotic zeros returned by besselj_zero_asymptotic for the J Bessel function of order nu = 725 superimposed over a line plot of the Bessel function:

The code branches so that the first 10 estimates are provided by the new "large nu/small n" asymptotic formula. It's clear that these are very near the actual zero locations. They thus comprise excellent starting points for the root-finder. But estimates starting at the 11th revert to using the original "small nu/large n" formula. You can see that these are very inaccurate, located in some cases almost at the extrema of the function. If we were to return these values to be substituted into your formula for the normalized basis functions, they would result in completely absurd values. So IMO it isn't sufficiently advantageous to return these incorrect estimates.

Using these starting estimates in the root finder, the function besselj_zero produces the following values as the first 15 zeros:

Although all 15 values are located at zeros of the Bessel function, the last five of these occur at the wrong zeros and are therefore incorrect (only 4 incorrect points are plotted because zeros 13 and 14 are both incorrectly found to be 874.7065310820391).

In the future, someone might implement functions for arbitrary Airy function zeros, or use some other method to locate particular zeros of the Bessel functions in the transition region where neither of the implemented asymptotic methods currently provide useful starting points for root finding.

[ChristopheBerthod]

Thank you. I was't aware that the asymptotic expansion is so imprecise at large order. I use these functions for a self-consistent eigenvalue problem (vortices in superconductors), where full convergence requires large order $\sim 500$. A previous code from 2022 gave no errors and I didn't notice a problem, because the final result looked good. But now you are telling me that some of my zeros may have had the wrong index. Good to know :-).

Trying with Mathematica, it gets it right, but its is quite slow.

[simonp0420]

Trying with Mathematica, it gets it right, but its is quite slow.

If you need all of the zeros up to some maximum index, then you could use the following simple method to find a sufficiently accurate seed for a given root from the previous two roots...

using Pkg: Pkg
Pkg.activate(temp=true)
Pkg.add(url="https://github.com/simonp0420/FunctionZeros.jl", rev="fix_overflow")
Pkg.add(["SpecialFunctions", "Roots"])
using FunctionZeros: besselj_zero
using SpecialFunctions: besselj
using Roots: find_zero


"""
    find_all_jzeros!(jzeros::Vector{<:AbstractFloat}, nu::Real) -> jzeros

Fill a preallocated vector `jzeros` with the first `length(jzeros)` nonvanishing zeros of the J Bessel function of order `nu`.

This requires the version of the FunctionZeros package introduced in PR #28, which can be found in the 
"fix_overflow" branch of the repository at https://github.com/simonp0420/FunctionZeros.jl
"""
function find_all_jzeros!(jzeros::Vector{<:AbstractFloat}, nu::Real) 
    nbuiltin = 10 # Number of zeros known to be correct from besselj_zero
    nmax = length(jzeros) # Number of desired zeros

    # Obtain known good zeros:
    for n in 1:min(nmax, nbuiltin)
        jzeros[n] = besselj_zero(nu, n)
    end

    for n in (nbuiltin + 1):nmax
        seed = 2 * jzeros[n-1] - jzeros[n-2]
        jzeros[n] = find_zero(x -> besselj(nu, x), seed)
    end

    return jzeros
end

I tried exercising this for nu = 500, seeking the first 500 zeros:

nu = 500
jzeros = zeros(500)
find_all_jzeros!(jzeros, nu)

It seems to produce the correct list of zeros, and timing it on my machine (Intel i7-14700 running Julia 1.12.2) using BenchmarkTools, it took about 2.5 msec.

If you try this, I'd be interested to see how the results compare to Mathematica for both accuracy and timing.

[ChristopheBerthod]

Thank you for the very good advice. I made some tests and found that your code gives the same zeros as Mathematica to machine accuracy in typically 1000 less time up to order 100, the speedup increasing to about 50'000 at order 750.

When I use your strategy in my self-consistent code, it gets stuck in some cases. I assume that this is an issue with find_zero getting trapped, which I have to investigate.

Nevertheless, I'll close the issue. Thanks again.

[simonp0420]

@ChristopheBerthod I suggest you reopen the issue, since the problem persists with the released version of FunctionZeros. Perhaps having several issues open that are solved by PR #28 will help motivate its prompt review.