最速下降法是以负梯度方向作为下降方向的极小化算法

#include<iostream.h>
#include<math.h>
#define n 2
double fun(double x[n],double f_xs[n+(n-1)*n/2])//输入变量为函数自变量初值
{
	int i,j;
	double f=0;
	for(i=0;i<n;i++)
	{
		f+=x[i]*f_xs[i];
	}
	for(i=0;i<n;i++)
	{
		for(j=i+1;j<n;j++)
		{
			f+=f_xs[n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];
		}
	}
	return f;
}

void Q_fun(double f_xs[n+(n-1)*n/2],double Q[n][n])
{
	int i,j;
	for(i=0;i<n;i++)
	{
		Q[i][i]=2*f_xs[i];
	}
	for(i=0;i<n;i++)
	{
		for(j=i+1;j<n;j++)
		{
			Q[j][i]=Q[i][j]=f_xs[n*(i+1)-i*(i+1)/2+j-i-1];
		}
	}
}
void D_fun(double x[n],double Q[n][n],double g0[n])//自变量初值，正定二次函数的各项系数，返回梯度的初值
{
	int i,j;
	for(i=0;i<n;i++)
	{
		g0[i]=0;//清零
		for(j=0;j<n;j++)
		{
			if(i==j)
				g0[i]+=Q[i][j]*x[j];
			else 
				g0[i]+=Q[i][j]*x[j];
		}
	}
}

void x_k_1(double x[n],double g0[n],double Q[n][n])//二次函数显式迭代公式
{
	double gg=0;
	double gqg[n];
	double k=0;
	for(int i=0;i<n;i++)
	{
		gqg[i]=0;
	}

	for(i=0;i<n;i++)
	{
		gg+=g0[i]*g0[i];
	}//求出gkt_gk
	for(i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			gqg[i]+=g0[j]*Q[j][i];
		}
	}//求出gkt_Q
	for(i=0;i<n;i++)
	{
		k+=gqg[i]*g0[i];
	}//求出gkt_Q_gk
	k=gg/k;//求出(gkt_gk)/(gkt_Q_gk)
	for(i=0;i<n;i++)
	{
		x[i]=x[i]-k*g0[i];
	}
}

int H(double g0[n],double c)
{
	double s=0;
	for(int i=0;i<n;i++)
	{
		s+=pow(g0[i],2);
	}
	if(sqrt(s)<c)
		return 1;
	else 
		return 0;
}


void main()
{
	double f_xs[n+(n-1)*n/2]={1,4,0};//正定二次函数的各项系数，第一个为：X1^2系数，第二个为：X2^2系数，第三个为：X1X2系数；
	double x[n]={1,1};//函数自变量初值
	double f0;//函数值
	double g0[n];//梯度的值
	double Q[n][n];//正定二次函数对应的实对称矩阵
	double c=0.001;//H准则限值
	int tap=0;//迭代次数

	Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵ok!
	f0=fun(x,f_xs);//求函数初值ok!
	D_fun(x,Q,g0);//返回梯度的初值
	
	do
	{
		x_k_1(x,g0,Q);//求出下一个迭代点，更新了x[n]的值
		tap++;

		f0=fun(x,f_xs);
		D_fun(x,Q,g0);
	}while(H(g0,c)==0);

	cout<<"函数f(x1,x2)=x1^2+4x2^2的极小点为："<<"f="<<f0<<endl;
	cout<<"自变量取值为："<<endl;
	for(int i=0;i<n;i++)
	{
		cout<<"x"<<i+1<<"="<<x[i]<<endl;
	}
	cout<<"迭代次数为："<<tap<<endl;
}