Hi. I find the NR books very useful. I have used some of the concepts, but implemented the functions where another function has to be passed in (e.g. the solvers) in a different way. I like to make a C++ class which includes, as member functions, both the solver, and the function being solved. Then, in each application, I define another derived class, where I supply a different function to be solved, overriding the function in the base class.

I find this useful, and easy to program. Perhaps
you might also find this technique useful.

This sample program I have enclosed prints out a table of square roots, just to demonstrate the solver.

Morris Pearl

# define TEST_MAIN 1

/* implements solving non-linear equation, using Ridder's formula */
/* from the Numerical Recipies in C book (2nd edition) */

/* To use this class, derive your own class from this one, and */
/* make a function named f, which is to be solved. See an example */
/* in the solver.cc file which creates square root tables */


class solver
{
public:

virtual double f(double x) = 0;

double lo;
double hi;
double g;
double tolx;
double toly;

double function_tolerance; // function values < this are considered zero

int ridder_solve(double *result);

solver() /* constructor */
{
tolx = 0.01;
toly = 1.0;
function_tolerance = 0;
}
};


// Module Name : solver.cc

# ifdef TEST_MAIN
# include <stdio.h>

main()
{
double answer;
int j,k;

class sqrt_solver : public solver
/* this is derived from the "solver" class, and redefines */
/* just the function (f) that we are trying to solve */
{
public:

double param;

double f(double x)
{
return (x * x - param);
}
} sq;


for (j = 2; j < 20; j++)
{
sq.param = j;

sq.lo = 0;
sq.g = j * 0.5;
sq.hi = j;

k = sq.ridder_solve(&answer);
printf("%3d %8.4f %8.4f %4d\n",j,answer,answer * answer, k);
}
}

# endif

# define abs(z) (((z) < 0) ? (-(z)) : (z))
# define min(x,y) ( ((x) < (y)) ? (x) : (y))
# define max(x,y) ( ((x) > (y)) ? (x) : (y))

# define EPSILON (1.0e-14)
# include <math.h>


int solver::ridder_solve(double *result)
{
double fhi;

double flo = f(lo);
double fg = f(g);
int count = 2;

/* lo will be a value such that f(lo) < 0 */
/* hi will be a value such that f(hi) > 0 */
/* flo is f(lo) fhi is f(hi) */
/* the tolerance criteria means that we will find some */
/* numbers such that either: */
/* f(lo) < 0, f(hi) > 0, |hi-lo| < tolx, and at least one of f(lo) or f(hi) < toly. */
/* (in which case we return a linear interpolation between the points lo,f(lo) and hi,f(hi) */
/* OR */
/* f(lo) < 0, f(hi) > 0, g is between lo and hi, and |f(g)| < function_tolerance */
/* in which case we will then return g */

if ((flo < 0) ^ (fg < 0))
{
hi = g;
fhi = fg;
}
else
{
fhi = f(hi);
count++;
}

if (flo > 0)
{
double temp = lo;
lo = hi;
hi = temp;
temp = flo;
flo = fhi;
fhi = temp;
}

if (flo > 0) return -1;
if (fhi < 0) return -1;
if ((toly < EPSILON) || (tolx < EPSILON)) return -3;

if (fg < -function_tolerance)
{
flo = fg;
lo = g;
}
else if (fg > function_tolerance)
{
fhi = fg;
hi = g;
}
else
{
*result = g;
return count;
}

while (count < 100)
{
double prev_g = g;
double m = (lo + hi) * 0.5;
double fm = f(m); /* bisection guess */
count++;

if ((abs(hi - lo) < (2.0 * tolx)) && (abs(fm) < toly))
{
/* If the remaining region is less then double the tolerance */
/* half of the region has to be less than the tolerance */
/* therefore, we can just fall through to the final */
/* linear interpolation step */
g = m;
fg = fm;
}
else
{
double s = sqrt(fm * fm - flo * fhi);
if (s < function_tolerance)
{
*result = m;
return count;
}

g = m + (lo - m) * fm / s;
// Ridders' formula, Numerical Recipes in C, 2nd edition, page 358

if (abs(prev_g - g) < 4 * tolx)
{
// If guess is going nowhere, move it a tiny bit towards m.
// This is to handle the case where the root is at one end
// of the interval, but we keep guessing on the same side

if (g + tolx < m)
g += tolx * 0.5;
else if (g - tolx > m)
g -= tolx * 0.5;
}

fg = f(g);
count++;
}

if ((fm > function_tolerance) && (fg < -function_tolerance))
{
lo = g;
flo = fg;
hi = m;
fhi = fm;
}
else if ((fm < -function_tolerance) && (fg > function_tolerance))
{
lo = m;
flo = fm;
hi = g;
fhi = fg;
}
else if ((fm < -function_tolerance) && (fg < -function_tolerance))
{
lo = g;
flo = fg;
}
else if ((fm > function_tolerance) && (fg > function_tolerance))
{
hi = g;
fhi = fg;
}
else
{
/* we found the answer, */
/* use whichever guess has a lower functional value */
*result = (abs(fm) < abs(fg)) ? m : g;
return count;
}
if ((abs(hi - lo) <= tolx) &&
((abs(flo) < toly) || (abs(fhi) < toly))) /* if tolerance is met */
{
/* do one last linear interpolation */
*result = lo + (hi - lo) * flo / (flo - fhi);
return count;
}