mnbrak.cpp

vc++常用数值算法集合

void mnbrak(double& ax, double& bx, double& cx, double& fa,
			double& fb, double& fc)
{
    double r,q,dum,gold = 1.618034;
    int glimit = 100;
    double u,ulim,fu,tiny = 1e-20;
    fa = func(ax);
    fb = func(bx);
    if (fb > fa)
	{
        dum = ax;
        ax = bx;
        bx = dum;
        dum = fb;
        fb = fa;
        fa = dum;
    }
    cx = bx + gold * (bx - ax);
    fc = func(cx);
	while (fb >= fc)
	{
        r = (bx - ax) * (fb - fc);
        q = (bx - cx) * (fb - fa);
        dum = q - r;
        if (fabs(dum) < tiny)
		{
			dum = tiny;
		}
        u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
        ulim = bx + glimit * (cx - bx);
        if ((bx - u) * (u - cx) > 0)
		{
            fu = func(u);
            if (fu < fc)
			{
                ax = bx;
                fa = fb;
                bx = u;
                fb = fu;
                return;
			}
            else
			{
				if (fu > fb)
				{
					cx = u;
					fc = fu;
					return;
				}
            }
            u = cx + gold * (cx - bx);
            fu = func(u);
		}
        else
		{
			if ((cx - u) * (u - ulim) > 0)
			{
				fu = func(u);
				if (fu < fc)
				{
					bx = cx;
					cx = u;
					u = cx + gold * (cx - bx);
					fb = fc;
					fc = fu;
					fu = func(u);
				}
            }
			else
			{
				if ((u - ulim) * (ulim - cx) >= 0)
				{
					u = ulim;
					fu = func(u);
				}
				else
				{
					u = cx + gold * (cx - bx);
					fu = func(u);
				}
			}
		}
		ax = bx;
		bx = cx;
		cx = u;
		fa = fb;
		fb = fc;
		fc = fu;
	}
}