d11r7.cpp

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

#include "iostream.h"
#include "math.h"

int ncom;
double pcom[51], xicom[51];

int sgn(double x)
{
	if (x>0)
	{
		return 1;
	}
	else
	{
		if(x<0)
		{
			return -1;
		}
	}
	return 0;
}

double func2(double x[], int n)
{
    double f = 0.0;
    f = f + (x[1] - 1) * (x[1] - 1);
    f = f + (x[2] - 1) * (x[2] - 1);
    f = f + (x[3] - 1) * (x[3] - 1);
    return f;
}

double f1dim(double x)
{
    double xt[51];
    for (int j = 1; j<=ncom; j++)
	{
        xt[j] = pcom[j] + x * xicom[j];
    }
    return func2(xt, ncom);
}

double func(double x)
{
    return f1dim(x);
}

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;
	}
}

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;
}

void linmin(double p[], double xi[], int n, double& fret)
{
	int j;
    double tol = 0.0001;
    ncom = n;
    for (j = 1; j<=n; j++)
	{
        pcom[j] = p[j];
        xicom[j] = xi[j];
    }
    double fa,fx,fb,bx,ax = 0.0;
    double xmin,xx = 1.0;
    mnbrak(ax, xx, bx, fa, fx, fb);
    fret = brent(ax, xx, bx, tol, xmin);
    for (j = 1; j<=n; j++)
	{
        xi[j] = xmin * xi[j];
        p[j] = p[j] + xi[j];
    }
}

void main()
{
    //program d11r7
    //driver for routine linmin
    int ndim = 3;
    double x,sr2,fret,pio2 = 1.5707963;
    double p[4], xi[4];
    cout<<endl;
    cout<<"minimum of a 3-d quadratic centered";
    cout<<endl;
    cout<<"at (1.0,1.0,1.0). minimum is found";
    cout<<endl;
    cout<<"along a series of radials."<<endl;
    cout<<"      x            y            z        minimum"<<endl;
    int i,j;
	cout.setf(ios::fixed|ios::left);
	cout.precision(6);
	for (i = 0; i<=10; i++)
	{
        x = pio2 * i / 10.0;
        sr2 = sqrt(2.0);
        xi[1] = sr2 * cos(x);
        xi[2] = sr2 * sin(x);
        xi[3] = 1.0;
        p[1] = 0.0;
        p[2] = 0.0;
        p[3] = 0.0;
        linmin(p, xi, ndim, fret);
        for (j = 1; j<=3; j++)
		{
			cout<<"  ";
			cout.width(11);
            cout<<p[j];
        }
		cout<<"  ";
		cout.width(11);
        cout<<fret;
		cout<<endl;
    }
}