Solve BesselJ (J0) function in R

A very brute-force (yet functional) approach would be to obtain the roots by minimization:

# tightly spaced grid
init_g <- seq(0,100,by = 1)

# objective function: sum of squares of the values of the Bessel function
J0 <- function(x) sapply(x, function(y) besselJ(nu = 0, y))
objective <- function(x) {sum(J0(x)^2)}

# minimize the objective functional with respect to the initial guess
# round the values to eliminate numerical uncertainty and output unique values
bessel_roots <- unique(round(nlm(objective, p=init_g, ndigit=4)$estimate, 4))

Output

> bessel_roots
 [1]  2.4048  5.5201  8.6537 11.7915 14.9309 18.0710 18.0711 21.2116 24.3525
[10] 27.4935 30.6346 33.7758 36.9171 40.0584 43.1997 43.1998 46.3412 49.4826
[19] 52.6240 55.7655 58.9070 62.0484 65.1899 68.3314 71.4730 71.4729 74.6145
[28] 77.7560 80.8975 84.0391 84.0390 87.1806 90.3221 93.4637 96.6052 99.7468

I'm a visual thinker, so here's how I would approach this. First of all, set up a sequence of x values between, say, 1 and 100. Then get the value of BesselJ with n = 0 at each point:

x <- seq(1, 100, 0.1)
y <- besselJ(x, 0)

Now put them in a data frame and plot them:

df <- data.frame(x, y)

library(ggplot2)

p <- ggplot(df, aes(x, y)) + geom_line()

p

Now we can see that wherever the roots are, there will be a change in sign of the y value, so let's identify the x values where the y values change sign. We can use these as points around which to look. We know, since our sequence is 0.1 apart, each root must fall between -0.1 and +0.1 of these points:

lower <- x[which(diff(sign(y)) != 0)] - 0.1

So we can use these ranges as inputs to uniroot, provided we use sapply to feed the ranges into uniroot one-by-one.

bes <- function(x) besselJ(x, 0)
roots <- sapply(lower, function(z) uniroot(bes, interval = c(0, 0.2) + z)$root)

We now have our roots between 1 and 100:

roots
#>  [1]  2.404822  5.520077  8.653728 11.791554 14.930917 18.071084 21.211636 24.352494
#>  [9] 27.493463 30.634588 33.775837 36.917098 40.058441 43.199791 46.341164 49.482615
#> [17] 52.624066 55.765511 58.906984 62.048469 65.189983 68.331441 71.472982 74.614519
#> [25] 77.756026 80.897543 84.039092 87.180630 90.322155 93.463719 96.605241 99.746820

And we can show these are correct by adding them to our plot:

p + geom_point(data = data.frame(x = roots, y = besselJ(roots, 0)),
               colour = "red")

In terms of speed, this will get you the first 1000 roots in under a 10th of a second:

bes1000 <- function()
{
  x <- seq(1, 3143, 0.1)
  y <- besselJ(x, 0)
  lower <- x[which(diff(sign(y)) != 0)] - 0.1
  bes <- function(x) besselJ(x, 0)
  roots <- sapply(lower, function(z) uniroot(bes, interval = c(0, 0.2) + z)$root)
}

microbenchmark::microbenchmark(bes1000())
#> Unit: milliseconds
#>       expr     min       lq     mean   median       uq      max neval
#>  bes1000() 78.1617 81.22905 88.62894 83.09845 90.67025 396.5271   100

^{Created on 2020-07-15 by the reprex package (v0.3.0)}