Approximation of a root using Halley's method in R

Question

Alright, so I'm working on a small R program in order to do approximation using Halley's method. Basically I need to be able to approximate to 9 decimal places the value of 59^(1/7) using Halley's method.

What I have for the first order recurrence relation of Halley's method is this:

Xn+1 = Xn - [ f(Xn) / ( f'(Xn) - (f(Xn)f''(Xn)/2f'(Xn)) ) ]

So far this is the code I have.

halleysMethodApprox = function(ftnH, x0, m0, k0, tol = 1e-8, max.iter=2) {
    x <- x0
    m <- m0
    k <- k0
    fx <- ftnH(x, m, k)
    iter <- 0
    b <- fx[1]/(fx[2] - (fx[1]*fx[3])/(2*fx[2]) )
    while( (abs(fx[1] - x) > tol) && (iter < max.iter) ) {
            # calculate X(n+1)
            b <- ( fx[1]/(fx[2] - ( (fx[1]*fx[3])/(2*fx[2]) ) ))
            x <- x - b
            fx <- ftnH(x, m-1, 0)
            iter <- iter + 1
            cat("At iteration", iter, "value of x is: ", x, "\n")
    }
    if( abs(x) > tol ) {
            cat("Algorithm failed to converge\n")
            return(NULL)
    } else {
            cat("Algorithm converged\n")
            return(x)
    }
}

and this function for generating a vector containing the function of x, and its derivatives.

ftnH = function (x, m, k) {
        fx <- x^m - k
        dfx <- (m*x^(m-1))
        ddfx <- ((m-1)*m)*x^(m-2)
        return (c(fx, dfx, ddfx))
}


print(halleysMethodApprox(ftnH, 59, (1/7), 0))

I'm not quite sure how I am supposed to numerically approximate 59^(1/7) using the above definition for Halley's method.

It's very unclear to me what your question is. What programming challenge are you facing? — MrFlick, Aug 07 '14 at 21:48
The programming challenge I am facing is rewriting that first order recurrence relation in a way that allows me to approximate the value of a function f(x) x^m - k where x = 59 and m = 1/7 — EthanLWillis, Aug 07 '14 at 21:51

score 0 · Answer 1 · answered Aug 12 '14 at 09:17

If you want to calculate x=59^(1/7), your function for using Halley's method is f(x) = x^7-59 = 0. And there were a few typos and minor errors in your code. Here's a corrected version:

halleysMethodApprox = function(ftnH, x0, m0, k0, tol = 1e-8, max.iter=100) {
  # x0: starting estimate, eqn: x = k0^(1/m0) 
  x <- x0
  m <- m0
  k <- k0
  fx <- ftnH(x0, m0, k0)
  iter <- 0
  b <- fx[1]/(fx[2] - (fx[1]*fx[3])/(2*fx[2]) )
  while( (abs(fx[1]) > tol) && (iter < max.iter) ) {
    # calculate X(n+1)
    b <- ( fx[1]/(fx[2] - ( (fx[1]*fx[3])/(2*fx[2]) ) ))
    x <- x - b
    fx <- ftnH(x, m0, k0)
    iter <- iter + 1
    cat("At iteration", iter, "value of x is: ", x, "\n")
  }
  if(abs(fx[1]) > tol) {
    cat("Algorithm failed to converge\n")
    return(NULL)
  } else {
    cat("Algorithm converged\n")
    return(x)
  }
}

ftnH = function (x, m, k) {
  fx <- x^m - k
  dfx <- (m*x^(m-1))
  ddfx <- ((m-1)*m)*x^(m-2)
  return (c(fx, dfx, ddfx))
}

halleysMethodApprox(ftnH, 1, 7, 59) 
 # At iteration 1 value of x is:  1.320442 
 # At iteration 2 value of x is:  1.655396 
 # At iteration 3 value of x is:  1.78716 
 # At iteration 4 value of x is:  1.790519 
 # At iteration 5 value of x is:  1.790519 
 # Algorithm converged
 # [1] 1.790519

Approximation of a root using Halley's method in R

1 Answers1