Here is a toy example of implementing Newton’s method in R. I found some old code that I had written a few years ago when illustrating the difference between convergence properties of various root-finding algorithms, and this example shows a couple of nice features of R.

Newtons method is an iterative root-finding algorithm where we start with an initial guess \(x_0\) and each successive guess is refined using the iterative process:

$$ x_{n+1} = x_n – \frac{f(x_n)}{f’(x_n)} $$

Here is a simple implementation of this in R:

newton <- function(f, tol=1E-12,x0=1,N=20) { h <- 0.001 i <- 1; x1 <- x0 p <- numeric(N) while (i<=N) { df.dx <- (f(x0+h)-f(x0))/h x1 <- (x0 - (f(x0)/df.dx)) p[i] <- x1 i <- i + 1 if (abs(x1-x0) < tol) break x0 <- x1 } return(p[1:(i-1)]) }

Note a couple of things:

* The step-size for numerical differentiation is hardcoded to be 0.001. This is arbitrary and should probably be a parameter.

* The algorithm will run until either the number of steps N has been reached, or the error tolerance \(\left|x_{n+1}-x_n\right|< \epsilon\), where \(\epsilon\) is defined as the tolerance parameter tol.

* The function returns a vector of the iterated x positions, which will be <= N.

Here is a simple example: for instance, to find the zeros of the function

$$ f(x) = x^3 + 4x^2 - 10 $$

which in R we can represent as:

f <- function(x) { x^3 + 4*x^2 -10 }

Let’s plot the function over the range [1,2]:

plot(x,f(x), type='l', lwd=1.5, main=expression(sin(x/2) * cos(x/4))) abline(h=0)

It seems obvious from the plot that the zero of f(x) in the range [1,2] is somewhere between x=1.3 and x=1.4.

This is made even more clear if we tabulate the x,y values over this range:

> xy <- cbind(x,f(x)) > xy x [1,] 1.0 -5.000 [2,] 1.1 -3.829 [3,] 1.2 -2.512 [4,] 1.3 -1.043 [5,] 1.4 0.584 [6,] 1.5 2.375 [7,] 1.6 4.336 [8,] 1.7 6.473 [9,] 1.8 8.792 [10,] 1.9 11.299 [11,] 2.0 14.000

Using the function defined earlier, we can run the root-finding algorithm:

> p <- newton(f, x0=1, N=10) > p [1] 1.454256 1.368917 1.365238 1.365230 1.365230 1.365230 1.365230

This returns 1.365230 as the root of f(x). Plotting the last value in the iteration:

> abline(v=p[length(p)])

In the iteration loop, we use a schoolbook-style numerical derivative:

$$ \frac{dy}{dx} = \frac{f(x+\epsilon)-f(x)}{\epsilon} $$

Its also worth noting that R has some support for symbolic derivatives, so we could use a symbolic function expression. A simple example of symbolic differentiation in R (note that R works with expression instances when using symbolic differentiation):

> e <- expression(sin(x/2)*cos(x/4)) > dydx <- D(e, "x") > dydx cos(x/2) * (1/2) * cos(x/4) - sin(x/2) * (sin(x/4) * (1/4))

In order for this to be useful, obviously we need to be able to evaluate the calculated derivative. We can do this using eval:

> z <- seq(-1,1,.1) > eval(dydx, list(x=z)) [1] 0.3954974 0.4146144 0.4320092 0.4475873 0.4612640 0.4729651 0.4826267 [8] 0.4901964 0.4956329 0.4989067 0.5000000 0.4989067 0.4956329 0.4901964 [15] 0.4826267 0.4729651 0.4612640 0.4475873 0.4320092 0.4146144 0.3954974

Note that we can bind the x parameter in the expression when calling `eval()`

.