如何避免停止循环的 uniroot 错误？

Question

我在循环中运行 uniroot 函数，但遇到错误，代码停止。代码如下;

func <-function(f) -a*b/c*0.5*d*e^2 + (d/f-1)*g*sin(h*(pi/180))-i
dat <- data.frame(a = c(0.99,0.99,0.99),
                  b = c(0.1986572,0.1986572,0.1986572),
                  c = c(237.5,237.5,237.5),
                  d = c(1028.372, 1028.711, 1028.372),
                  e = c(2.46261, 2.986461, 2.46261),
                  f = c(-1,-1,-1),
                  g = c(9.8,9.8,9.8),
                  h = c(-54.97964, -51.65978, -54.97964),
                  i = c(0.03699588, -0.0375189, 0.03699588))

for(j in 1:length(dat$a)){
   a <- dat$a[j]
   b <- dat$b[j]
   c <- dat$c[j]
   d <- dat$d[j]
   e <- dat$e[j]
   #f: this should be solved by uniroot
   g <- dat$g[j]
   h <- dat$h[j]
   i <- dat$i[j]
   sol <- uniroot(func,c(0, 2000),extendInt = "yes") 
   dat$f[j] <- sol$root
   print(j)
}

运行以上代码，出现以下错误：

[1] 1
Error in uniroot(func, c(0, 2000), extendInt = "yes") : 
      no sign change found in 1000 iterations

代码停在了j=1，并没有走到j=2 & 3。因此，dat$f 显示

> dat$f
[1] 1526.566   -1.000   -1.000

我的目标是当 uniroot 在给定的 j 中遇到错误时，将 NA 放入 dat$f[j]，并在结束时继续循环。

如果这有效，dat$f[1] 和 dat$f[3] 应该使用上述数据帧具有相同的值 (=1526.566)。

请告诉我如何处理 uniroot 错误。

Answer 1

如果下限设置为 1 而不是 0，则问题中的代码将起作用。问题是如果 f 为 0，则 func 未定义，因为 f 在分母中。

虽然这已经足够了，但建议您进行以下更改：

使用末尾注释中显示的数据，确保它不包括 f。
如图所示更紧凑地定义 func
使用try允许uniroot继续运行即使有错误（虽然在他的例子中有none）
使用下限 1，因为 func 在 0 处未定义
使用参数 uniroot
将结果放在 res 中（如果迭代失败则放在 NA 中）这样就不会混淆什么是输入什么是输出。

使用所示的更紧凑的 func 形式并将结果放入 res。

func <- function(f, data) with(data, -a*b/c*0.5*d*e^2 + (d/f-1)*g*sin(h*(pi/180))-i)
nr <- nrow(dat)
res <- numeric(nr)
for(j in 1:nr){
   sol <- try(uniroot(func, c(1, 2000), data = dat[j, ], extendInt = "yes") )
   res[j] <- if (inherits(sol, "try-error")) NA else sol$root
   print(j)
}
## [1] 1
## [1] 2
## [1] 3
res
## [1] 1526.566 2014.476 1526.566

备注

dat <- data.frame(a = c(0.99,0.99,0.99),
                  b = c(0.1986572,0.1986572,0.1986572),
                  c = c(237.5,237.5,237.5),
                  d = c(1028.372, 1028.711, 1028.372),
                  e = c(2.46261, 2.986461, 2.46261),
                  g = c(9.8,9.8,9.8),
                  h = c(-54.97964, -51.65978, -54.97964),
                  i = c(0.03699588, -0.0375189, 0.03699588))

如何避免停止循环的 uniroot 错误？

How can I avoid uniroot error that stops the loop?

r

try-catch

uniroot

备注