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:
[source lang=”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)])
}
[/source]
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:
[source lang="R"]f <- function(x) { x^3 + 4*x^2 -10 }[/source]
Let's plot the function over the range [1,2]:
[source lang=”R”]
plot(x,f(x), type=’l’, lwd=1.5, main=expression(sin(x/2) * cos(x/4)))
abline(h=0)
[/source]
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:
[source lang=”R”]
> 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
[/source]
Using the function defined earlier, we can run the root-finding algorithm:
[source lang=”R”]
> p <- newton(f, x0=1, N=10)
> p
[1] 1.454256 1.368917 1.365238 1.365230 1.365230 1.365230 1.365230
[/source]
This returns 1.365230 as the root of f(x). Plotting the last value in the iteration:
[source lang=”R”]
> abline(v=p[length(p)])
[/source]
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):
[source lang=”R”]
> 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))
[/source]
In order for this to be useful, obviously we need to be able to evaluate the calculated derivative. We can do this using eval:
[source lang=”R”]
> 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
[/source]
Note that we can bind the x parameter in the expression when calling eval().