sousuo.cpp

一维搜索法求极值的问题

#include<math.h>
#include<stdio.h>
#include<conio.h>

double f(double);
double golden_split(double(*)(),double,double);/*黄金分割法*/
double rcczf(double(*)(),double[],double[],double); /*二次插值法*/
void jtf(double(*)(),double,double *,double *,double *);/*进退法*/

void jtf(double(*objf)(),double a0,double *h0,double *a,double *y)
{
double h;
a[0]=a0;
h=*h0;
a[1]=a[0]+h,y[1]=objf(a[1]);
if(y[1]>y[0])
{
   h=-h;
   a[2]=a[0],y[2]=y[0];
   a[0]=a[1],y[0]=y[1];
   a[1]=a[2],y[1]=y[2];
}
a[2]=a[1]+h;
y[2]=objf(a[2]);
while(y[2]<y[1])
{
   h=2*h;
   a[0]=a[1],y[0]=y[1];
   a[1]=a[2],y[1]=y[2];
   a[2]=a[1]+h;
   y[2]=objf(a[2]);
}
if(a[2]<a[0])
{
   double temp;
   temp=a[2],a[2]=a[0],a[0]=temp;
      temp=y[2],y[2]=y[0],y[0]=temp;
}
    *h0=h; /*返回最后步长*/
}

double golden_split(double(*f)(),double a,double b)
{
    #define lamda 0.6180339
    #define eps 1.0e-8
double a1,a2,y1,y2;
a1=b-lamda*(b-a),y1=f(a1);
a2=a+lamda*(b-a),y2=f(a2);
while(fabs(b-a)>eps&&fabs(y1-y2)>eps)
{
   if(y1>=y2)
   {
     a=a1;
     a1=a2,y1=y2;
     a2=a+lamda*(b-a);
     y2=f(a2);
   }
   else
   {
    b=a2;
    a2=a1,y2=y1;
    a1=b-lamda*(b-a);
    y1=f(a1);
   }
}
return 0.5*(a+b);
}


double rcczf(double (*f)(),double a[3],double y[3],double h)
{
double c1,c2,ap,yp;
double accu=1.0e-6 ;
a[1]=0.4*(a[0]+a[2]);
while(1)
{
   c1=(y[2]-y[0])/(a[2]-a[0]);
   c2=((y[1]-y[0])/(a[1]-a[0])-c1)/(a[1]-a[2]);
   ap=0.5*(a[0]+a[2]-c1/c2);
   yp=f(ap);

   if(fabs(y[1]-yp)<accu)
   {
    if(y[1]<yp) return a[1];
    else return ap;
   }
   else
   {
    if((ap-a[1])*h>0)
    {
     if(y[1]>=ap){a[0]=a[1],y[0]=y[1],a[1]=ap,y[1]=yp;}
     else {a[2]=ap,y[2]=yp;}
    }
    else
    {
     if(y[1]>=yp){a[2]=a[1],y[2]=y[1],a[1]=ap,y[1]=yp;}
     else {a[0]=ap,y[0]=yp;}
    }
   }
}
}


double fun(double x)
{
double y;
y=x*x-2*x+5;
return(y);
}

void main()
{
   double X,Y;
   double X1,Y1;
   double a0,h0,a[3],y[3];
   int i;
   a0=0,h0=1;
   for(i=1;i<3;i++) a[i]=0,y[i]=0;
   jtf(fun,a0,&h0,a,y);
   X=golden_split(fun,a[0],a[2]);
   Y=fun(X);
   printf("\nthe otimial value is X=%10.4f,\tY=%10.4f\n\n",X,Y);
  
   a0=0.1,h0=1;
   for(i=1;i<3;i++) a[i]=0,y[i]=0;
   jtf(fun,a0,&h0,a,y);
   X1=rcczf(fun,a,y,h0);
   Y1=fun(X1);
   printf("\n\nthe otimial value is X1=%10.4f,\tY1=%10.4f\n",X1,Y1);

   printf("Press any key to end up!\n");
   getch();
}