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  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, . 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.

About these ads