extremum.inl

特征值和特征向量的计算,每种算法都用c++以函数形式实现

        for(i=1; i<2*n; i++)
          if(xx(n,i)>fr)
          {
			  r = i;
			  fr = xx(n,i);
		  }
        g = 0;
		j = 0;
		fg = xx(n,0);
        if(r==0)
        {
			g = 1;
			j = 1;
			fg = xx(n,1);
		}
        for(i=j+1; i<2*n; i++)
          if(i!=r)
            if(xx(n,i)>fg)
            {
				g = i;
				fg = xx(n,i);
			}
        for(i=0; i<n; i++)
        {
			xf[i] = 0.0;
            for(j=0; j<2*n; j++)
              if(j!=r)	xf[i]=xf[i]+xx(i,j)/(2.0*n-1.0);
            xt[i]=(1.0+alpha)*xf[i]-alpha*xx(i,r);
        }
        jt = 1;
        while(jt==1)
        {
			jt = 0;
            z = FunctionValueTarget(xt,n);
            while(z>fg)
            {
				for(i=0; i<n; i++)	xt[i]=(xt[i]+xf[i])/2.0;
                z = FunctionValueTarget(xt,n);
            }
            j=0;
            for(i=0; i<n; i++)
            { 
				if(a[i]>xt[i])
                {
					xt[i]=xt[i]+0.000001;
					j=1;
				}
				if(b[i]<xt[i])
                {
					xt[i]=xt[i]-0.000001;
					j=1;
				}
            }
            if(j!=0) jt=1;
            else
            {
				ConditionValue(n,m,xt,c,d,w);
                j=0;
				kt=1;
                while((kt==1)&&(j<m))
                {
					if((c[j]<=w[j])&&(d[j]>=w[j])) j=j+1;
                    else kt=0;
                }
                if(j<m)
                { 
					for(i=0; i<n; i++)	xt[i]=(xt[i]+xf[i])/2.0;
                    jt=1;
                }
            }
        }
        for(i=0; i<n; i++)	xx(i,r)=xt[i];
        xx(n,r) = z;
        fr = 0.0; 
		fg = 0.0;
        for(j=0; j<2*n; j++)
        {
			fj = xx(n,j);
            fr = fr + fj / (2.0*n);
            fg = fg + fj * fj;
        }
        fr = (fg - 2.0 * n * fr * fr) / (2.0 * n - 1.0);
        if(fr>=eps)
        {
			kk = kk + 1;
            if(kk<k) it = 1;
        }
    }
    for(i=0; i<n; i++)
    {
		x[i] = 0.0;
        for(j=0; j<2*n; j++)	x[i]=x[i]+xx(i,j)/(2.0*n);
    }
    z = FunctionValueTarget(x,n);
	x[n] = z;
    return(kk);
}

