12,552,576 members (27,634 online)
alternative version

20K views
8 bookmarked
Posted

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

, 26 Dec 2008 CPOL
 Rate this:
Templatized Newton-Raphson algorithm for SIN(X) + N and COS(X) + N where N is in [0, 2].

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

• Version 1.

## Share

programmer interested in cooking and dogs..

## You may also be interested in...

 Pro

 First Prev Next
 buggy invalid code. Bernard Gressing30-Dec-08 10:46 Bernard Gressing 30-Dec-08 10:46
 Re: buggy invalid code. jrivero7-Jan-09 11:19 jrivero 7-Jan-09 11:19
 My vote of 1 Ola Martin Lykkja23-Dec-08 11:23 Ola Martin Lykkja 23-Dec-08 11:23
 Re: My vote of 1 jrivero28-Dec-08 8:26 jrivero 28-Dec-08 8:26
 paste come code into article Adrian Pasik17-Dec-08 6:20 Adrian Pasik 17-Dec-08 6:20
 Re: paste come code into article jrivero28-Dec-08 8:27 jrivero 28-Dec-08 8:27
 Useful for small multidimensional problems also MJessick15-Dec-08 14:31 MJessick 15-Dec-08 14:31
 Re: Useful for small multidimensional problems also jrivero16-Dec-08 6:15 jrivero 16-Dec-08 6:15
 Last Visit: 31-Dec-99 18:00     Last Update: 24-Oct-16 11:41 Refresh 1