Increase robustness of zero-finding and document the valid range of usage by simonp0420 · Pull Request #28 · JuliaMath/FunctionZeros.jl

simonp0420 · 2025-10-17T22:07:37Z

Add methods to bessel_zero_asymptotic and bessel_deriv_zero_asymptotic that promote nu from Integer to float. Fixes Issue #27

simonp0420 · 2025-10-17T22:30:31Z

I submitted this PR a little too early. All tests ran successfully on my machine where Int === Int64, but I see some have failed in CI on x86 machines, where Int === Int32, and the @test_broken tests are actually no longer broken. I also now see from the comments preceding the "high nu" test set that Issue #27 overlaps with Issue #10. I'll modify this PR to replace the @test_broken tests to @test after which I think everything should work ok.

simonp0420 · 2025-10-17T23:18:10Z

Edit: Changing this to a WIP (Work In Progress) PR since the problem still occurs for even higher nu values. Please don't merge yet.

simonp0420 · 2025-10-19T23:59:46Z

Ok, I've implemented an asymptotic formula from DLMF that is appropriate for large nu and small n. Now the first ten zeros of any of the four treated functions are returned correctly for any positive nu value. All zeros are returned correctly for nu less than or equal to at least 150. This corrects the problem mentioned in this comment to Issue #10.

Here is an example of finding the correct first 25 zeros for the derivative of a Bessel Y function of order 150:

julia> using FunctionZeros, SpecialFunctions, Plots

julia> t = [bessely_deriv_zero(150, n) for n in 1:25]
25-element Vector{Float64}:
 159.85254354288477
 167.76036070522397
 174.32071688833557
 180.22537114717298
 185.71690058288073
   ⋮
 255.3186728810697
 259.1858488958748
 263.0241588781523
 266.8356390609487
 270.6221130091581

julia> bessely_deriv(nu, x) = 0.5 * (bessely(nu-1, x) - bessely(nu+1, x))
bessely_deriv (generic function with 1 method)

julia> maximum(x -> abs(bessely_deriv(150, x)), t)
8.847089727481716e-16

julia> plot(;legend=nothing); plot!(range(140, 272, 600), x -> bessely_deriv(150, x)); scatter!(t, zeros(length(t)))

And here is an example showing that the first 10 zeros are correctly found for a ridiculously large value of nu = 10_000:

julia> t = [bessely_deriv_zero(10_000, n) for n in 1:11]
11-element Vector{Float64}:
 10039.277992159914
 10069.79405626667
 10094.526103710305
 10116.360982858023
 10136.3273778352
 10154.945576301201
 10172.525286577882
 10189.26979494539
 10205.321627357598
 10220.78559381645
 11243.12462090691

julia> maximum(x -> abs(bessely_deriv(10_000, x)), t[1:10])
3.3608255067818504e-12

julia> plot(;legend=nothing); plot!(range(10038, 11245, 600), x -> bessely_deriv(10_000, x)); scatter!(t, zeros(length(t)))

Removed WIP from title of PR.
Ready for review

simonp0420 · 2025-12-30T16:55:28Z

