在R中计算when循环时需要TRUE/FALSE时缺少的值

Question

所以我试图实现一种算法，用牛顿的方法来求某个给定函数的根，我认为，除了在一个很小但非常令人沮丧的点外，我认为是相当成功的。

当我为x_0的小初始值(即小于或接近1，-1.1)计算函数时，我的函数似乎工作得很好，但是当我试图将它用于更大的值(甚至更高的值)时，我会得到以下错误消息：

Error in while (abs(x_1 - x_0) > 10^(-1 * k) && iter < max_iter) { : 
  missing value where TRUE/FALSE needed

下面是我针对这个问题的代码：

g_x <- function(x){
  return(x/(1+x^2)^(1/2))
}

dg_x <- function(x){
  # derivative of the given function of g
  return(1/(x^2+1)^(3/2))
}

newton <- function(x_0,k){
  # x_0 is initial guess from which you want to being zeroing in on the root
  # k is the parameter of the tolerance = 10^(-k)
  if(between(g_x(x_0),-10^(-1*k)/2,10^(-1*k)/2)){
    return(c(x_0,0))
  }
  
  iter <- 1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  
  while(abs(x_1-x_0) > 10^(-1*k)){
    x_0 <- x_1
    x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
    iter <- iter + 1
  }
  return(c((x_1+x_0)/2, iter))
}

我不知道该怎么做才能解决这个问题。我一直在控制台中逐行查看此过程，计算机能够将abs(x_1-x_0) > 10^(-k)表达式计算为对于第一对迭代是正确的。我不知道为什么这会成为一个问题。

由于该问题在算法的第一次迭代中就出现了，因此它不像计算机将x_1的某些超小/大的值概括为零或无穷大。

在这个问题中，我们给出了: x_0 =2和g(x) = x/sqrt(x^2+1)所以dg(x) = 1/(x^2+1)^3/2

因此，下一个要计算的点是x_1 = -8

所以计算机正在计算abs(-8-2) > 10^-k.为什么这个破了

Answer 1

在进行非线性优化问题之前，我通常会看一下我想用图求解的函数。因为您的函数只依赖于一个变量，所以这样做很简单。

xs = seq(-10,10,length=10000)
plot(xs,g_x(xs))

收益率

​

​

因此，‘根’在0(函数的值为0)，但导数“非常尖锐化”，这可能意味着牛顿方法会反复重复地射出解。

我修改了您的函数，使其为x设置了一个上界和下界，而且只要您开始相当接近解决方案，它就能工作：

newton <- function(x_0,k,lower=-10,upper=10){
  # x_0 is initial guess from which you want to being zeroing in on the root
  # k is the parameter of the tolerance = 10^(-k)
  if(between(g_x(x_0),-10^(-1*k)/2,10^(-1*k)/2)){
    return(c(x_0,0))
  }
  
  iter <- 1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  
  while( (abs(x_1-x_0) > 10^(-1*k)) & between(x_0,lower,upper)){
    x_0 <- x_1
    x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
    iter <- iter + 1
  }
  return(c((x_1+x_0)/2, iter))
}

#x_0=10
#k=10
newton(.3,10)

newton(5,10)

从0.3个开始工作(解决-3.812802e-15 )，从5开始不工作(出口在976500)。使用粗网格搜索为您的函数选择一个好的起始值是一种合理的方法。

这显示了如果你从2开始的话会发生什么：

x_0=2  
iter <- 1
x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))

for (inter in 1:10){
  x_0 <- x_1
  x_1 <- x_0 - (g_x(x_0)/dg_x(x_0))
  iter <- iter + 1
  print(x_1,iter)
}

返回：

[1] 512
[1] -1.34e+08
[1] 2.418e+24
[1] -1.4135e+73
[1] 2.82401e+219
[1] NaN
[1] NaN
[1] NaN
[1] NaN
[1] NaN