DFP算法C++

整理文档很辛苦,赏杯茶钱您下走!

免费阅读已结束,点击下载阅读编辑剩下 ...

阅读已结束,您可以下载文档离线阅读编辑

资源描述

//DFP算法2011/5/25习题五第六题#includeiostream#includemath.h#includeiomanipusingnamespacestd;intconstn=2;//正定二次函数的自变量个数doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2]);//输入变量为函数自变量初值voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1]);//Q[n][n]是二次函数的正定矩阵,但Q的第n+1列存储一次项的系数voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n]);//自变量初值,正定二次函数的各项系数,返回梯度的初值intH(doubleg0[n],doublec);//判别准则:返回结束,返回继续迭代voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3]);//t[3]中返回的是a,b,c的系数值.开始计算minf(Xk+tPk)时的步长t的值,由于这是n元二次函数所以minf(t)是关于t0的二次函数,先求二次方程a,b,c系数,用一阶导为零。doublet_c(doublet[3]);//二次函数一阶导为零计算t的值,t0voidmain(){doublef_xs[n+n+1+(n-1)*n/2]={4,1,-40,-12,136,0};//正定二次函数的各项系数,第一个为:X1^2系数,第二个为:X2^2系数,第三个为:X1系数,第四个为:X2系数,第五个为:常数项,第六个为:X1X2系数;//n元的系数存放类推doublex[n]={8,9};//函数自变量初值doublef0;//函数值doubleg0[n];//梯度的值doubleQ[n][n+1];//求梯度处值设置的中间变量,包含两部分:正定二次函数对应的实对称矩阵,还有一次项系数doublec=0.01;//H准则限值doublet[3];//返回求minf()时t的二次函数的a,b,c的系数值doublet_bc;//步长doublep[n];//保存下降方向doubleH0[n][n];//保存模拟Hesse矩阵的逆doubley[n];//y(k)=g0(k+1)-g0(k)doubles[n];//s(k)=X(k+1)-X(k)doubles_temp[n][n]={0};//计算保存矩阵doubles_temp2[n][n]={0};doubles_temp3[n][n]={0};doubles_tl[n]={0};doubletemp;//临时值inti,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;elseH0[i][j]=0;}}for(i=0;in;i++)p[i]=(-1)*g0[i];k=0;//step2;do{abc(x,p,f_xs,t);//开始计算minf(Xk+tPk)时的步长t的值,t_bc=t_c(t);//求一阶导来计算tfor(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);//step3;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//shavedone!temp=0;//初值for(i=0;in;i++)temp+=s[i]*y[i];for(i=0;in;i++){for(j=0;jn;j++){s_temp[i][j]=s[i]*s[j]/temp;}}//求出S(k)S(k)t/S(k)tykfor(i=0;in;i++)//初值s_tl[i]=0;for(i=0;in;i++){for(j=0;jn;j++){s_tl[i]+=y[j]*H0[j][i];}}temp=0;//初值for(i=0;in;i++){temp+=s_tl[i]*y[i];}//这时s_tl[n]和s_temp2[n][n]可以利用for(i=0;in;i++)//初值s_tl[i]=0;for(i=0;in;i++){for(j=0;jn;j++){s_tl[i]+=H0[i][j]*y[j];}}for(i=0;in;i++)//初值for(j=0;jn;j++)s_temp2[i][j]=0;for(i=0;in;i++){for(j=0;jn;j++){s_temp2[i][j]=s_tl[i]*y[j];}}for(i=0;in;i++)//初值for(j=0;jn;j++)s_temp3[i][j]=0;for(i=0;in;i++){for(j=0;jn;j++){for(intii=0;iin;ii++){s_temp3[i][j]+=s_temp2[i][ii]*H0[ii][j];}}}for(i=0;in;i++){for(j=0;jn;j++){H0[i][j]=H0[i][j]+s_temp[i][j]-s_temp3[i][j]/temp;}}for(i=0;in;i++){temp=0;for(j=0;jn;j++){temp+=H0[i][j]*g0[j];}p[i]=(-1)*temp;}k++;}else{f0=fun(x,f_xs);//这里x[n],g0[n]的值已经修改break;}}else{flag=1;break;}tap++;}while(H(g0,c)==0);if(flag==1)//符合梯度准则跳出{flag=0;break;}}while(k==n);coutDFP法endl;cout函数f(x1,x2)=4(X1-5)^2+(X2-6)^2.的极小点为:f=f0endl;cout自变量取值为:endl;for(i=0;in;i++){coutxi+1=x[i]endl;}cout迭代次数为:tapendl;system(pause);}doublefun(doublex[n],doublef_xs[n+n+1+(n-1)*n/2])//输入变量为函数自变量初值{inti,j;doublef=0;for(i=0;in;i++)//计算X^2部分{f+=pow(x[i],2)*f_xs[i];}for(;i2*n;i++)//计算X部分{f+=x[i%n]*f_xs[i];}f+=f_xs[i];//计算常数项部分for(i=0;in;i++)//计算XiXj部分{for(j=i+1;jn;j++){f+=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1]*x[i]*x[j];}}returnf;}voidQ_fun(doublef_xs[n+n+1+(n-1)*n/2],doubleQ[n][n+1])//Q[n][n]是二次函数的正定矩阵,但Q的第n+1列存储一次项的系数{inti,j;for(i=0;in;i++){Q[i][i]=2*f_xs[i];}for(i=0;in;i++){Q[i][n]=f_xs[n+i];}for(i=0;in;i++){for(j=i+1;jn;j++){Q[j][i]=Q[i][j]=f_xs[(n+1)+n*(i+1)-i*(i+1)/2+j-i-1];}}}voidD_fun(doublex[n],doubleQ[n][n+1],doubleg0[n])//自变量初值,正定二次函数的各项系数,返回梯度的初值{inti,j;for(i=0;in;i++){g0[i]=0;//清零for(j=0;jn;j++){if(i==j)g0[i]+=Q[i][j]*x[j];elseg0[i]+=Q[i][j]*x[j];}}for(i=0;in;i++){g0[i]+=Q[i][n];}coutendl梯度g0的值=endl;for(i=0;in;i++)coutsetw(10)g0[i]endl;coutendl;}intH(doubleg0[n],doublec)//判别准则:返回结束,返回继续迭代{doubles=0;for(inti=0;in;i++){s+=pow(g0[i],2);}if(sqrt(s)c)return1;elsereturn0;}voidabc(doublex[n],doublep[n],doublef_xs[n+n+1+(n-1)*n/2],doublet[3])//t[3]中返回的是a,b,c的系数值{inti,j;t[0]=t[1]=t[2]=0;t[0]=fun(p,f_xs)-f_xs[2*n];for(i=n;i2*n;i++){t[0]-=f_xs[i]*p[i%n];}for(i=0;in;i++){t[1]+=2*f_xs[i]*x[i]*p[i];}for(;i2*n;i++){t[1]+=f_xs[i]*p[i%n];}for(i=0;in;i++){for(j=i+1;jn;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]);}}t[2]=fun(x,f_xs);coutendl二次函数minf(Xk)的二次项系数:t[0]一次项系数:t[1]常数项系数:t[2]endlendl;}doublet_c(doublet[3]){cout用一阶导求得步长为:t=-t[1]/t[0]/2endl;return-t[1]/t[0]/2;}运行结果

1 / 12
下载文档,编辑使用

©2015-2020 m.777doc.com 三七文档.

备案号:鲁ICP备2024069028号-1 客服联系 QQ:2149211541

×
保存成功