posts - 81, comments - 257, trackbacks - 0

Modified False Position Method in C++ Accepting a Function Pointer


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


Print | posted on Wednesday, November 14, 2012 11:32 AM | Filed Under [ Technical Computing ]

Feedback

No comments posted yet.

Post Comment

Title  
Name  
Email
Url
Comment   
Please add 6 and 1 and type the answer here:

Powered by:
Powered By Subtext Powered By ASP.NET