brent.cpp

Visual C++ 常用数值算法集 源代码

double brent(double ax,double bx,double cx,double tol,double& xmin)
{
    int done,iter,itmax = 100;
    double d,fu,r,q,p,xm,tol1,tol2,a,b,cgold = 0.381966;
    double u,etemp,dum,v,w,x,e,fx,fv1,fw,zeps = 0.0000000001;
    a = ax;
    if (cx < ax)
	{
		a = cx;
	}
    b = ax;
    if (cx > ax)
	{
		b = cx;
	}
    v = bx;
    w = v;
    x = v;
    e = 0.0;
    fx = func(x);
    fv1 = fx;
    fw = fx;
    for (iter = 1; iter<=itmax; iter++)
	{
        xm = 0.5 * (a + b);
        tol1 = tol * fabs(x) + zeps;
        tol2 = 2.0 * tol1;
        if (fabs(x - xm) <= tol2 - 0.5 * (b - a))
		{
			break;
		}
        done = -1;
        if (fabs(e) > tol1)
		{
            r = (x - w) * (fx - fv1);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w) * r;
            q = 2.0 * (q - r);
            if (q > 0.0)
			{
				p = -p;
			}
            q = fabs(q);
            etemp = e;
            e = d;
            dum = fabs(0.5 * q * etemp);
            if (fabs(p) < dum && p > q * (a - x) && p < q * (b - x))
			{
                d = p / q;
                u = x + d;
                if (u - a < tol2 || b - u < tol2)
				{
                    d = fabs(tol1) * sgn(xm - x);
                }
                done = 0;
            }
        }
        if (done)
		{
            if (x >= xm)
			{
                e = a - x;
			}
            else
			{
                e = b - x;
            }
            d = cgold * e;
        }
        if (fabs(d) >= tol1)
		{
            u = x + d;
		}
        else
		{
            u = x + fabs(tol1) * sgn(d);
        }
        fu = func(u);
        if (fu <= fx)
		{
            if (u >= x)
			{
                a = x;
			}
            else
			{
                b = x;
            }
            v = w;
            fv1 = fw;
            w = x;
            fw = fx;
            x = u;
            fx = fu;
		}
        else
		{
            if (u < x)
			{
                a = u;
			}
            else
			{
                b = u;
            }
            if (fu <= fw || w == x)
			{
                v = w;
                fv1 = fw;
                w = u;
                fw = fu;
			}
            else
			{
				if (fu <= fv1 || v == x || v == w)
				{
					v = u;
					fv1 = fu;
				}
            }
        }
    }
    if (iter > itmax)
	{
		cout<<" brent exceed maximum iterations."<<endl;
	}
    xmin = x;
    return fx;
}