Sunday, March 8, 2015

R34. Newton's method in R

Newton's method is based on the approximating a function with only its first Taylor Series term:


y(x+dx) = y(x) + dx y'(x)


We are interested in dx, which will result in y(x + dx) = 0.


0 = y(x) + dx y'(x).


-y(x) = dx y'(x).


dx = -y(x)/y'(x).


A function newton is written implementing this, and used to calculate cube root of 3. Later, we will use it for more complicated cases.

# ex34.R
fit <- function(x) {
  3-x^3
}

newton <- function(f, tol=1E-10, x0=1, N=50, h=0.2) {
  for (i in 1:N) {
    df.dx <- (f(x0+h/2)-f(x0-h/2))/h
    x1 <- x0 - f(x0)/df.dx
    error <- abs(x1-x0)
    cat("i = ",i,", x1 = ",x1,", error = ",error,'\n')
    if (error < tol) break
    x0 <- x1
  }
  x1
}

cube.root.3 <- newton(fit)
cat('cube.root.3 =',cube.root.3,
    '\nfit(cube.root.3) = ',fit(cube.root.3))
# i =  1 , x1 =  1.664452 , error =  0.6644518 
# i =  2 , x1 =  1.470826 , error =  0.1936256 
# i =  3 , x1 =  1.442844 , error =  0.02798194 
# i =  4 , x1 =  1.442251 , error =  0.0005935402 
# i =  5 , x1 =  1.44225 , error =  1.193565e-06 
# i =  6 , x1 =  1.44225 , error =  1.910613e-09 
# i =  7 , x1 =  1.44225 , error =  3.056666e-12 
# cube.root.3 = 1.44225 
# fit(cube.root.3) =  -3.153033e-14

No comments:

Post a Comment