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:
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));
#ifdef _DEBUG
std::cout << "Delta : " << delta << std::endl;
#endif
if (0.0 > delta)
delta *= -1.0;
return delta < Error;
}
};
struct Sine : public Funct
{
double Function( double x, int n = 0);
bool CalculateDiff(double fX);
};
struct Cosine : public Funct
{
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:
double Sine::Function( double x, int n)
{
return (sin(x) + n)/cos(x);
}
bool Sine::CalculateDiff(double fX)
{
return Funct::CalculateDiff<sin>(fX);
}
double Cosine::Function( double x, int n)
{
return (cos(x)+n)/-sin(x);
}
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 T , int N = 0>
class NR
{
double _x; double _fX; bool _success;
size_t _i;
NR(const NR& nr);
NR& operator=(const NR& nr);
inline void PrintCalculations();
inline void PrintResults();
inline void Function();
public:
static size_t MAXITER; 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);
template < class T , int N>
void NR<T, N>::Function()
{
T t;
_fX = _x - t.Function(_x,N);
_success = t.CalculateDiff(_fX);
}
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;
}
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