Newton's Method in R

Question

I have an issue when trying to implement the code for Newton's Method for finding the value of the square root (using iterations). I'm trying to get the function to stop printing the values once a certain accuracy is reached, but I can't seem to get this working. Below is my code.

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

Why won't the if-statement work?

Wendiwendie

Change it to abs(GUESS^2-x) – Nannette

@Wendiwendie You should really separate the task of getting user input (readline) from the task of doing the calculation. Functions should do one task, or your code gets confusing. – Strafe

Why are you converting GUESS to int? It just goes back to a double first time thru the loop anyway. – Swashbuckler

Answer 1

UPDATE:

Please see @RichieCotton's comment to @agstudy's answer. I agree with Richie, and in fact it makes more sense to use @agstudy's approach.

---

Original answer:

Your function is fine, your math is off.
GUESS and x should not (necessarilly) be close, but GUESS * GUESS and x should be.

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
}

Ralston

"it makes more sense to use @agstudy's approach". And for non-homework problems, it makes more sense to use a robust, well-tested function built into R like nlm or optimize or optim. – Strafe

Answer 2

This should work:

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

Earnest

It took a bit of puzzling to work out the difference between this answer and Ricardo's. This checks convergence by comparing the current iteration to the last one (abs(GUESS-myvector[i]) < eps), whereas Ricardo's answer compares the current answer squared to the input (abs((GUESS*GUESS)-x) < eps). Both techniques work; I think this is slightly more standard. – Strafe

@RichieCotton well, there's always a risk in comparing successive elements in a sequence. Some functions converge very slowly indeed, and you may well get false positives this way. Comparing w/ x^2 will always get you a correct evaluation of the convergence. – Swashbuckler

@CarlWitthoft For sqrt it doesn't matter. But breaking on "converging to something" may be better than "converging on the right answer" for cases where, e.g., you get stuck at a stationary point. Better to stop calculating, rather than repeating the same incorrect number over and over until you reach itmax. – Strafe

@RichieCotton I see your point. I guess I'd prefer to use a "sensible" value of itmax and examine the results if convergence failure is reported. – Swashbuckler

@RichieCotton, you make a great point wrt '.. But breaking on "converging to something" may be better than "converging on the right answer"'. In that context, this approach makes more sense. – Ralston