//确定1维函数极小值点所在区间
template <class _Ty>
void MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, _Ty& fa, _Ty& fb, _Ty& fc)
{
    _Ty r,q,dum,gold = 1.0 + GoldNo;
    _Ty u,ulim,fu,tiny = 1e-20;
	int glimit = 100;

    fa = FunctionValue(ax);
    fb = FunctionValue(bx);
    if (fb > fa)
	{
        dum = ax;
        ax = bx;
        bx = dum;
        dum = fb;
        fb = fa;
        fa = dum;
    }
    cx = bx + gold * (bx - ax);
    fc = FunctionValue(cx);
	while (fb >= fc)
	{
        r = (bx - ax) * (fb - fc);
        q = (bx - cx) * (fb - fa);
        dum = q - r;
        if (Abs(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 = FunctionValue(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 = FunctionValue(u);
		}
        else
		{
			if ((cx - u) * (u - ulim) > 0)
			{
				fu = FunctionValue(u);
				if (fu < fc)
				{
					bx = cx;
					cx = u;
					u = cx + gold * (cx - bx);
					fb = fc;
					fc = fu;
					fu = FunctionValue(u);
				}
            }
			else
			{
				if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
				{
					u = ulim;
					fu = FunctionValue(u);
				}
				else
				{
					u = cx + gold * (cx - bx);
					fu = FunctionValue(u);
				}
			}
		}
		ax = bx;
		bx = cx;
		cx = u;
		fa = fb;
		fb = fc;
		fc = fu;
	}
}

//确定一维函数极小值点所在区间(重载)
template <class _Ty>
int MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, int& sign, 
									_Ty b, _Ty& fa, _Ty& fb, _Ty& fc)
{
    _Ty r,q,dum,gold = GoldNo + 1.0;
    _Ty u,ulim,fu,tiny = 1e-20;
	int glimit = 100;

    fa = FunctionValue(ax);
    fb = FunctionValue(bx);
    if (fb > fa)
	{
		sign = -1;
		fa *= sign;
		fb *= sign;
    }
    cx = bx + gold * (bx - ax);
	fc = sign * FunctionValue(cx);
	while (fb >= fc)
	{
        r = (bx - ax) * (fb - fc);
        q = (bx - cx) * (fb - fa);
        dum = q - r;
        if (Abs(dum) < tiny)
		{
			dum = tiny;
		}
        u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
        ulim = bx + glimit * (cx - bx);
		if(u > b)	//超出f(x)右边界
		{
			cout<<"  The MinimizerInterval exceed right border. Program is over!"<<endl;
			return 0;
		}
        if ((bx - u) * (u - cx) > 0)
		{
            fu = sign * FunctionValue(u);
            if (fu < fc)
			{
                ax = bx;
                fa = fb;
                bx = u;
                fb = fu;
                return 1;
			}
            else
			{
				if (fu > fb)
				{
					cx = u;
					fc = fu;
					return 1;
				}
            }
            u = cx + gold * (cx - bx);
            fu = sign * FunctionValue(u);
		}
        else
		{
			if ((cx - u) * (u - ulim) > 0)
			{
				fu = sign * FunctionValue(u);
				if (fu < fc)
				{
					bx = cx;
					cx = u;
					u = cx + gold * (cx - bx);
					fb = fc;
					fc = fu;
					fu = sign * FunctionValue(u);
				}
            }
			else
			{
				if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
				{
					u = ulim;
					fu = sign * FunctionValue(u);
				}
				else
				{
					u = cx + gold * (cx - bx);
					fu = sign * FunctionValue(u);
				}
			}
		}
		ax = bx;
		bx = cx;
		cx = u;
		fa = fb;
		fb = fc;
		fc = fu;
	}
	return 1;
}

//1维函数极小值的黄金分割法
template <class _Ty>
_Ty ExtremumGold1D(_Ty ax, _Ty bx, _Ty cx, int sign, _Ty tol, _Ty& xmin)
{
    _Ty x0,x1,x2,x3,f0,f1,f2,f3,r = GoldNo;
    _Ty c = 1.0 - GoldNo;
    x0 = ax;
    x3 = cx;
    if (Abs(cx - bx) > Abs(bx - ax))
	{
		x1 = bx;
        x2 = bx + c * (cx - bx);
	}
    else
	{
        x2 = bx;
        x1 = bx - c * (bx - ax);
    }
    f1 = sign * FunctionValue(x1);
    f2 = sign * FunctionValue(x2);
    while (Abs(x3 - x0) > tol * (Abs(x1) + Abs(x2)))
	{
        if (f2 < f1)
		{
            x0 = x1;
            x1 = x2;
            x2 = r * x1 + c * x3;
            f0 = f1;
            f1 = f2;
            f2 = sign * FunctionValue(x2);
		}
        else
		{
            x3 = x2;
            x2 = x1;
            x1 = r * x2 + c * x0;
            f3 = f2;
            f2 = f1;
            f1 = sign * FunctionValue(x1);
        }
    }
    if (f1 < f2)
	{
        xmin = x1;
		if (sign == -1)
			f1 = sign * f1;
        return f1;
	}
    else
	{
        xmin = x2;
		if (sign == -1)
			f2 = sign * f2;
        return f2;
    }
}

//1维函数极小值点的不用导数布伦特(Brent)法
template <class _Ty>
_Ty ExtremumBrentNonDerivative1D(_Ty ax, _Ty bx, _Ty cx, int sign, 
												_Ty tol, _Ty& xmin)
{
    int done,iter,itmax = 100;
    _Ty d,fu,r,q,p,xm,tol1,tol2,a,b,cgold = 1.0 - GoldNo;
    _Ty u,etemp,dum,v,w,x,e,fx,fv1,fw,zeps = DOUBLEERROR;
    a = ax;
    if (cx < ax)
	{
		a = cx;
	}
    b = ax;
    if (cx > ax)
	{
		b = cx;
	}
    v = bx;
    w = v;
    x = v;
    e = 0.0;
    fx = sign * FunctionValue(x);
    fv1 = fx;
    fw = fx;
    for (iter = 1; iter<=itmax; iter++)
	{
        xm = 0.5 * (a + b);
        tol1 = tol * Abs(x) + zeps;
        tol2 = 2.0 * tol1;
        if (Abs(x - xm) <= tol2 - 0.5 * (b - a))
		{
			break;
		}
        done = -1;
        if (Abs(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 = Abs(q);
            etemp = e;
            e = d;
            dum = Abs(0.5 * q * etemp);
            if (Abs(p) < dum && p > q * (a - x) && p < q * (b - x))
			{
                d = p / q;
                u = x + d;
                if (u - a < tol2 || b - u < tol2)
				{
                    d = Abs(tol1) * Sgn(xm - x);
                }
                done = 0;
            }
        }
        if (done)
		{
            if (x > xm || FloatEqual(x,xm))
			{
                e = a - x;
			}
            else
			{
                e = b - x;
            }
            d = cgold * e;
        }
        if (Abs(d) > tol1 || FloatEqual(d,tol1))
		{
            u = x + d;
		}
        else
		{
            u = x + Abs(tol1) * Sgn(d);
        }
        fu = sign * FunctionValue(u);
        if (fu < fx || FloatEqual(fu,fx))
		{
            if (u > x || FloatEqual(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 || FloatEqual(fu,fw) || FloatEqual(w,x))
			{
                v = w;
                fv1 = fw;
                w = u;
                fw = fu;
			}
            else
			{
				if (fu < fv1 || FloatEqual(fu,fv1) ||  FloatEqual(v,x) ||  FloatEqual(v,w))
				{
					v = u;
					fv1 = fu;
				}
            }
        }
    }
    if (iter > itmax)
	{
		cout<<" ExtremumBrentNonDerivative1D exceed maximum iterations."<<endl;
	}
    xmin = x;
    if(sign == -1) 
	{
		fx *= sign;
	}
	return fx;
}

#endif		// _EXTREMUM_INL