Templatized Newton-Raphson algorithm for SIN(X) + N and COS(X) + N where N is in [0, 2]

• Download source code - 14.3 KB

Introduction

Hello, and thanks for stopping by. If you are here, you know what the Newton-Raphson algorithm does. In short, it finds the solution to f(x) of an equation of the form:

f(x) = ax + c

by iteratively using derivatives until there is a point where the result reaches 0 on the y-axis. There is an initial value for x to start the iterations (the seed). In this case, the functions are SIN and COS, just to make it more interesting (and simple, I think). If you want to learn more about Newton-Raphson, there is a ton of info on the net. Please search for it using your favorite Google search engine! ;).

Background

This code was given to me as a C++ Templates test, as well as to show some understanding of basic Calculus, and I thought that maybe someone might find it useful.

After some research, I found that the choice of 'seed' to start the iterations must be sufficiently close to the desired root. For non-trigonometric functions, the choice of a 'seed' close to a value that will result in f(x) = 0 is advisable.

However, for trigonometric functions, their cyclical nature will make them always vary from 0 to 2PI. In other words, the result of the iteration (the seed for the next iteration) has to be reduced to a value on the circumference -in radians- to figure out if we are on or close enough to the root.

Finding the root for the function

In order to know if the result from the iteration is approximating to 0, we apply the corresponding trigonometric function (either SIN or COS) to the result of the iteration to figure out if the value is 0 or close to 0. This value is stored in the variable delta.

If delta is within an 'acceptable' value close to 0, the iterations can be stopped. To achieve this, it is needed to set an acceptable minimum error threshold. 1E-5 was chosen as an acceptable error threshold to stop the iterations.

Therefore, if delta is less than error, the iterations are stopped and the result is displayed, avoiding infinite iterations for functions that do not converge to 0. A maximum number of iterations were set (given by the static variable size_t MAXITER). Otherwise, for non-convergent graphs, the program could enter in an infinite loop without ever finding the result.

Using the code

This is a sample of the code in main():

switch(n) 
{ 
    default:
    case 0:
    {
        NR <Sine> ns(x);
        ns.Process();
        break;
    }
    case 1:
    {
        NR <Sine,1> ns(x);
        ns.Process();
        break;
    }
    case 2:
    {
        NR <Sine,2> ns(x);
        ns.Process(); 
        break;
    }
}

As you can see, the template to be instantiated needs to be known at runtime. In other words, if you try to do this:

NR <Sine, n> ns(x);

...the compiler will yelp. If anybody has an idea of how to make this happen, I am open to suggestions and learning.

The code that performs the calculations is:

// interface class to homogenize the 
// behavior of our  function classes
struct Funct
{
    virtual double Function( double x, int n = 0) = 0;
    virtual bool CalculateDiff(double fX) = 0;
    template < double f(double) >
    bool CalculateDiff(double fX)
    {
        double delta(f(fX));
        // f can be SIN or COS functions from math.h header
#ifdef _DEBUG
        std::cout << "Delta : " << delta << std::endl;
#endif
        if (0.0 > delta)
            delta *= -1.0;
        return delta < Error;
    }

};
// class Sine: performs calculations for SIN(X) function
struct Sine : public Funct
{
    // Function sin
    double Function( double x, int n = 0);
    bool CalculateDiff(double fX);
};
// class Cosine: performs calculations for COS(X) function
struct Cosine : public Funct
{
    // Function sin
    double Function( double x, int n = 0);
    bool CalculateDiff(double fX);
};

As you can see, we utilize a class for each function, using the same interface. To enforce this, we derive the func classes from an abstract class (or protocol class, or interface). The implementation of our class functions is as follows:

// Newton-Raphson for sin
double Sine::Function( double x, int n)
{
    return (sin(x) + n)/cos(x);
}
////////////////////////////////////////////////////////
// Calculates differences for sine
bool Sine::CalculateDiff(double fX)
{
    return Funct::CalculateDiff<sin>(fX);
}
////////////////////////////////////////////////////////
// Newton-Raphson for cosine
double Cosine::Function( double x, int n)
{
    return (cos(x)+n)/-sin(x);
}
////////////////////////////////////////////////////////
// Calculate differences for Cosine
bool Cosine::CalculateDiff(double fX)
{
    return Funct::CalculateDiff<cos>(fX);
}

And finally, the wrapper that allows us to utilize all of the above magic is:

// template class to work as a wrapper
// for the NewtonRapshon algorithm
template < class T , int N = 0>
class NR
{
    double _x; // seed to start iterations from. 
    double _fX; // function to approximate to 0. 
    // to mark if the iterations where successful
    // after the max number of iterations
    bool _success;
    size_t _i; // iterator

    // Copy constructor and assignment made
    // private to disallow copy and assignment
    NR(const NR& nr);
    NR& operator=(const NR& nr);
    // methods
    inline void PrintCalculations();
    inline void PrintResults();
    inline void Function();
public:
    static size_t MAXITER; // Maximum number of iterations
    explicit NR(double x) : _x(x), _fX(0.0) , _success(false), _i(0)
    {
    }
    void Process();
};

template  < class T, int N > size_t NR < T, N > ::MAXITER(50);

// execute function
template  < class T , int N>
void NR<T, N>::Function()
{
    T t;
    _fX = _x - t.Function(_x,N);
    _success = t.CalculateDiff(_fX);
}
// Print results
template  < class T, int N >
void NR<T, N>::PrintResults()
{
    using namespace std;
    if (_success)
        cout << "The root is: [" << _fX << "]. Result reached in " 
             << _i+1 << " Iterations." << endl;
    else
        cout << "** Result NOT reached in " << _i 
             << " Iterations." << endl;
}

template  < class T, int N >
void NR<T,N>::PrintCalculations()
{
    std::cout << _i+1 << "] x: " << _x << " - (" << sin(_x) 
              << "/" << cos(_x) << ") = f(x) : " << _fX << std::endl;
}
// Executes algorithm
// Iterates until: 
// a) result is found 
// or 
//b) MAXITER has been reached without a result found
template  < class T, int N >
void NR<T,N>::Process()
{
    for (; _i < NR::MAXITER ; _i++)
    {
        Function();
#ifdef _DEBUG
        PrintCalculations();
#endif
        if (_success)
        {
            break;
        }
        _x = _fX;
    }
    PrintResults();
}

The code included has been written and tested using MS Visual Studio 2005 and the GNU C++ compiler 3.4.4. To compile on Linux, use:

c++ -D_DEBUG  -Wall -pedantic -o NewRaph NewRaph.cpp StrHelper.cpp

…if you want to see the iterations and the deltas, or

c++ -Wall -pedantic -o NewRaph NewRaph.cpp StrHelper.cpp 

…if you don't.

Sample executions

Below is a sample of executing the program without parameters:

NewRaph
Usage: NewRaph SEED {S|C|B} [SHIFT] where:
SEED: Integer. Seed to start the iterations from.
{S|C|B}: S to calculate SIN(x), C for COS(X), B for Both
[SHIFT]: n value to establish the function f(x) + n. n varies from 0 to 2

Points of interest

For SIN(X) and COS(X), the algorithm very quickly converges to f(x) = 0. For SIN(X) + N and COS(X) + N, the algorithm converges towards -1.

There might be a better way to implement this using templates (using template metaprogramming or something like that), but I am not aware of it. Feedback is welcomed.

History

• Version 1.