共轭梯度法cpp2.txt

我编写的一个关于解决非线性方程的优化问题的共轭梯度法

#define N 2 
#define EPS 0.001  //判断收敛用的估算值    
#include <iostream.h>
#include <math.h>

double *prmchg(double *x,double ambda,double *p);
double funct(double *x);
double *grad(double *x);
double linear_search(double *start,double *direction);
double mole(double *m);
void output(double *x,double *grad);
int cnjgrd(double *x);

int k=0;

void main()
{
	double x[N]={0,3};
    int result;
	cout<<"the process of iteration:"<<endl;
	result=cnjgrd(x);
	cout<<"the final result:"<<endl<<"fmin="<<funct(x)<<",x[1]="<<x[0]<<",x[2]="<<x[1]<<endl;	
}

double *prmchg(double *x,double ambda,double *p)
{
	int j=0;
	for(j=0;j<N;j++)
		x[j]=x[j]+ambda*p[j];
	return x;
}

double funct(double *x)//求在给定点的函数值
{
	//return x[0]*x[0]+x[0]*x[1]+x[1]*x[1]; 
	return (x[0]-2)*(x[0]-2)*(x[0]-2)*(x[0]-2)+(x[0]-2*x[1])*(x[0]-2*x[1]);
}

double *grad(double *x)//求在给定点的梯度值
{
	double *g;
	g=new double[N];

	//g[0]=2*x[0]+x[1];
	//g[1]=x[0]+2*x[1];

    g[0]=4*(x[0]-2)*(x[0]-2)*(x[0]-2)+2*(x[0]-2*x[1]);
	g[1]=-4*(x[0]-2*x[1]);

	return g;
}

double mole(double *m)//求向量的模
{
	int i;
	double sum=0.0;
    for(i=1;i<=N;i++)
         sum=sum+m[i-1]*m[i-1];
	sum=sqrt(sum);
	return sum;
}
void output(double *x,double *grad)//输出函数
{
   int i;
   cout<<"k="<<k;

   for(i=1;i<=N;i++)
         cout<<",x["<<i<<"]="<<x[i-1];
	  
   cout<<",f(x)="<<funct(x)<<",mgrad="<<mole(grad)<<endl;
   k++;
}

double linear_search(double *start,double *direction)//三次样条插值求步长 
{
	int i;
	double step=0.001;
	double a=0,ya,va; //ya—a点的值；va—a点的倒数
	double b,yb,vb;
	double t,yt,vt;
	double s,z,w;
	double *g,*temp_point;

	g=new double[N];
	temp_point=new double[N];
    
	g=grad(start);
	va=0;

	for(i=1;i<=N;i++)//计算a点导数值；
		va=va+g[i-1]*direction[i-1];
	
	//do—while循环用于确定搜索区间；
	do
	{
		b=a+step;
		for(i=1;i<=N;i++)
			temp_point[i-1]=start[i-1]+b*direction[i-1];
		g=grad(temp_point);
		vb=0;
		for(i=1;i<=N;i++)
			vb=vb+g[i-1]*direction[i-1];
		if(fabs(vb)<1E-10)
		{
			delete[] g;
			delete[] temp_point;
			return b;
		}
		if(vb<-1E-15)
		{
			a=b;
			va=vb;
			step=2*step;
		}		
	}
	while(vb<=1E-15);

	//计算a点函数值；
   for(i=1;i<=N;i++)
		temp_point[i-1]=start[i-1]+a*direction[i-1];
   ya=funct(temp_point);
   
   //计算b点函数值；
   for(i=1;i<=N;i++)
		temp_point[i-1]=start[i-1]+b*direction[i-1];
   yb=funct(temp_point);

     //do—while循环用于确定步长t；
     do
   {
	   s=3*(yb-ya)/(b-a);
	   z=s-va-vb;
	   w=sqrt(z*z-va*vb);
	   t=a+(w-z-va)*(b-a)/(vb-va+2*w);
	   yb=funct(temp_point);
	   for(i=1;i<=N;i++)
		   temp_point[i-1]=start[i-1]+t*direction[i-1];
       yt=funct(temp_point);
       g=grad(temp_point);
	   vt=0;
       for(i=1;i<=N;i++)
		   vt=vt+g[i-1]*direction[i-1];
	   if(vt>1E-6)
	   {
		   b=t;
		   yb=yt;
		   vb=vt;
	   }
	   else if(vt<-1E-6)
	   {
		   a=t;
		   ya=yt;
		   va=vt;
	   }
	   else break;
   }
   while(fabs(vt)>=1E-6&&(fabs(b-a)>1E-6));
   
   delete[] g;
   delete[] temp_point;
   
   return t;
}

int cnjgrd(double *x)//共轭梯度法子程序
{
	int i;
	double ambda,beta;
	double *g,*p,*oldgrad;
    
	g=new double[N];
	oldgrad=new double[N];
    p=new double[N];
goon:		
	    g=grad(x);	    
		for(i=1;i<=N;i++)
	        oldgrad[i-1]=g[i-1];	    
		output(x,g);	    
		for(i=1;i<=N;i++) //计算初始方向向量；
		    p[i-1]=-1*g[i-1];	    
		ambda=linear_search(x,p);//计算搜索步长；		 
		x=prmchg(x,ambda,p);	//计算新点；
		
		int l=1;
		do
		{
		    g=grad(x);
	        output(x,g);
			beta=(mole(g)*mole(g))/(mole(oldgrad)*mole(oldgrad));//计算方向乘因子；			
			for(i=1;i<=N;i++)
	            oldgrad[i-1]=g[i-1];		    
			for(i=1;i<=N;i++) //计算方向向量；
				p[i-1]=-1*g[i-1]+beta*p[i-1];
			ambda=linear_search(x,p);
			x=prmchg(x,ambda,p);
		}
		while(l++!=N-1);

	if(mole(g)>EPS)
		goto goon;//迭代次数大于向量维数后，把当前点再作为初始点，重新开始迭代；
	else 
		output(x,g);    
    return 0;
}