拟牛顿法(变尺度法)DFP算法的cc源码

1、拟牛顿法(变尺度法)DFP算法的C/C+源码#include "iostream.h"#inClude "math.h"void Comput_grad(double (*pf)(double *x), int n,double *point,double *grad);/ 计算梯度double line_searCh1(double(*pf)(double*x), intn, double*start,double*direCtion); /0.618法线搜索double line_searCh(double(*pf)(double*x), intn,

2、 double*start,double*direCtion);/解析法线搜索double DFP(double (*pf)(double*x),intn, 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)(double *

4、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=diver_a+gr

5、adi-1*directioni-1;dob=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=diver_b

6、;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) );/!t

7、=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;i<=n;i+)diver_t=diver_t+gradi-1*directioni-1;if(diver_t>1E-6)b=t;value_b=value_t;diver

8、_b=diver_t;else if(diver_t<-1E-6)a=t;value_a=value_t;diver_a=diver_t;else break;while( (fabs(diver_t)>=1E-6) && (fabs(b-a)>1E-6) );delete grad;delete temp_point;return(t);/无约束变尺度法 DFP函数声明/n 变量个数， min_point 接受初始点、存放结果/ 参数： pf 指向目标函数的指针，/ 返回：极小点处函数值/double DFP(double (*pf)(double *x),

9、 int n,double *min_point)int i,j;int k=0;double e=1E-5;double *g0=new doublen; / double *g1=new doublen;double *dg=new doublen;double *p=new doublen; / double t; /梯度搜索方向 =-g一维搜索步长double *x0=new doublen; double *x1=new doublen; double *dx=new doublen;double *H=new double*n;for (i=0; i<n; i+)Hi = n

10、ew doublen;double *tempH=new double*n;for (i=0; i<n; i+)tempHi = new doublen;double *gH=new doublen; double *Hg=new doublen; double num1;double num2;for(i=0;i<n;i+)for(j=0;j<n;j+) if(i=j) Hij=1.0; / H0=I else Hij=0.0; tempHij=0.0;for(i=0;i<n;i+) x0i=min_pointi;comput_grad(pf,n,x0,g0);for

11、(i=0;i<n;i+) g_norm=g_norm+g0i*g0i;g_norm=sqrt(g_norm);if (g_norm<e)for(i=0;i<n;i+) min_pointi=x0i;delete Hi;delete tempHi;delete g0; delete g1; delete dg; delete p; delete x0; delete x1; delete dx; for (i=0; i<n; i+) delete H;for (i=0; i<n; i+) delete tempH; delete gH; delete Hg;retu

12、rn pf(min_point);for(i=0;i<n;i+) pi=-g0i;dot=line_search(pf,n,x0,p); for(i=0;i<n;i+) x1i=x0i+t*pi; comput_grad(pf,n,x1,g1);for(i=0;i<n;i+) g_norm=g_norm+g1i*g1i; g_norm=sqrt(g_norm);/cout<<k<<" "<<x00<<"<<"n"if (g_norm<e)for(i=0;i<

13、;n;i+) min_pointi=x1i;"<<x01<<"delete g0; delete g1; delete dg;delete p;delete x0;delete x1;delete dx; for (i=0; i<n; i+) delete H;for (i=0; i<n; i+) delete tempH; delete gH; delete Hg;return pf(min_point); for(i=0;i<n;i+) dxi=x1i-x0i; dgi=g1i-g0i;/delete Hi;delete tem

14、pHi;求 Hk+1 的矩阵运算/g*H,H*gfor(i=0;i<n;i+)gHi=0.0;Hgi=0.0;for(i=0;i<n;i+)for(j=0;j<n;j+)gHi=gHi+dgj*Hji;/Hgi=Hgi+Hij*dgj;Hgi=gHi;/num1,num2num1=0.0;num2=0.0;for(i=0;i<n;i+) num1=num1+dxi*dgi; num2=num2+gHi*dgi;/tempHijfor(i=0;i<n;i+)for(j=0;j<n;j+)tempHij=0.0;for(i=0;i<n;i+)for(j=0

15、;j<n;j+)tempHij=tempHij+Hij; tempHij=tempHij+dxi*dxj/num1; tempHij=tempHij-Hgi*gHj/num2;for(i=0;i<n;i+) for(j=0;j<n;j+) Hij=tempHij;/Pfor(i=0;i<n;i+) pi=0.0;for(i=0;i<n;i+)for(j=0;j<n;j+)pi=pi-Hij*g1j;for(i=0;i<n;i+)g0i=g1i;x0i=x1i; k=k+1;while(g_norm>e);for(i=0;i<n;i+) mi

16、n_pointi=x1i;delete g0;delete g1;delete dg;delete p;delete x0;delete x1;delete dx;for (i=0; i<n; i+)delete Hi;delete H;for (i=0; i<n; 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 main()int

17、 n=2;double min_point2=-5,10;double min_value=DFP(fun,n,min_point); cout<<min_point0<<" and "<<min_point1<<" "<<min_value;/0.618 法线搜索/ 参数：指向目标函数的指针，变量个数，出发点，搜索方向/ 返回：最优步长/double line_search1(double (*pf)(double *x),int n,double *start,double *directi

18、on)int i;int k;double l,a,b,c,u,lamda,t,e;double *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;dok+;c=a+l;for(i=0;i<n;i+)xci=starti+c*directioni; xai=starti+a*directioni;l=l/3;while( pf(xc) >= pf(xa) );/ ?k=0;dok+; l=2*l;b=c+l;for(i=0;i<n;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.00000