6

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?

SuperPrograman
  • 1,716
  • 15
  • 23
user2884679
  • 129
  • 3
  • 9

2 Answers2

3

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
}
Ricardo Saporta
  • 51,025
  • 13
  • 129
  • 166
  • 1
    "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`. – Richie Cotton Oct 16 '13 at 14:08
3

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
agstudy
  • 113,354
  • 16
  • 180
  • 244
  • 1
    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. – Richie Cotton Oct 16 '13 at 09:59
  • 1
    @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. – Carl Witthoft Oct 16 '13 at 12:05
  • 1
    @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`. – Richie Cotton Oct 16 '13 at 13:10
  • @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. – Carl Witthoft Oct 16 '13 at 13:27
  • @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. – Ricardo Saporta Oct 16 '13 at 13:37