Click here to Skip to main content
12,348,674 members (53,394 online)
Click here to Skip to main content

Stats

48K views
2.2K downloads
48 bookmarked
Posted

Three Methods for Root-finding in C#

, 26 Jun 2014 BSD
Three numerical algorithms for solving equations, each implemented in C#
RootFindingDemo
RootFindingDemo.suo
Properties
RootFindingDemo.csproj.user
using System;

namespace DemoApplication
{
    public delegate double FunctionOfOneVariable(double x);

    class RootFinding
    {
        const int maxIterations = 50;

        public static double Bisect
        (
            FunctionOfOneVariable f,
            double left,
            double right,
            double tolerance = 1e-6,
            double target = 0.0
        )
        {
            // extra info that callers may not always want
            int iterationsUsed;
            double errorEstimate;

            return Bisect(f, left, right, tolerance, target, out iterationsUsed, out errorEstimate);
        }

        public static double Bisect
        (
            FunctionOfOneVariable f,
            double left,
            double right,
            double tolerance,
            double target,
            out int iterationsUsed,
            out double errorEstimate
        )
        {
            if (tolerance <= 0.0)
            {
                string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
                throw new ArgumentOutOfRangeException(msg);
            }

            iterationsUsed = 0;
            errorEstimate = double.MaxValue;

            // Standardize the problem.  To solve f(x) = target,
            // solve g(x) = 0 where g(x) = f(x) - target.
            FunctionOfOneVariable g = delegate(double x) { return f(x) - target; };


            double g_left = g(left);  // evaluation of f at left end of interval
            double g_right = g(right);
            double mid;
            double g_mid;
            if (g_left * g_right >= 0.0)
            {
                string str = "Invalid starting bracket. Function must be above target on one end and below target on other end.";
                string msg = string.Format("{0} Target: {1}. f(left) = {2}. f(right) = {3}", str, g_left + target, g_right + target);
                throw new ArgumentException(msg);
            }

            double intervalWidth = right - left;

            for
            (
                iterationsUsed = 0;
                iterationsUsed < maxIterations && intervalWidth > tolerance;
                iterationsUsed++
            )
            {
                intervalWidth *= 0.5;
                mid = left + intervalWidth;

                if ((g_mid = g(mid)) == 0.0)
                {
                    errorEstimate = 0.0;
                    return mid;
                }
                if (g_left * g_mid < 0.0)           // g changes sign in (left, mid)    
                    g_right = g(right = mid);
                else                            // g changes sign in (mid, right)
                    g_left = g(left = mid);
            }
            errorEstimate = right - left;
            return left;
        }

        public static double Brent
        (
            FunctionOfOneVariable f,
            double left,
            double right,
            double tolerance = 1e-6,
            double target = 0.0
        )
        {
            // extra info that callers may not always want
            int iterationsUsed;
            double errorEstimate;

            return Brent(f, left, right, tolerance, target, out iterationsUsed, out errorEstimate);
        }

