newton法求n元二次正定函数最小值通用程序（四章2题）.txt

Newton法的本质就是不断用切线来近似曲线

#include<iostream.h>
#include<math.h>
#include<iomanip.h>
#define n 2//正定二次函数的自变量个数
double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值
{
	int i,j;
	double f=0;
	for(i=0;i<n;i++)//计算X^2部分
	{
		f+=pow(x[i],2)*f_xs[i];
	}
	for(;i<2*n;i++)//计算X部分
	{
		f+=x[i%n]*f_xs[i];
	}
	f+=f_xs[i];//计算常数项部分
	for(i=0;i<n;i++)//计算XiXj部分
	{
		for(j=i+1;j<n;j++)
		{
			f+=f_xs[(n+1)+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-1)*n/2],double Q[n][n+1])
{
	int i,j;
	for(i=0;i<n;i++)
	{
		Q[i][i]=2*f_xs[i];
	}
	for(i=0;i<n;i++)
	{
		Q[i][n]=f_xs[n+i];
	}
	for(i=0;i<n;i++)
	{
		for(j=i+1;j<n;j++)
		{
			Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];
		}
	}
}
void D_fun(double x[n],double Q[n][n+1],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];
		}
	}
	for(i=0;i<n;i++)
	{
		g0[i]+=Q[i][n];
	}
	cout<<"梯度g0的值="<<endl;
	for(i=0;i<n;i++)
		cout<<setw(10)<<g0[i]<<endl;
}

void Hesse_fun(double Q[n][n+1],double Hesse[n][n])
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			Hesse[i][j]=Q[i][j];
		}
	}
	cout<<"Hesse 矩阵:"<<endl;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			cout<<setw(10)<<Hesse[i][j]<<" ";
		}
		cout<<endl;
	}
}

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 inv(double a[n][n],double c[n][n])//求Hesse矩阵的逆矩阵
{
	//double a[n][n]={{2,0},{0,8}};//线性方程组的系数矩阵  
	//double c[n][n];//求逆矩阵
	double m[n];//辅助乘数
	double T;//存储换行临时变量
	double temp=0;
	double t;//最大列主元
	int tap1=0,tap2=0;//最大列主元下标

	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i==j)
			{
				c[i][j]=1;
			}
			else
			{
				c[i][j]=0;
			}
		}
	}
	//for(i=0;i<n;i++){for(int j=0;j<n;j++){ cout<<setw(3)<<c[i][j];}cout<<endl;}
	/*cout<<"原线性方程组系数矩阵为:"<<endl;
	for(i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{ 
			cout<<setw(3)<<a[i][j]<<" ";			
		}
		cout<<endl;
	}
	cout<<endl;*/

	for(int s=0;s<n;s++)
	{
		t=a[s][s];//赋初值
		for(int p=s;p<n;p++)//选取最大列主元
		{
			if(fabs(a[p][s])>t)
			{
				t=a[p][s];
				tap1=p;
				tap2=s;
			}
		}
		//cout<<"t="<<t<<" tap1="<<tap1<<" tap2="<<tap2<<endl;
		if(t==0)
		{cout<<"列主元为零，失败！"<<endl;break;}
		if(tap1!=tap2)//换行
		{
			for(int i=0;i<n;i++)////////////////////////////////s
			{
				T=a[tap1][i];
				a[tap1][i]=a[s][i];
				a[s][i]=T;
				T=c[tap1][i];//逆矩阵换行
				c[tap1][i]=c[s][i];
				c[s][i]=T;
			}
		}
		tap1=tap2=0;//置零
		/*cout<<"换行后："<<endl;for(i=0;i<n;i++){for(int j=0;j<n;j++){cout<<setw(3)<<a[i][j]<<" ";	}cout<<endl;}
		cout<<endl;*/
		for(int j=0;j<n;j++)//单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响)
		{
			a[s][j]=a[s][j]/t;
			c[s][j]=c[s][j]/t;
		}		
		for(int i=0;i<n;i++)//消元
		{
			if(i!=s)
			{
				m[i]=a[i][s]/a[s][s];
			    for(int j=0;j<n;j++)//注意：逆矩阵每列都消元，而原矩阵则不受影响
				{
					a[i][j]=a[i][j]-m[i]*a[s][j];
					c[i][j]=c[i][j]-m[i]*c[s][j];
				}
			}
		}
		/*cout<<"消元后："<<endl;for(i=0;i<n;i++){for(int j=0;j<n;j++){ cout<<setw(3)<<a[i][j]<<" ";		}cout<<endl;}cout<<endl;*/
	}
	/*cout<<"消元后的系数矩阵为:"<<endl;
	for(i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			cout<<setw(3)<<a[i][j]<<" ";
		}
	cout<<endl;
	}*/

	cout<<endl<<endl;
	cout<<"Hesse矩阵的逆矩阵为:"<<endl;
	for(i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{ 
			cout<<setw(10)<<c[i][j]<<" ";
		}
		cout<<endl;
	}

}

void xiang_cheng(double a[n][n],double b[n],double c[n])//n*n矩阵乘以n*1矩阵
{
	//double a[n][n]={{0.5,0},{0,0.125}};
	//double b[n]={2,8};
	//double c[n];//结果
	int i,j,k;
	/*cout<<"a矩阵为:"<<endl;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			cout<<setw(10)<<a[i][j];
		}
		cout<<endl;
	}
	cout<<"b矩阵为:"<<endl;
	for(i=0;i<n;i++)
	{
		cout<<setw(10)<<b[i]<<endl;
	}*/
	if(n==n)
	{
		cout<<"两矩阵阶数相符可以相乘！"<<endl;
		for(i=0;i<n;i++)
		{
			c[i]=0;
		}
		for(k=0;k<n;k++)//c矩阵的第k行
		{
			for(j=0;j<n;j++)
			{
				c[k]+=a[k][j]*b[j];
			}
		}
		cout<<"Hesse矩阵的逆矩阵和梯度向量相乘，两矩阵相乘结果为:"<<endl;
		for(i=0;i<n;i++)
		{
			cout<<setw(10)<<c[i]<<endl;
		}
	}
	else
		cout<<"两个矩阵阶数不符，不能相乘"<<endl;
}


void main()
{
	double f_xs[n+n+1+(n-1)*n/2]={1,1,-10,-4,60,-1};//正定二次函数的各项系数，第一个为：X1^2系数，第二个为：X2^2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数；
	//n元的系数存放类推
	double x[n]={0,0};//函数自变量初值
	double f0;//函数值
	double g0[n];//梯度的值
	double Q[n][n+1];//正定二次函数对应的实对称矩阵
	double c=0.01;//H准则限值
	double Hesse[n][n];//Hesse矩阵
	double inv_Hesse[n][n];//Hesse矩阵的逆
	double temp[n];//保存矩阵相乘结果

	int i,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]的值
	

		Hesse_fun(Q,Hesse);//返回Hesse矩阵的值？？一个二次函数的Hesse矩阵是变的还是不变
	    inv(Hesse,inv_Hesse);//返回Hesse矩阵的逆矩阵

	    xiang_cheng(inv_Hesse,g0,temp);

	    for(i=0;i<n;i++)
		{
			x[i]=x[i]-temp[i];
		}
			tap++;

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

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