icedb::zeros Namespace Reference

Functions
double	findzero (double a, double b, const std::function< double(double) > &evaltarget)
	Zero-finding implementation - Brent's method. More...
template<class T , class U >
U	secantMethod (const T &f, U guess_a, U guess_b, double eps=0.000001, size_t maxIter=50)

Function Documentation

findzero()

double icedb::zeros::findzero	(	double	a,
		double	b,
		const std::function< double(double) > &	f
	)

Zero-finding implementation - Brent's method.

Definition at line 7 of file zeros.cpp.

References ICEDB_throw, and icedb::error::xBadInput.

        {
            using namespace std;
            /* Define the basic testing variables.
               eval func is the function's evaluator.
               It typically links to the polynomial eval func.
               The rest of the vars are used in Brent's method.
               */
            //evalfunction *f = evaltarget;
            //double (*f)(double) = (double (*)(double)) evalfunc;
            double convint = 1E-7;
            double fcint = 1E-13;
            double fa, fb, fc, fd, fs = 1.0;
            unsigned int i=0;
            fa = f(a);
            fb = f(b);

            // Bad starting points if fa*fb > 0
            if (fa*fb > 0) 
                ICEDB_throw(icedb::error::error_types::xBadInput)
                    .add<std::string>("Reason",
                        "Bad selection of a,b. f(a) and f(b) need to have opposite signs.")
                .add<double>("a",a).add<double>("b",b)
                .add<double>("fa",fa).add<double>("fb",fb);

            if ( abs(fa) < abs(fb) )
            {
                // Swap the two here. Junk only declared in block.
                double junk = b;
                b = a;
                a = junk;
                fa = f(a);
                fb = f(b);
            }

            double c = a; // A testing point
            double s = 0, d = 0;  // Initialize these to stop VS whining
            fc = f(c);
            bool mflag = true; // Flag for method stuff
            while (fb != 0 || fs != 0 || abs(b-a) > convint)
            {
                if (fa != fc && fb != fc)
                {
                    // Perform inverse quadratic interpolation
                    double i,j,k;
                    i = a * fb * fc / ( (fa - fb)*(fa - fc));
                    j = b * fa * fc / ( (fb - fa)*(fb - fc));
                    k = c * fa * fb / ( (fc - fa)*(fc - fb));
                    s = i + j + k;
                } else {
                    // Use Secant Rule
                    s = b - fb * (b - a) / ( fb - fa );
                }

                if (
                        ( (3.0*a+b)/4.0 < s && s < b) == false ||
                        (mflag && abs(s-b) >= (abs(b-c)/2.0) ) ||
                        (!mflag && abs(s-b) >= (abs(c-d)/2.0) ) ||
                        (mflag && abs(b-c) < convint) ||
                        (!mflag && abs(c-d) < convint)
                        )
                {
                    // Bisection Method
                    mflag = true;
                    s = 0.5 * (a + b);
                } else {
                    mflag = false;
                }

                fs = f(s);

                d = c;
                c = b;

                fd = f(d);
                fc = f(c);

                if (fa * fs < 0)
                {
                    b = s;
                    fb = f(b);
                } else {
                    a = s;
                    fa = f(a);
                }


                if ( abs(fa) < abs(fb) )
                {
                    // Do a swap again
                    double junk = b;
                    b = a;
                    a = junk;
                    fa = f(a);
                    fb = f(b);
                }
                i++;
                //if (i > 100) break;
                if ( abs(b-a) < convint && abs(fb) < fcint && abs(fs) < fcint ) break;
            }
            return s;
        }

secantMethod()

template<class T , class U >

U icedb::zeros::secantMethod	(	const T &	f,
		U	guess_a,
		U	guess_b,
		double	eps = 0.000001,
		size_t	maxIter = 50
	)

Definition at line 17 of file zeros.hpp.

References ICEDB_throw, and icedb::error::xModelOutOfRange.

Referenced by icedb::refract::bruggeman(), icedb::refract::guessTemp(), and icedb::refract::sihvola().

        {
            // Secant method is defined by recurrance
            // xn = x_(n-1) - f(x_(n-1)) * (x_(n-1) - x_(n-2)) / (f(x_(n-1)) - f(x_(n-2)))
            using namespace std;
            U zero;
            U xn = guess_a, xn1 = guess_b, xn2;
            U fxn1, fxn2;
            size_t i=0;
            do {
                xn2 = xn1;
                xn1 = xn;

                fxn1 = f(xn1);
                fxn2 = f(xn2);

                xn = xn1 - fxn1 * (xn1 - xn2) / (fxn1 - fxn2);
            } while ( (abs(xn-xn1) > eps) && (i++ < maxIter));

            if (i >= maxIter) ICEDB_throw(icedb::error::error_types::xModelOutOfRange)
                .add<size_t>("maxIter", maxIter)
                .add<double>("eps", eps)
                .template add<U>("guess_a", guess_a)
                .template add<U>("guess_b", guess_b)
                .template add<U>("xn",xn)
                .template add<U>("xn1",xn1)
                .template add<U>("xn2",xn2)
                .template add<U>("fxn1",fxn1)
                .template add<U>("fxn2",fxn2);
            zero = xn;
            return zero;
        }

Here is the caller graph for this function: