R中的牛顿法

Posted 2023-03-12

【中文标题】R中的牛顿法

【英文标题】：Newton's Method in R

【发布时间】：2013-10-24 00:21:59

【问题描述】：

我在尝试实现用于查找平方根值的牛顿法的代码（使用迭代）时遇到问题。一旦达到一定的准确性，我试图让函数停止打印值，但我似乎无法让它工作。下面是我的代码。

MySqrt <- function (x, eps = 1e-6, itmax = 100, verbose = TRUE)
  i <- 1
  myvector <- integer(0)
  GUESS <- readline(prompt="Enter your guess: ")
  GUESS <- as.integer(GUESS)
  while(i <= itmax)
      GUESS <- (GUESS + (x/GUESS)) * 0.5
      myvector <- c(myvector, GUESS)
      if (abs(GUESS-x) < eps) break
      i <- i + 1
  

  myvector

为什么 if 语句不起作用？

【问题讨论】：

改成abs(GUESS^2-x)

@user2884679 您应该真正将获取用户输入 (readline) 的任务与进行计算的任务分开。函数应该完成一项任务，否则您的代码会变得混乱。

为什么要将GUESS 转换为int？无论如何，它只是通过循环回到第一次的两倍。

【参考方案1】：

这应该可行：

MySqrt <- function (x, eps = 1e-6, itmax = 100, verbose = TRUE)
  i <- 1
  myvector <- vector(mode='numeric',itmax)  ## better to allocate memory
  GUESS <- readline(prompt="Enter your guess: ")
  GUESS <- as.numeric(GUESS)
  myvector[i] <- GUESS
  while(i <= itmax)
    GUESS <- (GUESS + (x/GUESS)) * 0.5
    if (abs(GUESS-myvector[i]) < eps) break
    i <- i + 1
    myvector[i] <-  GUESS
  
  myvector[seq(i)]


MySqrt(2)
Enter your guess: 1.4
[1] 1.400000 1.414286 1.414214

【讨论】：

解决这个答案和里卡多的答案之间的区别有点令人费解。这通过将当前迭代与最后一个迭代 (abs(GUESS-myvector[i]) < eps) 进行比较来检查收敛性，而 Ricardo 的答案将当前答案与输入 (abs((GUESS*GUESS)-x) < eps) 的平方进行比较。两种技术都有效；我认为这更标准一些。

@RichieCotton 好吧，比较序列中的连续元素总是有风险的。有些函数确实收敛得很慢，这样你很可能会得到误报。与x^2 进行比较总能让您正确评估收敛性。

@CarlWitthoft 对于 sqrt 没关系。但是，在例如你卡在一个固定点的情况下，打破“收敛到某事”可能比“收敛到正确的答案”更好。最好停止计算，而不是一遍又一遍地重复相同的错误数字，直到达到itmax。

@RichieCotton 我明白你的意思。我想我更喜欢使用 itmax 的“合理”值，并在报告收敛失败时检查结果。

@RichieCotton，你在'.. But breaking on "converging to something" may be better than "converging on the right answer"' 中提出了很好的观点。在这种情况下，这种方法更有意义。

【参考方案2】：

更新：

请参阅@RichieCotton 对@agstudy 回答的评论。我同意 Richie 的观点，事实上使用@agstudy 的方法更有意义。

---

原答案：

你的函数很好，你的数学是关闭的。GUESS 和 x 不应该（必然）接近，但 GUESS * GUESS 和 x 应该接近。

MySqrt <- function (x, eps = 1e-6, itmax = 100, verbose = TRUE)
  i <- 1
  myvector <- integer(0)
  GUESS <- readline(prompt="Enter your guess: ")
  GUESS <- as.integer(GUESS)
  while(i <= itmax)
      GUESS <- (GUESS + (x/GUESS)) * 0.5
      myvector <- c(myvector, GUESS)
      browser(expr=i == 10 || abs(GUESS-x) < eps)
      if (abs((GUESS*GUESS)-x) < eps) break    ###  <~~~~  SEE HERE
      i <- i + 1
  

  myvector

【讨论】：

“使用@agstudy 的方法更有意义”。对于非家庭作业问题，使用 R 内置的强大、经过充分测试的函数更有意义，例如 nlm 或 optimize 或 optim。