CODEWITHSUNDEEP
c++algorithm
user3215609
user3215609

Reputation: 535

Newton-Raphson algorithm: find all roots include negative

How can I use Newton's method to find all roots of polynomial, not only unique?

Link to code: http://quantcorner.wordpress.com/2012/09/14/an-implementation-of-the-newton-raphson-algorithm-in-c-cpp-and-r/

As an example I have this equation: x^2-12x+34=0, when I use this formula I get only one root 4.5857, how can I get second root - 7.4142?

Upvotes: 1

Views: 10817

Answers (1)

user3116431
user3116431

Reputation: 283

Restart your Newton Raphson from a different starting position.

I have pasted here some code I have implemented in developing a library of financial instrument pricing. See here my code for Newton Raphson. You see as one of the input parameters a start position double start. You can start from different positions on a grid e.g. and compare solutions with what you have already found as solution.

#include "math.h"
#include "ExceptionClass.h"

template <class T, double (T::*Value)(const double) const,
                   double (T::*Derivative)(const double) const>
        double NewtonRaphson(double Target, double start, double Tolerance,
                             size_t max_count, const T& Object)
    {
    size_t count = 0;
    double diff = Target-(Object.*Value)(start);
    double x = start;
    double derivative = (Object.*Derivative)(start);

    do{
        count++;
        if(fabs(derivative) < 1.E-6)
            throw DivideByZeroException("Dividing by a derivative smaller than: 1.E-6!");

        x = x - diff/(-derivative);

        diff = Target-(Object.*Value)(x);
        derivative = (Object.*Derivative)(x);
    } while((fabs(diff)>Tolerance)&&(count < max_count));

    if(count >= max_count)
        throw("Newton-Raphson did not converge in the defined number of steps!");
    else return x;
    }

T is a class where you have defined functions to evaluate your (quadratic) equation (here referred to by double (T::*Value)(const double) const in the template) and the derivative of this equation (here referred to by double (T::*Derivative)(const double) const in the template). Note I have included an exception class as Newton Raphson has some issues.

Use the bisection algorithm for a more stable algorithm. In practice you can use numbers smaller than 1.E-6.

Value should be here a pointer to a function that evaluates your quadratic function, the target should be set to 0 and the derivative should be a pointer to a function that evaluates the derivative of your function.

In your case my template code can be simplified: Replace the code diff = Target-(Object.*Value)(x); with diff = (Object.*Value)(x); Replace x = x - diff/(-derivative); with x = x - diff/(derivative); This code will be more easily recognized as the Newton Raphson algorithm. If you want to increase the speed of convergence use Halley's (yes the guy from the comet) algorithm or Householder's algorithm.

See numerical recipes in c++ for alternative implementation. Numerical recipes in C++ is the book to go for to address these kind of questions.

Upvotes: 2

Related Questions