修正牛顿法最优化C++程序

}
for(;i<2*n;i++)
{
t[1]+=f_xs[i]*p[i%n];
}
for(i=0;i<n;i++)
{
for(j=i+1;j<n;j++)
{
t[1]+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*(x[i]*p[j]+x[j]*p[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矩阵
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(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]<<" ";
}
}
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];
p[i]*=(-1);
abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值，
t_bc=t_c(t);//求一阶导来计算t
for(i=0;i<n;i++)
{
x[i]=x[i]+t_bc*p[i];
}
tap++;
}
}
}
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++)
{
int i,j,k;
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++)
{
double inv_Hesse[n][n];//Hesse矩阵的逆
double t[3];//返回求minf()时t的二次函数的a,b,c的系数值
double t_bc;//步长
double p[n];//保存矩阵相乘结果_加负号后变成下降方向
int i,tap=0;//迭代次数
{
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])
}
}
t[2]=fun(x,f_xs);
cout<<endl<<"二次函数minf(Xk)的二次项系数："<<t[0]<<" 一次项系数："<<t[1]<<" 常数项系数："<<t[2]<<endl<<endl;
}
//二次函数一阶导为零计算t的值，t>0
double t_c(double t[3])
}
cout<<endl;
}
}
int H(double g0[n],double c)//判别准则：返回1结束，返回0继续迭代
{
double s=0;
for(int i=0;i<n;i++)
{
s+=pow(g0[i],2);
}
if(sqrt(s)<c)
{
int i,j;
t[0]=t[1]=t[2]=0;
t[0]=fun(p,f_xs)-f_xs[2*n];
for(i=n;i<2*n;i++)
{
t[0]-=f_xs[i]*p[i%n];
}
for(i=0;i<n;i++)
{
t[1]+=2*f_xs[i]*x[i]*p[i];
#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])//输入变量为函数自变量初值
{
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<<endl<<"梯度g0的值="<<endl;
a[s][i]=T;
T=c[tap1][i];//逆矩阵换行
c[tap1][i]=c[s][i];
c[s][i]=T;
}
}
tap1=tap2=0;//置零
for(int j=0;j<n;j++)//单位化,(注意：逆矩阵每列都单位化，而原矩阵则不受影响)
Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵ok!
f0=fun(x,f_xs);//求函数初值ok!
D_fun(x,Q,g0);//返回梯度的初值
while(H(g0,c)==0)
{
//x_k_1(x,g0,Q);//求出下一个迭代点，更新了x[n]的值
for(i=0;i<n;i++)
cout<<setw(10)<<g0[i]<<endl;
cout<<endl;
}
void Hesse_fun(double Q[n][n+1],double Hesse[n][n])
{
int i,j;
for(i=0;i<n;i++)
{
for(i=0;i<n;i++)
{
cout<<"x"<<i+1<<"="<<x[i]<<endl;
}
cout<<"迭代次数为："<<tap<<endl;
//n元的系数存放类推
double x[n]={0,0};//函数自变量初值
double f0;//函数值