DFP算法C++..docVIP

//DFP算法2011/5/25 #includeiostream #includemath.h #includeiomanip using namespace std; int const n=2;//正定二次函数的自变量个数 double fun(double x[n],double f_xs[n+n+1+(n-1)*n/2]);//输入变量为函数自变量初值 void Q_fun(double f_xs[n+n+1+(n-1)*n/2],double Q[n][n+1]);//Q[n][n]是二次函数的正定矩阵，但Q的第n+1列存储一次项的系数 void D_fun(double x[n],double Q[n][n+1],double g0[n]);//自变量初值，正定二次函数的各项系数，返回梯度的初值 int H(double g0[n],double c);//判别准则：返回结束，返回继续迭代 void abc(double x[n],double p[n],double f_xs[n+n+1+(n-1)*n/2],double t[3]);//t[3]中返回的是a,b,c的系数值.开始计算minf(Xk+tPk)时的步长t的值，由于这是n元二次函数所以minf(t)是关于t0的二次函数，先求二次方程a,b,c系数,用一阶导为零。 double t_c(double t[3]);//二次函数一阶导为零计算t的值，t0 void main() { double f_xs[n+n+1+(n-1)*n/2]={4,1,-40,-12,136,0};//正定二次函数的各项系数，第一个为：X1^2系数，第二个为：X2^2系数，第三个为：X1系数，第四个为:X2系数，第五个为：常数项，第六个为：X1X2系数； //n元的系数存放类推 double x[n]={8,9};//函数自变量初值 double f0;//函数值 double g0[n];//梯度的值 double Q[n][n+1];//求梯度处值设置的中间变量，包含两部分：正定二次函数对应的实对称矩阵，还有一次项系数 double c=0.01;//H准则限值 double t[3];//返回求minf()时t的二次函数的a,b,c的系数值 double t_bc;//步长 double p[n];//保存下降方向 double H0[n][n];//保存模拟Hesse矩阵的逆 double y[n];//y(k)=g0(k+1)-g0(k) double s[n];//s(k)=X(k+1)-X(k) double s_temp[n][n]={0};//计算保存矩阵 double s_temp2[n][n]={0}; double s_temp3[n][n]={0}; double s_tl[n]={0}; double temp;//临时值 int i,j,k,flag=0,tap=0;//迭代次数 Q_fun(f_xs,Q);//计算正定二次函数对应的实对称矩阵 f0=fun(x,f_xs);//求函数初值 D_fun(x,Q,g0);//返回梯度的初值 do { for(i=0;in;i++)//给H0[n][n]的处值赋单位矩阵 { for(j=0;jn;j++) { if(i==j) H0[i][j]=1; else H0[i][j]=0; } } for(i=0;in;i++) p[i]=(-1)*g0[i]; k=0;//step 2; do { abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值， t_bc=t_c(t);//求一阶导来计算t for(i=0;in;i++) { x[i]=x[i]+t_bc*p[i]; s[i]=t_bc*p[i];//保存计算之值X(k+1)-X(k) } for(i=0;in;i++) y[i]=g0[i];//保存之类的 f0=fun(x,f_xs); D_fun(x,Q,g0);//step 3; for(i=0;in;i++) y[i]=g0[i]-y[i];//保存计算g0(k+1)-g0(k) if(H(g0,c)==0)//即不满足小于c { if(k!=n) { //y //s have done! temp=0;//初值 for(i=0;in;