// Newton-Raphson numeric algorithm: find a solution for f( x ) = 0
// Note: func_x maybe a function object or a pointer to function
#ifndef NEWTONRAPHSON_H
#define NEWTONRAPHSON_H
#include <cmath>
template<class T> double NewtonRaphson(
T & func_x, double x0, double step_x,
double x_abs_err, int max_iter) {
enum {MAX_ITER_EXCEEDED, DIVISION_BY_ZERO, X_IN_TOLERANCE}
NRResult;
double y0 = func_x(x0);
double x1 = x0 + step_x;
double y1 = func_x(x1);
double x_jump;
NRResult = MAX_ITER_EXCEEDED;
for(int i = 0; i < max_iter; ++i) {
double x1_keep = x1;
double y1_keep = y1;
if ( y1 != y0 ) {
x_jump = y0 * (x1 - x0) / (y1 - y0);
x1 = x0 - x_jump;
}
else {
x_jump = x1 - x0;
NRResult = DIVISION_BY_ZERO;
}
double err_x = fabs( x_jump );
if ( err_x < x_abs_err ) NRResult = X_IN_TOLERANCE;
if ( NRResult != MAX_ITER_EXCEEDED ) break;
y1 = func_x(x1);
x0 = x1_keep;
y0 = y1_keep;
}
double nr_value;
switch ( NRResult ) {
case X_IN_TOLERANCE:
nr_value = x1;
break;
default:
nr_value = 0.0;
}
return nr_value;
}
#endif // NEWTONRAPHSON_H
End of Listing