        public static double Brent
        (
            FunctionOfOneVariable g,
            double left,
            double right,
            double tolerance,
            double target,
            out int iterationsUsed,
            out double errorEstimate
        )
        {
            if (tolerance <= 0.0)
            {
                string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
                throw new ArgumentOutOfRangeException(msg);
            }

            errorEstimate = double.MaxValue;

            // Standardize the problem.  To solve g(x) = target,
            // solve f(x) = 0 where f(x) = g(x) - target.
            FunctionOfOneVariable f = delegate(double x) { return g(x) - target; };

            // Implementation and notation based on Chapter 4 in
            // "Algorithms for Minimization without Derivatives"
            // by Richard Brent.

            double c, d, e, fa, fb, fc, tol, m, p, q, r, s;

            // set up aliases to match Brent's notation
            double a = left; double b = right; double t = tolerance;
            iterationsUsed = 0;

            fa = f(a);
            fb = f(b);

            if (fa * fb > 0.0)
            {
                string str = "Invalid starting bracket. Function must be above target on one end and below target on other end.";
                string msg = string.Format("{0} Target: {1}. f(left) = {2}. f(right) = {3}", str, target, fa + target, fb + target);
                throw new ArgumentException(msg);
            }

        label_int:
            c = a; fc = fa; d = e = b - a;
        label_ext:
            if (Math.Abs(fc) < Math.Abs(fb))
            {
                a  =  b;  b =  c;  c =  a;
                fa = fb; fb = fc; fc = fa;
            }

            iterationsUsed++;

            tol = 2.0 * t * Math.Abs(b) + t;
            errorEstimate = m = 0.5 * (c - b);
            if (Math.Abs(m) > tol && fb != 0.0) // exact comparison with 0 is OK here
            {
                // See if bisection is forced
                if (Math.Abs(e) < tol || Math.Abs(fa) <= Math.Abs(fb))
                {
                    d = e = m;
                }
                else
                {
                    s = fb / fa;
                    if (a == c)
                    {
                        // linear interpolation
                        p = 2.0 * m * s; q = 1.0 - s;
                    }
                    else
                    {
                        // Inverse quadratic interpolation
                        q = fa / fc; r = fb / fc;
                        p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));
                        q = (q - 1.0) * (r - 1.0) * (s - 1.0);
                    }
                    if (p > 0.0)
                        q = -q;
                    else
                        p = -p;
                    s = e; e = d;
                    if (2.0 * p < 3.0 * m * q - Math.Abs(tol * q) && p < Math.Abs(0.5 * s * q))
                        d = p / q;
                    else
                        d = e = m;
                }
                a = b; fa = fb;
                if (Math.Abs(d) > tol)
                    b += d;
                else if (m > 0.0)
                    b += tol;
                else
                    b -= tol;
                if (iterationsUsed == maxIterations)
                    return b;

                fb = f(b);
                if ((fb > 0.0 && fc > 0.0) || (fb <= 0.0 && fc <= 0.0))
                    goto label_int;
                else
                    goto label_ext;
            }
            else
                return b;
        }

        public static double Newton
        (
            FunctionOfOneVariable f,
            FunctionOfOneVariable fprime,
            double guess,
            double tolerance = 1e-6,
            double target = 0.0
        )
        {
            // extra info that callers may not always want
            int iterationsUsed;
            double errorEstimate;

            return Newton(f, fprime, guess, tolerance, target, out iterationsUsed, out errorEstimate);
        }

        public static double Newton
        (
            FunctionOfOneVariable f,
            FunctionOfOneVariable fprime,
            double guess,
            double tolerance,
            double target,
            out int iterationsUsed,
            out double errorEstimate
        )
        {
            if (tolerance <= 0)
            {
                string msg = string.Format("Tolerance must be positive. Recieved {0}.", tolerance);
                throw new ArgumentOutOfRangeException(msg);
            }

            iterationsUsed = 0;
            errorEstimate = double.MaxValue;

            // Standardize the problem.  To solve f(x) = target,
            // solve g(x) = 0 where g(x) = f(x) - target.
            // Note that f(x) and g(x) have the same derivative.
            FunctionOfOneVariable g = delegate(double x) { return f(x) - target; };

            double oldX, newX = guess;
            for
            (
                iterationsUsed = 0;
                iterationsUsed < maxIterations && errorEstimate > tolerance;
                iterationsUsed++
            )
            {
                oldX = newX;
                double gx = g(oldX);
                double gprimex = fprime(oldX);
                double absgprimex = Math.Abs(gprimex);
                if (absgprimex > 1.0 || Math.Abs(gx) < double.MaxValue * absgprimex)
                {
                    // The division will not overflow
                    newX = oldX - gx / gprimex;
                    errorEstimate = Math.Abs(newX - oldX);
                }
                else
                {
                    newX = oldX;
                    errorEstimate = double.MaxValue;
                    break;
                }
            }
            return newX;
        }
    }
}

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

License

This article, along with any associated source code and files, is licensed under The BSD License

Share

About the Author

John D. Cook
Singular Value Consulting
United States United States
I am an independent consultant in software development and applied mathematics. I help companies learn from their data to make better decisions.

Check out my blog or send me a note.

 


You may also be interested in...

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.160621.1 | Last Updated 26 Jun 2014
Article Copyright 2010 by John D. Cook
Everything else Copyright © CodeProject, 1999-2016
Layout: fixed | fluid