拟牛顿法Word格式文档下载.docx
《拟牛顿法Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《拟牛顿法Word格式文档下载.docx（17页珍藏版）》请在冰点文库上搜索。
temp[i-1]-=h;
grad[i-1]-=4*pf（temp）/（3*h）;
temp[i-1]+=（3*h/2）;
grad[i-1]-=（pf（temp）/（6*h））;
temp[i-1]-=（2*h）;
grad[i-1]+=（pf（temp）/（6*h））;
delete[]temp;
}
//一维搜索模块
指向目标函数的指针，变量个数，出发点，搜索方向
//返回：
最优步长
doubleline_search（
double（*pf）（double*x）,
double*start,
double*direction）
doublestep=0.001;
doublea=0,value_a,diver_a;
doubleb,value_b,diver_b;
doublet,value_t,diver_t;
doubles,z,w;
double*grad,*temp_point;
grad=newdouble[n];
temp_point=newdouble[n];
comput_grad（pf,n,start,grad）;
diver_a=0;
diver_a=diver_a+grad[i-1]*direction[i-1];
do
b=a+step;
for（i=1;
temp_point[i-1]=start[i-1]+b*direction[i-1];
comput_grad（pf,n,temp_point,grad）;
diver_b=0;
diver_b=diver_b+grad[i-1]*direction[i-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;
step=2*step;
}while（diver_b<
=1E-15）;
temp_point[i-1]=start[i-1]+a*direction[i-1];
value_a=（*pf）（temp_point）;
temp_point[i-1]=start[i-1]+b*direction[i-1];
value_b=（*pf）（temp_point）;
s=3*（value_b-value_a）/（b-a）;
z=s-diver_a-diver_b;
w=sqrt（fabs（z*z-diver_a*diver_b））;
//////////////////!
!
t=a+（w-z-diver_a）*（b-a）/（diver_b-diver_a+2*w）;
value_b=（*pf）（temp_point）;
temp_point[i-1]=start[i-1]+t*direction[i-1];
value_t=（*pf）（temp_point）;
diver_t=0;
diver_t=diver_t+grad[i-1]*direction[i-1];
if（diver_t>
1E-6）
b=t;
value_b=value_t;
diver_b=diver_t;
elseif（diver_t<
-1E-6）
a=t;
value_a=value_t;
diver_a=diver_t;
elsebreak;
}while（（fabs（diver_t）>
=1E-6）&
&
（fabs（b-a）>
1E-6））;
delete[]grad;
delete[]temp_point;
return（t）;
//无约束变尺度法DFP函数声明
//
pf指向目标函数的指针，n变量个数，min_point接受初始点、存放结果
极小点处函数值
doubleDFP（
double（*pf）（double*x）,
intn,
double*min_point
）
inti,j;
intk=0;
doublee=1E-5;
doubleg_norm;
double*g0=newdouble[n];
//梯度
double*g1=newdouble[n];
double*dg=newdouble[n];
double*p=newdouble[n];
//搜索方向=-g
doublet;
//一维搜索步长
double*x0=newdouble[n];
double*x1=newdouble[n];
double*dx=newdouble[n];
double**H=newdouble*[n];
for（i=0;
i<
n;
i++）H[i]=newdouble[n];
double**tempH=newdouble*[n];
i++）tempH[i]=newdouble[n];
double*gH=newdouble[n];
double*Hg=newdouble[n];
doublenum1;
doublenum2;
for（i=0;
for（j=0;
j<
j++）
if（i==j）H[i][j]=1.0;
//H0=I
elseH[i][j]=0.0;
tempH[i][j]=0.0;
x0[i]=min_point[i];
comput_grad（pf,n,x0,g0）;
g_norm=0.0;
i++）g_norm=g_norm+g0[i]*g0[i];
g_norm=sqrt（g_norm）;
if（g_norm<
e）
for（i=0;
i++）min_point[i]=x0[i];
delete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
for（i=0;
i++）delete[]H[i];
delete[]H;
i++）delete[]tempH[i];
delete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
i++）p[i]=-g0[i];
t=line_search（pf,n,x0,p）;
i++）x1[i]=x0[i]+t*p[i];
comput_grad（pf,n,x1,g1）;
g_norm=0.0;
i++）g_norm=g_norm+g1[i]*g1[i];
g_norm=sqrt（g_norm）;
//cout<
<
k<
"
"
x0[0]<
x0[1]<
g_norm<
\n"
;
if（g_norm<
{
for（i=0;
i++）min_point[i]=x1[i];
delete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
for（i=0;
delete[]H;
delete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
}
dx[i]=x1[i]-x0[i];
dg[i]=g1[i]-g0[i];
//////////////////求Hk+1的矩阵运算
//g*H,H*g
gH[i]=0.0;
Hg[i]=0.0;
for（j=0;
{
gH[i]=gH[i]+dg[j]*H[j][i];
//Hg[i]=Hg[i]+H[i][j]*dg[j];
Hg[i]=gH[i];
}
//num1,num2
num1=0.0;
num2=0.0;
num1=num1+dx[i]*dg[i];
num2=num2+gH[i]*dg[i];
//tempH[i][j]
tempH[i][j]=0.0;
for（j=0;
{
tempH[i][j]=tempH[i][j]+H[i][j];
tempH[i][j]=tempH[i][j]+dx[i]*dx[j]/num1;
tempH[i][j]=tempH[i][j]-Hg[i]*gH[j]/num2;
}
H[i][j]=tempH[i][j];
/////////////////////////////
//P
i++）p[i]=0.0;
p[i]=p[i]-H[i][j]*g1[j];
}
g0[i]=g1[i];
x0[i]=x1[i];
k=k+1;
}while（g_norm>
e）;
delete[]g0;
delete[]g1;
delete[]dg;
delete[]p;
delete[]x0;
delete[]x1;
delete[]dx;
delete[]H;
delete[]tempH;
delete[]gH;
delete[]Hg;
returnpf（min_point）;
/////////////
doublefun（double*x）
return100*（x[1]-x[0]*x[0]）*（x[1]-x[0]*x[0]）+（1-x[0]）*（1-x[0]）;
voidmain（）
intn=2;
doublemin_point[2]={-5,10};
doublemin_value=DFP（fun,n,min_point）;
cout<
min_point[0]<
and"
min_point[1]<
min_value;
//0.618法线搜索
doubleline_search1（
double（*pf）（double*x）,
intn,
double*start,
double*direction）
intk;
doublel,a,b,c,u,lamda,t,e;
double*xa=newdouble[n];
double*xb=newdouble[n];
double*xc=newdouble[n];
double*xl=newdouble[n];
double*xu=newdouble[n];
//确定初始搜索区间
l=0.001;
a=0;
k=0;
k++;
c=a+l;
for（i=0;
xc[i]=start[i]+c*direction[i];
xa[i]=start[i]+a*direction[i];
l=l/3;
}while（pf（xc）>
=pf（xa））;
//?
?
l=2*l;
b=c+l;
xb[i]=start[i]+b*direction[i];
a=c;
c=b;
}while（pf（xb）<
=pf（xc））;
b=0.1;
//寻优
t=0.618;
e=0.000001;
lamda=b-t*（b-a）;
u=a+t*（b-a）;
xl[i]=start[i]+lamda*direction[i];
xu[i]=start[i]+u*direction[i];
if（pf（xl）<
pf（xu））
b=u;
u=lamda;
lamda=b-t*（b-a）;
else
a=lamda;
lamda=u;
u=t*（b-a）;
}while（b-a>
=e）;
lamda=（a+b）/2;
delete[]xa;
delete[]xb;
delete[]xc;
delete[]xl;
delete[]xu;
returnlamda;
升级会员 