dbrent.cpp

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

double dbrent(double ax,double bx,double cx,double tol,double& xmin)
{

    int done,iter,itmax = 100;
    double d,fu,xm,tol1,tol2,a,b,cgold = 0.381966;
    double u,v,w,x,e,fx,fv1,fw,zeps = 0.0000000001;
	double dx,dv,dw,d1,d2,u1,u2,olde,du;
	bool ok1,ok2;
    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;
    dx = deriv(x);
    dv = dx;
    dw = dx;
    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)))
		{
			done = -1;
			break;
		}
        else
		{
          done = 0;
		}
        if (fabs(e) > tol1)
		{
            d1 = 2.0 * (b - a);
            d2 = d1;
            if (dw != dx)
			{
				d1 = (w - x) * dx / (dx - dw);
			}
            if (dv != dx)
			{
				d2 = (v - x) * dx / (dx - dv);
			}
            u1 = x + d1;
            u2 = x + d2;
            ok1 = ((a - u1) * (u1 - b) > 0.0) && (dx * d1 <= 0.0);
            ok2 = ((a - u2) * (u2 - b) > 0.0) && (dx * d2 <= 0.0);
            olde = e;
            e = d;
            if (ok1 || ok2)
			{
                if (ok1 && ok2)
				{
                    d = d1;
				}
                else
				{
                    d = d2;
                }
			}
            else
			{
				if (ok1)
				{
					d = d1;
				}
				else
				{
					d = d2;
				}
            }
            if (fabs(d) > fabs(0.5 * olde))
			{
				u = x + d;
				if (u - a < tol2 || b - u < tol2)
				{
					d = fabs(tol1) * sgn(xm - x);
				}
            }
		}
        if (dx >= 0.0)
		{
            e = a - x;
		}
        else
		{
            e = b - x;
        }
        d = 0.5 * e;
        if (fabs(d) >= tol1)
		{
            u = x + d;
            fu = func(u);
		}
        else
		{
            u = x + fabs(tol1) * sgn(d);
            fu = func(u);
            if (fu > fx)
			{
				done = -1;
				break;
			}
            else
			{
				done = 0;
            }
        }
        du = deriv(u);
        if (fu <= fx)
		{
            if (u >= x)
			{
                a = x;
			}
            else
			{
                b = x;
            }
            v = w;
            fv1 = fw;
            dv = dw;
            w = x;
            fw = fx;
            dw = dx;
            x = u;
            fx = fu;
            dx = du;
		}
        else
		{
            if (u < x)
			{
                a = u;
			}
            else
			{
                b = u;
			}
            if (fu <= fw || w == x)
			{
                v = w;
                fv1 = fw;
                dv = dw;
                w = u;
                fw = fu;
                dw = du;
			}
            else
			{
				if (fu <= fv1 || v == x || v == w)
				{
					v = u;
					fv1 = fu;
					dv = du;
				}
			}
		}
	}	
    if (!done)
	{
        cout<<"dbrent exceeded maximum iterations."<<endl;
		exit(1);
	}
    else
	{
        xmin = x;
        return fx;
    }
}