拟牛顿法.docx

1、拟牛顿法拟牛顿法（变尺度法）DFP算法的c/c+源码#include iostream.h#include math.hvoid comput_grad(double (*pf)(double *x), int n, double *point, double *grad); /计算梯度double line_search1(double (*pf)(double *x), int n, double *start, double *direction); /0.618法线搜索double line_search(double (*pf)(double *x), int n, double *

2、start, double *direction); /解析法线搜索double DFP(double (*pf)(double *x), int n, double *min_point); /无约束变尺度法/梯度计算模块/参数：指向目标函数的指针，变量个数，求梯度的点，结果void comput_grad(double (*pf)(double *x), int n, double *point, double *grad) double h=1E-3; int i; double *temp; temp = new doublen; for(i=1;i=n;i+) tempi-1=poi

3、nti-1; for(i=1;i=n;i+) tempi-1+=0.5*h; gradi-1=4*pf(temp)/(3*h); tempi-1-=h; gradi-1-=4*pf(temp)/(3*h); tempi-1+=(3*h/2); gradi-1-=(pf(temp)/(6*h); tempi-1-=(2*h); gradi-1+=(pf(temp)/(6*h); tempi-1=pointi-1; delete temp;/一维搜索模块/参数：指向目标函数的指针，变量个数，出发点，搜索方向/返回：最优步长double line_search( double (*pf)(doubl

4、e *x), int n, double *start, double *direction) int i; double step=0.001; double a=0,value_a,diver_a; double b,value_b,diver_b; double t,value_t,diver_t; double s,z,w; double *grad,*temp_point; grad=new doublen; temp_point=new doublen; comput_grad(pf,n,start,grad); diver_a=0; for(i=1;i=n;i+) diver_a

5、=diver_a+gradi-1*directioni-1; do b=a+step; for(i=1;i=n;i+) temp_pointi-1=starti-1+b*directioni-1; comput_grad(pf,n,temp_point,grad); diver_b=0; for(i=1;i=n;i+) diver_b=diver_b+gradi-1*directioni-1; if( fabs(diver_b)1E-10 ) delete grad; delete temp_point; return(b); if( diver_b-1E-15 ) a=b; diver_a=

6、diver_b; step=2*step; while( diver_b=1E-15 ); for(i=1;i=n;i+) temp_pointi-1=starti-1+a*directioni-1; value_a=(*pf)(temp_point); for(i=1;i=n;i+) temp_pointi-1=starti-1+b*directioni-1; value_b=(*pf)(temp_point); do s=3*(value_b-value_a)/(b-a); z=s-diver_a-diver_b; w=sqrt( fabs(z*z-diver_a*diver_b) );

7、/! t=a+(w-z-diver_a)*(b-a)/(diver_b-diver_a+2*w); value_b=(*pf)(temp_point); for(i=1;i=n;i+) temp_pointi-1=starti-1+t*directioni-1; value_t=(*pf)(temp_point); comput_grad(pf,n,temp_point,grad); diver_t=0; for(i=1;i1E-6) b=t; value_b=value_t; diver_b=diver_t; else if(diver_t=1E-6) & (fabs(b-a)1E-6) )

8、; delete grad; delete temp_point; return(t); /无约束变尺度法DFP函数声明/参数：pf指向目标函数的指针，n变量个数，min_point接受初始点、存放结果/返回：极小点处函数值/double DFP( double (*pf)(double *x), int n, double *min_point ) int i,j; int k=0; double e=1E-5; double g_norm; double *g0=new doublen; /梯度 double *g1=new doublen; double *dg=new doublen;

9、 double *p=new doublen; /搜索方向 =-g double t; /一维搜索步长 double *x0=new doublen; double *x1=new doublen; double *dx=new doublen; double *H=new double*n; for (i=0; in; i+) Hi = new doublen; double *tempH=new double*n; for (i=0; in; i+) tempHi = new doublen; double *gH=new doublen; double *Hg=new doublen;

10、double num1; double num2; for(i=0;in;i+) for(j=0;jn;j+) if(i=j) Hij=1.0; / H0=I else Hij=0.0; tempHij=0.0; for(i=0;in;i+) x0i=min_pointi; comput_grad(pf,n,x0,g0); g_norm=0.0; for(i=0;in;i+) g_norm=g_norm+g0i*g0i; g_norm=sqrt(g_norm); if (g_norme) for(i=0;in;i+) min_pointi=x0i; delete g0; delete g1;

11、delete dg; delete p; delete x0; delete x1; delete dx; for (i=0; in; i+) delete Hi; delete H; for (i=0; in; i+) delete tempHi; delete tempH; delete gH; delete Hg; return pf(min_point); for(i=0;in;i+) pi=-g0i; do t=line_search(pf,n,x0,p); for(i=0;in;i+) x1i=x0i+t*pi; comput_grad(pf,n,x1,g1); g_norm=0.

12、0; for(i=0;in;i+) g_norm=g_norm+g1i*g1i; g_norm=sqrt(g_norm); /coutk x00 x01 g_normn; if (g_norme) for(i=0;in;i+) min_pointi=x1i; delete g0; delete g1; delete dg; delete p; delete x0; delete x1; delete dx; for (i=0; in; i+) delete Hi; delete H; for (i=0; in; i+) delete tempHi; delete tempH; delete g

13、H; delete Hg; return pf(min_point); for(i=0;in;i+) dxi=x1i-x0i; dgi=g1i-g0i; /求Hk+1的矩阵运算 /g*H,H*g for(i=0;in;i+) gHi=0.0; Hgi=0.0; for(i=0;in;i+) for(j=0;jn;j+) gHi=gHi+dgj*Hji; /Hgi=Hgi+Hij*dgj; Hgi=gHi; /num1,num2 num1=0.0; num2=0.0; for(i=0;in;i+) num1=num1+dxi*dgi; num2=num2+gHi*dgi; /tempHij fo

14、r(i=0;in;i+) for(j=0;jn;j+) tempHij=0.0; for(i=0;in;i+) for(j=0;jn;j+) tempHij=tempHij+Hij; tempHij=tempHij+dxi*dxj/num1; tempHij=tempHij-Hgi*gHj/num2; for(i=0;in;i+) for(j=0;jn;j+) Hij=tempHij; / /P for(i=0;in;i+) pi=0.0; for(i=0;in;i+) for(j=0;jn;j+) pi=pi-Hij*g1j; for(i=0;ie); for(i=0;in;i+) min_

15、pointi=x1i; delete g0; delete g1; delete dg; delete p; delete x0; delete x1; delete dx; for (i=0; in; i+) delete Hi; delete H; for (i=0; in; i+) delete tempHi; delete tempH; delete gH; delete Hg; return pf(min_point);/double fun(double *x) return 100*( x1-x0*x0 )*( x1-x0*x0 ) + (1-x0)*(1-x0);void ma

16、in() int n=2; double min_point2=-5,10; double min_value=DFP(fun,n,min_point); coutmin_point0 and min_point1 min_value; /0.618法线搜索/参数：指向目标函数的指针，变量个数，出发点，搜索方向/返回：最优步长/double line_search1( double (*pf)(double *x), int n, double *start, double *direction) int i; int k; double l,a,b,c,u,lamda,t,e; double

17、 *xa=new doublen; double *xb=new doublen; double *xc=new doublen; double *xl=new doublen; double *xu=new doublen; /确定初始搜索区间 l=0.001; a=0; k=0; do k+; c=a+l; for(i=0;i= pf(xa) ); / ? k=0; do k+; l=2*l; b=c+l; for(i=0;in;i+) xci=starti+c*directioni; xbi=starti+b*directioni; a=c; c=b; while( pf(xb) = pf(xc) ); a=0; b=0.1; /寻优 t=0.618; e=0.000001; lamda=b-t*(b-a); u=a+t*(b-a); for(i=0;in;i+) xli=starti+lamda*directioni; xui=starti+u*directioni; k=0; do k+; if( pf(xl)=e ); lamda=(a+b)/2; delete xa; delete xb; delete xc; delete xl; delete xu; return lamda ;