C#，数值计算——抛物线插值与Brent方法（Parabolic Interpolation and Brent‘s Method）的计算方法与源程序

 

using System;

namespace Legalsoft.Truffer
{
    /// <summary>
    /// 抛物线插值与Brent方法
    /// Parabolic Interpolation and Brent's Method
    /// </summary>
    public class Brent : Bracketmethod
    {
        public double xmin { get; set; }
        public double fmin { get; set; }
        public double tol { get; set; }

        public Brent(double toll = 3.0e-8)
        {
            this.tol = toll;
        }

        public double minimize(UniVarRealValueFun func)
        {
            const int ITMAX = 100;
            const double CGOLD = 0.3819660;
            double ZEPS = float.Epsilon * 0.001;
            double d = 0.0;
            double e = 0.0;
            double a = (ax < cx ? ax : cx);
            double b = (ax > cx ? ax : cx);
            double x = bx;
            double w = bx;
            double v = bx;
            double fx = func.funk(x);
            double fw = fx;
            double fv = fx;
            for (int iter = 0; iter < ITMAX; iter++)
            {
                double xm = 0.5 * (a + b);
                double tol1 = tol * Math.Abs(x) + ZEPS;
                double tol2 = 2.0 * (tol1);
                if (Math.Abs(x - xm) <= (tol2 - 0.5 * (b - a)))
                {
                    fmin = fx;
                    return xmin = x;
                }
                if (Math.Abs(e) > tol1)
                {
                    double r = (x - w) * (fx - fv);
                    double q = (x - v) * (fx - fw);
                    double p = (x - v) * q - (x - w) * r;
                    q = 2.0 * (q - r);
                    if (q > 0.0)
                    {
                        p = -p;
                    }
                    q = Math.Abs(q);
                    double etemp = e;
                    e = d;
                    if (Math.Abs(p) >= Math.Abs(0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x))
                    {
                        d = CGOLD * (e = (x >= xm ? a - x : b - x));
                    }
                    else
                    {
                        d = p / q;
                        double uu = x + d;
                        if (uu - a < tol2 || b - uu < tol2)
                        {
                            d = Globals.SIGN(tol1, xm - x);
                        }
                    }
                }
                else
                {
                    d = CGOLD * (e = (x >= xm ? a - x : b - x));
                }
                double u = (Math.Abs(d) >= tol1 ? x + d : x + Globals.SIGN(tol1, d));
                double fu = func.funk(u);
                if (fu <= fx)
                {
                    if (u >= x)
                    {
                        a = x;
                    }
                    else
                    {
                        b = x;
                    }
                    shft3(ref v, ref w, ref x, u);
                    shft3(ref fv, ref fw, ref fx, fu);
                }
                else
                {
                    if (u < x)
                    {
                        a = u;
                    }
                    else
                    {
                        b = u;
                    }
                    if (fu <= fw || w == x)
                    {
                        v = w;
                        w = u;
                        fv = fw;
                        fw = fu;
                    }
                    else if (fu <= fv || v == x || v == w)
                    {
                        v = u;
                        fv = fu;
                    }
                }
            }
            throw new Exception("Too many iterations in brent");
        }
    }
}