# An implementation of the Newton-Raphson algorithm in C/C++ and R

Today, we write a small piece of C/C++ code that implements the well-known Newton-Raphson algorithm (see, Mathworld). We also provide the R code.

Exercise: Find the unique root of the function $f(x) = e^{-x^2} - x^3 + 2$ using the Newton-Raphson method.

Notice that we choose a function that has a unique root (that is the value of x at f(0)) for simplicity.

We first draw a graph to see the general shape of the function. If we don’t get it right, that will tell us on the domain of the function. moreover, we will have a rough idea of where the root lies on the x-axis. That will then be helpful for the initial guess. In our case the guess is a value less than the actual root.

# Édouard Tallent @ TaGoMa.Tech
# September 2012
# Plots the function exp(-x^2) - x^3 + 2
x <- seq (-3, 3, 0.1)
y <- exp(-x^2) - x^3 + 2
plot(x, y, type = "l", col = "seagreen1",
lwd = 2,
lty = 5,
ylim = c(-2, 5),
ylab = expression(paste(f(x), " = ", e^-x^2 - x^3 + 2)))
abline(v = 0 , h = 0)

We then need the first derivative of f(x) = f'(x). Simply enough, $f\prime(x)=-2xe^{-x^2}-3x^2$. Then, we are ready for the implementation of the Newton-Raphson algorithm that will give us the root of f(x).

// Édouard Tallent @ TaGoMa.Tech
// September 2012
// This code implements the Newton-Raphson algorithm

#include<iostream>
#include<cmath>
using namespace std;

double newtonRaphson (double x)
{
return x - (exp(-pow(x, 2)) - pow(x, 3) + 2) / (-2 * exp(-pow(x, 2)) - 3 * pow(x, 2));
}

int main ()
{
double x = 0.5;         // the initial guess close and < to the actual root
double tol = 0.0001;    // 0.0001 is the error level we wish
double old;

do
{
old = x;
x = newtonRaphson(x);
}
while (abs(old - x) > tol);     // while loop because the
// number of iterations is not
// known beforehand
cout << "Root: " << x << endl;

return 0;
}

The compiler returns 1,29775 that is a value close to the actual root. And, we see it converged quickly towards the solution as only few iterations were required.

To check this result, one can write a piece of code in R as follows:

# Édouard Tallent @ TaGoMa.Tech
# September 2012
# This code implements the Newton-Raphson algorithm
newtonRaphson <- function (x){
x - ((exp(-x^2) - x^3 + 2) / (-2 * exp(-x^2) - 3 * x^2))
}

x <- 0.5	# initial guess
old <- 0	# sets old to 0 before the first loop
tol <- 0.0001

while (abs(old - x) > tol){
old <- x
x <- newtonRaphson(x)
}

print (paste("Root: ", x))

Finally, the result given in R is approximately the same as the one we obtained with the C/C++ snippet. We are done.