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

Leave a Reply

Your email address will not be published. Required fields are marked *