?? cpp1.cpp
字號:
#include<iostream.h>
#include<math.h>
#include<time.h>
#define N1 40
//cF54 DIXMAANE
//*
//* Initial point x0=[2.,2.,2...,2.].
//*
void CALF(int n,double x[N1],double &f){
double alpha,beta,gamma,delta;
int k1,k2,k3,k4;
alpha=1.0;
beta=0.0;
gamma=0.1250;
delta=0.1250;
k1=1;
k2=0;
k3=0;
k4=1;
int m = n/3;
f=1.0;
for(int i=0;i<n;i++)
f=f+alpha*x[i]*x[i]*pow((float(i+1)/float(n)),k1);
for(i=0;i<n-1;i++)
f=f+beta*x[i]*x[i]*(x[i+1]+x[i+1]*x[i+1])*pow((float(i+1)/float(n)),k2);
for(i=0;i<2*m;i++)
f=f+gamma*x[i]*x[i]*(pow(x[i+m],4))*pow((float(i+1)/float(n)),k3);
for(i=0;i<m;i++)
f=f+delta*x[i]*x[i+2*m]*pow((float(i+1)/float(n)),k4);
}
//*-- Gradient
void CALG(int n,double x[N1],double g[N1]){
double alpha,beta,gamma,delta;
int k1,k2,k3,k4;
alpha=1.0;
beta=0.0;
gamma=0.1250;
delta=0.1250;
k1=1;
k2=0;
k3=0;
k4=1;
int m=n/3;
g[0]=2.0*alpha*x[0]*pow((float(1)/float(n)),k1)+2.0*beta*x[0]*pow((x[0]+x[1]*x[1]),2)*pow((float(1)/float(n)),k2)+2.0*gamma*x[0]*pow(x[m],4)*pow((float(1)/float(n)),k3)+delta*x[2*m]*pow((float(1)/float(n)),k4);
for(int i=1;i<m;i++)
g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+delta*x[i+2*m]*pow((float(i+1)/float(n)),k4);
for(i=m;i<2*m;i++)
g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)*pow((float(i-m+1)/float(n)),k3);
for(i=2*m;i<n-1;i++)
g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i)/float(n)),k2)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)* pow((float(i-m+1)/float(n)),k3)+delta*x[i-2*m]*pow((float(i-2*m+1)/float(n)),k4);
g[n-1]=2.0*alpha*x[n-1]+2.0*beta*x[n-2]*x[n-2]*(x[n-1]+x[n-1]*x[n-1])*(1.0+2.0*x[n-1])*pow((float(n-1)/float(n)),k2)+4.0*gamma*x[2*m-1]*x[2*m-1]*x[n-1]*x[n-1]*x[n-1]*pow((float(2*m)/float(n)),k3)+delta*x[m-1]*pow((float(m)/float(n)),k4);
}
int conclude1(double G[N1])
{double a=0.00;
for(int i=0;i<N1;i++)
{a+=G[i]*G[i];}
if(sqrt(a)>=0.0001) {return 1;}
else return 0;
}
double DDOT1(double p[N1],double Q[N1])
{double t=0.00;
for(int i=0;i<N1;i++)
{t+=p[i]*Q[i];}
return t;
}
void main()
{clock_t start,end;
start=clock();
double k=1,kesi=0.0001;
double beta,m=0.00,n=0.00;
double X[N1],G[N1],WORK[N1],D[N1],P[N1],G1[N1] ;
double F,ALPHA,IPRINT=1,C1=0.01,C2=0.05;
int NLSMAX ;
double FO ;
int I;
double DGO;
double A1;
double F1;
double DG1;
double A2;
double A21;
double RATIO1;
double RATIO2;
int NCALFN,NCALG,NCALF;
double DD,XX;
double TRY,F2,DG,DG2;
double T1,T2,T3,T4;
double INFO=1;
for(int i=0;i<N1;i++)
{X[i]=2;}
CALG(N1,X,G);
for(i=0;i<N1;i++)
{G1[i]=G[i];}
for(i=0;i<N1;i++)
{P[i]=-G[i];}
//................................................................................
step1: CALG(N1,X,G);
//cout<<x[1]<<'\t';
for(i=0;i<N1;i++)
{D[i]=-G[i];}
step2: if(conclude1(G)){
if(k==1)
{beta=0;
for(i=0;i<N1;i++)
{ G1[i]=G[i];}
}
else {for(i=0;i<N1;i++)
{m+=G[i]*G[i];
n+=G1[i]*G1[i];
}
beta=m/n;
for(i=0;i<N1;i++)
{ G1[i]=G[i];}
}
for(i=0;i<N1;i++)
{P[i]=-G[i]+beta*P[i];
}
goto step3;
}
else goto step7;
//========= step 1 設初始值 ===========/
step3: NLSMAX=50;
FO=F;
for( I=0;I<N1;I++)
WORK[I]=X[I];
DGO=DDOT1(D,G);
A1=0.0;
F1=F;
DG1=DGO;
A2=-1.0;
A21=0;
RATIO1=1.0;
RATIO2=0.1;
NCALFN=0,NCALG=1,NCALF=1;
ALPHA=1.0;
//INFO;
if(DGO>=0.0) goto enter_1100;
if(IPRINT!=0) goto enter_40;
cout<<"ENTER LINE SEARCH"<<endl;
////////////////////////////////////////////////////
//========== step 2 計算新的試驗點 ===========//
enter_40 : for(I=0;I<N1;I++)
X[I]=WORK[I]+ALPHA*D[I];
CALF(N1,X,F);//目標函數值
NCALFN=NCALFN+1;
if(NCALFN>NLSMAX) goto enter_1300;
if(F<=FO+C1*ALPHA*DGO) goto enter_200;
//========== Step 3 如果第一個等式不成立========//
//========== 通過二次修正降低步長 ==========//
RATIO1=RATIO1*0.1;
if(DGO/(1.0+fabs(FO))<-1.0e-16) goto enter_80;
DD=DDOT1(D,D);
XX=DDOT1(X,X);
if(ALPHA*sqrt(DD)/(1.0+sqrt(XX))<=1.0e-12) goto enter_1400;
enter_80 : A2=ALPHA;
F2=F;
A21=0;
enter_90 : TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
TRY=-DG1/TRY/2.0;
goto enter_240;
//========= Step 4 如果第一個滿足測試第二個不等式=========//
enter_200 : CALG(N1,X,G);//函數的梯度
NCALG=NCALG+1;
DG=DDOT1(D,G);
if(fabs(DG)<=-C2*DGO) goto enter_1500;
if(DG>0.0) goto enter_210;
if(A2<0.0) goto enter_260;
RATIO2=RATIO2*0.1;
A1=ALPHA;
F1=F;
DG1=DG;
RATIO1=0.1;
if(!A21) goto enter_90;
goto enter_220;
//========== Step 5 使用三次插值尋找===========//
// 在[A1,A2]內的新的點 //
enter_210 : RATIO1=0.1;
A2=ALPHA;
F2=F;
DG2=DG;
A21=1;
enter_220 : T1=(F2-F1)/(A2-A1)-DG1;
T2=DG2-DG1;
T3=3.0*T1-T2;
T4=T2-T1-T1;
TRY=T3*T3-3.0*T4*DG1;
if(TRY<0.0) goto enter_230;
if(sqrt(TRY)<=T3) goto enter_230;
TRY=-DG1*(A2-A1)/(sqrt(TRY)+T3);
goto enter_240;
enter_230 : TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
TRY=-DG1/TRY/2.0;
//======== 如果必要去掉這步==========//
enter_240 : if(TRY<RATIO1*(A2-A1)) TRY=RATIO1*(A2-A1);
if(TRY>(1.0-RATIO2)*(A2-A1)) TRY=(1.0-RATIO2)*(A2-A1);
ALPHA=TRY+A1;
goto enter_40;
//========= Step 6 ALPHA是最大的可行點,通過使用三次插值找到另一個可行的值使STEP>ALPHA=======//
enter_260 : T1=(F-F1)/(ALPHA-A1)-DG;
T2=DG1-DG;
T3=3.0*T1-T2;
T4=T2-T1-T1;
TRY=T3*T3-3.0*T4*DG;
if(TRY<0.0) goto enter_275;
if(sqrt(TRY)<=T3) goto enter_275;
TRY=-DG*(ALPHA-A1)/(sqrt(TRY)-T3);
if(TRY<=1.0*10*(ALPHA-A1)) goto enter_280;
enter_275 : TRY=1.0*10*(ALPHA-A1);
enter_280 : if(A2<0.0) goto enter_290;
if(TRY<0.1*(A2-A1)) TRY=0.1*(A2-A1);
if(TRY>0.9*(A2-A1)) TRY=0.9*(A2-A1);
enter_290 : A1=ALPHA;
F1=F;
DG1=DG;
ALPHA=A1+TRY;
goto step6;
goto enter_40;
//========== 使用上升方向搜索或者參量是不正確的 ===========//
enter_1100 : //cout<<"AN UPHILL SEARCH DIRECTION IS USED!"<<endl;
INFO=1;
goto enter_1500;
enter_1300 : //cout<<"LINE SEARCH FAILS AFTER NLSMAX TRIAL STEPSIZES!"<<endl;
INFO=2;
goto enter_1500;
enter_1400 : CALG(N1,X,G);
NCALG=NCALG+1;
//cout<<"STEP IS TOO SMALL!"<<endl;
INFO=3;
enter_1500 : NCALF=NCALF+NCALFN;
//return;
step6:
k=k+1;
goto step1;
step7: CALF(N1,X,F);
cout<<"問題的規模:"<<N1<<endl<<"迭代次數:"<<k<<endl<<"精度為:"<<kesi<<endl<<"F="<<F<<endl;
end=clock();
cout<<"計算用時:"<<(double)(end-start)/CLOCKS_PER_SEC<<"ms"<<endl;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -