dfp.cpp

DFP变尺度法

//dfp.cpp
#include"dfp.h"

CMatrix CDfp::operator()(double(*p)(CMatrix&),CMatrix& X0,double e)
{
	pfun=p,n=X0.getrow();                
	X.set(n,1),S.set(n,1);
	X=X0, error=e;
	CMatrix B(n,n),B0(n,n);				 //B--变尺度矩阵
	CMatrix TiDu(n,1),TiDu0(n,1);	     //TiDu--函数梯度
	CMatrix DG(n,1),DGT(1,n),DX(n,1),DXT(1,n);
	CMatrix E(n,n);                      //E--校正矩阵 
	CDiffer differ;
	double a;                            //a--迭代步长
	double dfm;						     //dfm--梯度的模
	int flag(0);
	int k(1),c(1);

	for(int i=1;i<=n;i++)
		B(i,i)=1; 
	B0=B;
	TiDu=differ(pfun,X,error/10);
	dfm=TiDu.getmod();
	ofstream out(".//result.txt",ios::ate); //用于输出过程值
	for(k;k<=n;k++)
	{
		S=(-1)*B*TiDu;
		a=find_a(TiDu);
		X=X+a*S;
		TiDu0=TiDu,	TiDu=differ(pfun,X,error);
		dfm=TiDu.getmod ();
		out<<"k="<<c++<<"     "<<"[X]'=["<<X<<"]       "<<"[TiDu]'=["
		   <<TiDu<<"]     "<<"f(x)="<<pfun(X)<<"       a="<<a<<endl;
		if(dfm<error) 
			return X;	  
		else
		{
			if(k=n)	   
				k=1,   B=B0; 
			else
			{
				DG=TiDu-TiDu0,  DGT=DG.invert();
				DX=a*S,		DXT=DX.invert();
				E=(DX*DXT)/DX.mult(DGT)-(B*DG*DGT*B)/DGT.mult(B*DG);
				B=B+E,		k+=1;
			}
		}
	}
	return X;
}

double CDfp::find_a(CMatrix TiDu)
{
	double a=explore(1,TiDu);
	return a;
}

double CDfp::f(double a)
{
	return pfun(X+a*S);
}

double CDfp::explore0 (double a,double b)
{
	double x1,x2;
	x1=a+0.382*(b-a);
	x2=a+0.618*(b-a);
	while((x2-x1)>=error)
	{
		if(f(x1)<f(x2))	
			b=x2;
		else 
			a=x1;
		x1=a+0.382*(b-a);
		x2=a+0.618*(b-a);
	}
	return f(x1)<f(x2) ? x1:x2;
}

double CDfp::explore(double a0,CMatrix TiDu)
{
	double h;
	double t=-2*pfun(X)/(TiDu.invert().mult(S));
	if(t>0 && t<1/S.getmod())
		h=t;
	else
		h=1/S.getmod();
	double a1(a0),a2(a0+h);
	double f1,f2;
	f1=f(a1), f2=f(a2);
	while(f2<f1)
	{
		h*=2, a1=a2, a2+=h;
		f1=f2,     f2=f(a2);	
		if(f2>f1)
			return explore0(a1,a2);
	}
	h/=4;
	a2=a1,a1-=h;
	f2=f1,f1=f(a1);
	while(f2>f1)
	{
		h*=2, a2=a1, a1-=h;
		f2=f1,    f1=f(a1);
	}
	return explore0(a1,a2);
}