Repository navigation
besselj_zero failing to return some zeros (v1.1.0) #30
Description
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
end90 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).
Activity
simonp0420 commentedon Jan 7, 2026
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:
Thus, the above code produces no output. When the ranges are extended to 1000 we get:
ChristopheBerthod commentedon Jan 8, 2026
Thank you for working on this.
If I may, let me give one reason why these cases are problematic. I suppose that returning
$$\phi_{ni}(x)=\frac{\sqrt{2}}{RJ_{n+1}(\alpha_{ni})}J_n\left(\alpha_{ni}\frac{x}{R}\right)$$ $\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$ .
0.0is 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 functionswith
Therefore, if the result is anyway wrong in case of failure, why not returning the asymptotic zero rather than
0.0?simonp0420 commentedon Jan 9, 2026
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.
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_asymptoticfor 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_zeroproduces 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 commentedon Jan 9, 2026
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 commentedon Jan 9, 2026
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...
I tried exercising this for nu = 500, seeking the first 500 zeros:
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 commentedon Jan 19, 2026
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_zerogetting trapped, which I have to investigate.Nevertheless, I'll close the issue. Thanks again.
simonp0420 commentedon Jan 19, 2026
@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.