conjugategradientmethod.txt

使用黄金分割法和进退法进行区间搜索的共轭梯度法。

#include<stdio.h>
#include<math.h>
#define N 10
#define eps pow(10,-6)

double f(double x[],double p[],double t)
{
double s;
s=pow(x[0]+t*p[0],2)+25*pow(x[1]+t*p[1],2);
return s;
}

/*以下是进退法搜索区间源程序*/
void sb(double *a,double *b,double x[],double p[])
{
double t0,t1,t,h,alpha,f0,f1;
int k=0;
t0=2.5;      /*初始值*/
h=1;        /*初始步长*/
alpha=2;   /*加步系数*/

f0=f(x,p,t0);
t1=t0+h;
f1=f(x,p,t1);
while(1)
{
   if(f1<f0)
   {
    h=alpha*h; t=t0;
    t0=t1;   f0=f1;
    k++;
   }
   else
   {
    if(k==0)
    {h=-h;t=t1;}
    else
    {   
     *a=t<t1?t:t1;
     *b=t>t1?t:t1;
     break;
    }
   }
   t1=t0+h;
   f1=f(x,p,t1);
}
}
/*以下是黄金分割法程序源代码*/
double hjfg(double x[],double p[])
{
double beta,t1,t2,t;
double f1,f2;
double a=0,b=0;
double *c,*d;
c=&a,d=&b;
sb(c,d,x,p);/*调用进退法搜索区间*/
printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf",x[0],x[1],p[0],p[1]);
printf("\n[a,b]=[%lf,%lf]",a,b);
beta=(sqrt(5)-1.0)/2;
t2=a+beta*(b-a); f2=f(x,p,t2);
t1=a+b-t2;      f1=f(x,p,t1);
while(1)
{
   if(fabs(t1-t2)<eps)
   break;
   else 
   {
    if(f1<f2)
    {
     t=(t1+t2)/2;
     b=t2;     t2=t1;
       f2=f1;     t1=a+b-t2;
             f1=f(x,p,t1); 
    }
      else
      {
      a=t1;     t1=t2;
       f1=f2;
       t2=a+beta*(b-a);
       f2=f(x,p,t2);
      }
   }
}
t=(t1+t2)/2;
return t;
}
/*以下是共轭梯度法程序源代码*/
void gtd()
{
double x[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0;
int i,k,n;
printf("请输入函数的元数值n=");
scanf("%d",&n);
printf("\n请输入初始值:\n");
for(i=0;i<n;i++)
scanf("%lf",&x[i]);
f0=f(x,g,t);
g[0]=2*x[0]; g[1]=50*x[1];
mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/
if(mod1>eps)
{
   p[0]=-g[0]; p[1]=-g[1]; k=0;

   while(1)
   {
    t=hjfg(x,p);/*调用黄金分割法求t的值*/
    printf("\np1=%lf,p2=%lf,t=%lf",p[0],p[1],t);
   
    x[0]=x[0]+t*p[0]; x[1]=x[1]+t*p[1];
    g[0]=2*x[0];       g[1]=50*x[1];
         /*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf",x[0],x[1],g[0],g[1]);*/
         mod2=sqrt(pow(g[0],2)+pow(g[1],2)); /*求梯度的长度*/
        
         if(mod2<=eps)    break;
         else
         {
         if(k+1==n)
         { 
            g[0]=2*x[0]; g[1]=50*x[1];
            p[0]=-g[0]; p[1]=-g[1]; k=0;
         }
         else
         {
            nanda=pow(mod2,2)/pow(mod1,2);
            printf("\nnanda=%lf,mod=%lf",nanda,mod2);
            p[0]=-g[0]+nanda*p[0];
            p[1]=-g[1]+nanda*p[1];
            mod1=mod2;
            k++;
         }
         }
         printf("\n--------------------------");
   }
}
printf("\n最优解为x1=%lf,x2=%lf",x[0],x[1]);
printf("\n最终的函数值为%lf",f(x,g,t));
}
main()
{
gtd();
}

运行结果如下：

请输入函数的元数值n=2

请输入初始值:
2 2

x1=2.000000,x2=2.000000,p1=-4.000000,p2=-100.000000
[a,b]=[-4.500000,1.500000]
p1=-4.000000,p2=-100.000000,t=0.020030
nanda=0.001474,mod=3.842730
--------------------------
x1=1.919879,x2=-0.003022,p1=-3.845655,p2=0.003665
[a,b]=[-4.500000,1.500000]
p1=-3.845655,p2=0.003665,t=0.499240
--------------------------
x1=-0.000026,x2=-0.001192,p1=0.000052,p2=0.059610
[a,b]=[-4.500000,1.500000]
p1=0.000052,p2=0.059610,t=0.020000
nanda=0.000000,mod=0.000050
--------------------------
x1=-0.000025,x2=-0.000000,p1=0.000050,p2=0.000001
[a,b]=[-4.500000,1.500000]
p1=0.000050,p2=0.000001,t=0.495505
--------------------------
x1=-0.000000,x2=0.000000,p1=0.000000,p2=-0.000023
[a,b]=[-4.500000,1.500000]
p1=0.000000,p2=-0.000023,t=0.020007
最优解为x1=-0.000000,x2=-0.000000
最终的函数值为0.000000