The modified false position method is fairly straightforward to implement. The method is described in most numerical methods texts in easily translatable pseudo code. The method requires evaluating the function during the solution process, requiring in a naïve implementation hard coding of the function being solved.
Function pointers can be used to pass a general function to the method, so that one implementation can be used to solve various functions. In particular, member function pointers can be used to root solve a function for one of many independent variables. Consider for example, f(x,a), where you would like to find x such that f(x,a) = 0 for a particular value of a.
One solution is to define a function g(x) = f(x,a) in a class, using a from a class member, then use the solver to solve g(x)=0. The class is used to contain the fixed independent variables so that a function pointer with one argument can be passed to the solver.
Being somewhat new to C++ myself, getting this to work took some research into member function pointers. My implementation is included below. To use it, one would, in a class, define a function which has an argument of double and returns a double using member variables to evaluate the function given a particular x. Then in the class create a member function pointer to the function being solved.
For example one could use the solver in a class "TranscendentalSolver" with a double fTranscendental(double x) member function like so:
double (TranscendentalSolver::*f)(double) = &TranscendentalSolver::fTranscendental;
// use a modified false position solver to solve the transcendental
b = ModFalsePos(f, *this, xLeft, xRight, 1E-15, 1E-15, 100, verbose);
My implementation of the modified false position method is listed below. Note: I modified it to use bisection in the case when the new root is out of bounds (an issue I had solving atan).
#ifndef __MODIFIEDFALSEPOSITION_H
#define __MODIFIEDFALSEPOSITION_H
#include
<cstdio>
#include
<math.h>
//
// Modified False Position Method, for details
// see p 129 in Numerical Methods by Chapra.
// f - function pointer to the function to find a root
// xl - lower x bound
// xu - upper x bound
// es - relative stopping error
// et - stopping error (abs(f(xr)) < et)
// imax - maximum iterations
// verbose - whether to print information to cout
//
// This is a templated Function so that a member function pointer
// may be passed. This is very useful when the function has
// several arguments, but only one is being varied. Then the
// caller may define a member function, using class members to
// find the other parameters.
template<typename F>
double ModFalsePos(double (F::*f)(double), F & obj,
double xl, double xu, double es, double et, int imax,
bool verbose)
{
// the current iteration
int iter = 0;
// lower and upper iteration counts,
// used to detect when one of the bounds
// is 'stuck'.
int il = 0;
int iu = 0;
// the current root
double xr = 0;
// the current error estimate
double ea = 0;
// the lower function value
double fl = (obj.*f)(xl);
// the upper function value
double fu = (obj.*f)(xu);
if (verbose)
{
printf("xl = %12.4e, xu = %12.4e\n", xl, xu);
printf("%6s%12s%12s%12s%12s%12s\n",
"Iter.", "x_l", "x_u", "x_r", "f(x_r)", "ea");
}
while(true)
{
// new root estimate and function value
double xrold = xr;
xr = xu - fu * (xl - xu) / (fl - fu);
// if out of bounds - false position
// has failed, use bisection
if ( xr > xu || xr < xl )
{
xr = (xl+xu)/2;
}
double fr = (obj.*f)(xr);
// iteration counter increment
iter += 1;
// error estimate (relative)
if (xr != 0 )
ea = std::abs((xr-xrold)/xr)*100;
// test signs
double test = fl*fr;
if (test < 0)
{
xu = xr;
fu = (obj.*f)(xu);
iu = 0;
il += 1;
if (il >= 2)
fl = fl / 2;
}
else
if (test > 0)
{
xl = xr;
fl = (obj.*f)(xl);
il = 0;
iu += 1;
if (iu >= 2)
fu = fu /2;
}
else
{
ea = 0;
}
if(verbose)
{
printf("%6i%12.4e%12.4e%12.4e%12.4e%12.4e\n",
iter, xl, xu, xr, fr, ea);
}
if (ea < es || std::abs(fr) < et || iter >= imax)
break;
}
return xr;
}
#endif