Friendly bump (maybe reviewing this would make a good new year's resolution? 😉).

simonp0420 · 2026-01-31T22:18:39Z

@jlapeyre Could you please review? Or is there anyone else who could be asked to review? TIA

simonp0420 · 2026-04-25T13:04:14Z

After multiple friendly bumps over several months, I've gotten zero responses. @DilumAluthge Is there anyone else in JuliaMath who could possibly review this PR?

DilumAluthge · 2026-04-25T15:49:24Z

Hi @simonp0420, I'd be happy to help. Let me see if I can recruit someone to review.

DilumAluthge · 2026-04-25T15:49:47Z

In the meantime, could you update the PR title to be a bit more descriptive?

simonp0420 · 2026-04-25T16:55:11Z

To whomever reviews this PR: Please also take a look at Issue #30 which is addressed by this PR, and which contains a lot of useful background info in the comments. Note that Issue #29 is also fixed by this PR.

jlapeyre · 2026-04-25T17:10:57Z

Thanks, @simonp0420, for all the work. And thanks @DilumAluthge for the ping.

oscardssmith · 2026-04-25T17:14:37Z

 function bessel_deriv_zero_asymptotic(nu_in::Real, n::Integer, kind=1)
-    # Reference: https://dlmf.nist.gov/10.21.E20
    nu = abs(nu_in)
+    if (nu ≥ 33 && n ≤ 10) || (nu ≥ 30 && n ≤ 9) || (nu ≥ 26 && n ≤ 8) || (nu ≥ 25 && n ≤ 7)


could this be if nu>3n? or something similar?

Yes, thanks. Will replace with (nu ≥ 25) && (nu ≥ 3.25 * n)

oscardssmith · 2026-04-25T17:18:16Z

+   package. Each field in the named tuple consists of a tuple of `n` increasingly negative values. Here `n`
+   is a small integer, currently 10. 
+"""
+@inline function airy_zeros()


where do these numbers come from?

I'm embarassed to admit that were obtained by first plotting the airy functions, then using Roots.find_zero with a seed value obtained for each root by visually inspecting the plot. Here is some validation:

julia> using FunctionZeros: FunctionZeros as fz [ Info: Precompiling FunctionZeros [b21f74c0-b399-568f-9643-d20f4fa2c814](cache misses: incompatible header (9)) Precompiling FunctionZeros finished. 7 dependencies successfully precompiled in 8 seconds. 19 already precompiled. julia> (; ai, bi, aiprime, biprime) = fz.airy_zeros(); julia> using SpecialFunctions julia> maximum(x -> abs(airyai(x)), ai) 2.0278943050169084e-15 julia> maximum(x -> abs(airybi(x)), bi) 1.2350493860255089e-15 julia> maximum(x -> abs(airyaiprime(x)), aiprime) 2.905936767899529e-15 julia> maximum(x -> abs(airybiprime(x)), biprime) 4.585753687468294e-15

We should definitely do the calculation of the roots in BigFloat precision (and then round that down to Float64). (possibly using mathematica if Special Functions doesn't support it). These functions aren't that useful if they aren't hitting the true zeros.

I'll push back gently on this...

I had considered obtaining these to higher precision when I wrote the function, but since this function is not exported nor externally documented, and its used only for obtaining approximate starting values for finding Bessel function roots, I didn't think that it was necessary. I figured that eventually, someone will write a proper package to compute any number of airy function zeros to arbitrary precision. Perhaps it's sufficient to state in the docstring that the zeros are only supplied to Float64 accuracy.

If you still think it's important to supply BigFloat precision for the airy function zeros here then I'll implement it.

in that case, this is probably fine. I thought that these values would be directly affecting the output.

I would have trouble implementing it anyway, since I don't have Mathematica and I just found out that SpecialFunctions.airybi is erroring for negative BigFloat arguments.

oscardssmith · 2026-04-25T17:19:11Z

+               -7.940178689168584, -9.019583358794248, -10.037696334908555, -11.00646266771229,
+               -11.934261645014844, -12.82725830917722)
+    #return (; ai, bi, aiprime, biprime) # Not compatible with Julia 1.0
+    return (ai=ai, bi=bi, aiprime=aiprime, biprime=biprime) # Compatible with Julia 1.0


given that users of this function only use 1 set of these, the interface seems slightly awkward.

To obtain the zeros of airyai by writing airy_zeros().ai doesn't seem too onerous to me, and since it's a tuple there aren't any allocations involved in the call. That said, I'm open to suggestions for a more elegant interface.

jlapeyre · 2026-04-25T17:20:14Z

There is some trailing whitespace in src/FunctionZeros.jl that should be removed.

oscardssmith · 2026-04-25T17:21:48Z

+    for alphak in (alpha1, alpha2, alpha3, alpha4, alpha5)
+        xpk *= x
+        nextterm = alphak * xpk
+        abs(nextterm) > abs(lastterm) && break # Asymptotic series starting to diverge


is this early break actually faster? this loop has a maximum of 5 iters

The early break isn't intended to reduce execution time but rather the error in the partial sum. Since it's an asymptotic series, the "truncation error" can start to grow with additional terms, unlike a convergent series.

* Remove trailing white spaces

simonp0420 · 2026-04-25T21:23:51Z

@jlapeyre Sorry about the white space--don't know why it keeps crawling in there, but I think it's now gone and I've configured VS Code to eliminate it from now on.
@oscardssmith I simplified the test for which asymptotic series to use, iaw your suggestion.

simonp0420 · 2026-05-02T16:40:58Z

Any other changes you'd like to see?

oscardssmith · 2026-05-02T21:00:44Z

lgtm

simonp0420 · 2026-05-10T23:14:29Z

@jlapeyre If there are no other changes you desire, would you consider merging this PR? TIA!

DilumAluthge · 2026-05-22T19:21:06Z

@oscardssmith Is this good to merge from your point of view?

DilumAluthge · 2026-05-22T19:22:46Z

If so, I can merge this and make a new release.

oscardssmith · 2026-05-22T19:33:19Z

sgtm

simonp0420 · 2026-05-22T19:43:49Z

Thanks, @DilumAluthge and @oscardssmith!

Fix Issue JuliaMath#27

494816f

Convert @test_broken to @test

4410a5e

simonp0420 changed the title ~~Fix Issue #27~~ WIP: Fix Issue #27 Oct 17, 2025

simonp0420 added 4 commits October 19, 2025 15:47

Add Limitations section to README.md

fa52c85

Add additional tests for large nu values

890a138

Implement asymptotic formulas valid for large nu and small n

2c26ca7

Fix named tuple constructor for Julia 1.0 compatibility

2ddbd98

simonp0420 changed the title ~~WIP: Fix Issue #27~~ Fix Issue #27 Oct 20, 2025

simonp0420 added 2 commits October 20, 2025 10:58

Increase size of table for airy function zeros to 10

c0d7ba3

Adjust test for which asymptotic formula to use.

d585e2d

simonp0420 mentioned this pull request Dec 25, 2025

When nu is large enough, the root finding fails. #29

Open

simonp0420 mentioned this pull request Jan 7, 2026

besselj_zero failing to return some zeros (v1.1.0) #30

Open

simonp0420 changed the title ~~Fix Issue #27~~ Increase robustness of zero-finding and document the valid range of usage Apr 25, 2026

DilumAluthge requested a review from oscardssmith April 25, 2026 17:06

oscardssmith reviewed Apr 25, 2026

View reviewed changes

* Simplify test for which asymptotic series to use

2bf03f0

* Remove trailing white spaces

Remove trailing spaces

b8f69f3

DilumAluthge merged commit 6086d6d into JuliaMath:master May 22, 2026
15 checks passed

simonp0420 deleted the fix_overflow branch May 22, 2026 19:42

Conversation

simonp0420 commented Oct 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simonp0420 commented Oct 17, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simonp0420 commented Oct 17, 2025

Uh oh!

simonp0420 commented Oct 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

simonp0420 commented Dec 30, 2025

Uh oh!

simonp0420 commented Jan 31, 2026

Uh oh!

simonp0420 commented Apr 25, 2026

Uh oh!

DilumAluthge commented Apr 25, 2026

Uh oh!

DilumAluthge commented Apr 25, 2026

Uh oh!

simonp0420 commented Apr 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jlapeyre commented Apr 25, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

simonp0420 Apr 25, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jlapeyre commented Apr 25, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

simonp0420 commented Apr 25, 2026

Uh oh!

simonp0420 commented May 2, 2026

Uh oh!

oscardssmith commented May 2, 2026

Uh oh!

simonp0420 commented May 10, 2026

Uh oh!

DilumAluthge commented May 22, 2026

Uh oh!

DilumAluthge commented May 22, 2026

Uh oh!

oscardssmith commented May 22, 2026

Uh oh!

Uh oh!

simonp0420 commented May 22, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

simonp0420 commented Oct 17, 2025 •

edited

Loading

simonp0420 commented Oct 17, 2025 •

edited

Loading

simonp0420 commented Oct 19, 2025 •

edited

Loading

simonp0420 commented Apr 25, 2026 •

edited

Loading

simonp0420 Apr 25, 2026 •

edited

Loading