解决R中的BesselJ（J0）函数

debugcn 发表于 Dev

奥罗努斯一世

我正在寻找一种解决R中BesselJ函数的方法。我知道有无限的可能性，但比方说，我想知道100甚至1000个第一个值……可能吗？

艾伦·卡梅伦

我是一个视觉思想家，所以这是我的处理方法。首先，设置一个介于1到100之间的x值序列。然后在每个点上获取n = 0的BesselJ值：

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

现在将它们放在数据框中并绘制它们：

df <- data.frame(x, y)

library(ggplot2)

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

p

现在我们可以看到，无论根在哪里，y值的符号都会发生变化，因此让我们确定y值更改符号的x值。我们可以将它们用作观察点。我们知道，因为我们的序列相距0.1，所以每个根必须落在这些点的-0.1和+0.1之间：

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

因此，我们可以使用这些范围作为输入uniroot，只要我们使用sapply到的范围，送入uniroot一个接一个。

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

现在我们的根在1到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

我们可以通过将它们添加到绘图中来证明它们是正确的：

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

在速度方面，这将使您在不到十分之一秒的时间内获得前1000个根：

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

^{由reprex软件包（v0.3.0）创建于2020-07-15}

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。