Powell法优化子程序

风花雪月

本程序包含5个C文件
mpowell.c
powell.c
funct.c
jtf.c
hjfgf.c


本程序由heizi友情提供

风花雪月

mpowell.c代码如下：

1. #include "powell.c"
   
2. main()
   
3. {double p[]={1,2};
   
4. double ff,x[2];
   
5. ff=powell(p,0.3,0.001,0.0001,2,x);
   
6.   printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff);
   
7.   getch();
   
8. }

复制代码

powell.c代码如下：

1. #include "hjfgf.c"
   
2. double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
   
3. {double *a,*b,ff;
   
4. a=(double *)malloc(n*sizeof(double));
   
5. b=(double *)malloc(n*sizeof(double));
   
6. jtf(x0,h0,s,n,a,b);
   
7. ff=gold(a,b,epsg,n,x);
   
8. free(a);
   
9. free(b);
   
10. return (ff);
    
11. }
    
12. double powell(double p[],double h0,double eps,double epsg,int n,double x[])
    
13. {int i,j,m;
    
14. double *xx[4],*ss,*s;
    
15. double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
    
16. ss=(double *)malloc(n*(n+1)*sizeof(double));
    
17. s=(double *)malloc(n*sizeof(double));
    
18. for(i=0;i<n;i++)
    
19.   {for(j=0;j<=n;j++)
    
20.    *(ss+i*(n+1)+j)=0;
    
21.    *(ss+i*(n+1)+i)=1;
    
22.   }
    
23.   for(i=0;i<4;i++)
    
24.   xx[i]=(double *)malloc(n*sizeof(double));
    
25.   for(i=0;i<n;i++)
    
26.   *(xx[0]+i)=p[i];
    
27.   for(;;)
    
28.   {for(i=0;i<n;i++)
    
29.     {*(xx[1]+i)=*(xx[0]+i);
    
30.      x[i]=*(xx[1]+i);
    
31.     }
    
32.    f0=f1=objf(x);
    
33.    dlt=-1;
    
34.    for(j=0;j<n;j++)
    
35.     {for(i=0;i<n;i++)
    
36.       {*(xx[0]+i)=x[i];
    
37.        *(s+i)=*(ss+i*(n+1)+j);
    
38.       }
    
39.      f=oneoptim(xx[0],s,h0,epsg,n,x);
    
40.      df=f0-f;
    
41.      if(df>dlt)
    
42.       {dlt=df;
    
43.        m=j;
    
44.       }
    
45.     }
    
46.     sdx=0;
    
47.     for(i=0;i<n;i++)
    
48.      sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
    
49.     if(sdx<eps)
    
50.      {free(ss);
    
51.       free(s);
    
52.       for(i=0;i<4;i++)
    
53.       free(xx[i]);
    
54.       return(f);
    
55.      }
    
56.     for(i=0;i<n;i++)
    
57.      *(xx[2]+i)=x[i];
    
58.     f2=f;
    
59.     for(i=0;i<n;i++)
    
60.     {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
    
61.      x[i]=*(xx[3]+i);
    
62.     }
    
63.     fx=objf(x);
    
64.     f3=fx;
    
65.     q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
    
66.     d=0.5*dlt*(f1-f3)*(f1-f3);
    
67.     if((f3<f1)||(q<d))
    
68.      {if(f2<=f3)
    
69.        for(i=0;i<n;i++)
    
70.         *(xx[0]+i)=*(xx[2]+i);
    
71.       else
    
72.        for(i=0;i<n;i++)
    
73.         *(xx[0]+i)=*(xx[3]+i);
    
74.      }
    
75.     else
    
76.     {for(i=0;i<n;i++)
    
77.      {*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
    
78.       *(s+i)=*(ss+(i+1)*(n+1));
    
79.      }
    
80.    f=oneoptim(xx[0],s,h0,epsg,n,x);
    
81.    for(i=0;i<n;i++)
    
82.     *(xx[0]+i)=x[i];
    
83.    for(j=m+1;j<=n;j++)
    
84.    for(i=0;i<n;i++)
    
85.     *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    
86.     }
    
87.   }
    
88. }

复制代码

风花雪月

funct.c代码如下：

1. #include "stdio.h"
   
2. #include "stdlib.h"
   
3. #include "math.h"
   
4. double objf(double x[])
   
5. {double ff;
   
6. ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
   
7. return(ff);
   
8. }

复制代码

jtf.c代码如下：

1. #include "funct.c"
   
2. void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
   
3. {int i;
   
4. double *x[3],h,f1,f2,f3;
   
5. for(i=0;i<3;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7.   h=h0;
   
8. for(i=0;i<n;i++)
   
9.   *(x[0]+i)=x0[i];
   
10. f1=objf(x[0]);
    
11. for(i=0;i<n;i++)
    
12.   *(x[1]+i)=*(x[0]+i)+h*s[i];
    
13.   f2=objf(x[1]);
    
14. if(f2>=f1)
    
15.   {h=-h0;
    
16.     for(i=0;i<n;i++)
    
17.     *(x[2]+i)=*(x[0]+i);
    
18.    f3=f1;
    
19.     for(i=0;i<n;i++)
    
20.     {*(x[0]+i)=*(x[1]+i);
    
21.      *(x[1]+i)=*(x[2]+i);
    
22.     }
    
23.    f1=f2;
    
24.    f2=f3;
    
25.    }
    
26.    for(;;)
    
27.    {h=2*h;
    
28.      for(i=0;i<n;i++)
    
29.      *(x[2]+i)=*(x[1]+i)+h*s[i];
    
30.    f3=objf(x[2]);
    
31.    if(f2<f3) break;
    
32.    else
    
33.     { for(i=0;i<n;i++)
    
34.        {*(x[0]+i)=*(x[1]+i);
    
35.         *(x[1]+i)=*(x[2]+i);
    
36.        }
    
37.       f1=f2;
    
38.       f2=f3;
    
39.     }
    
40.    }
    
41.    if(h<0)
    
42.     for(i=0;i<n;i++)
    
43.     {a[i]=*(x[2]+i);
    
44.      b[i]=*(x[0]+i);
    
45.     }
    
46.    else
    
47.     for(i=0;i<n;i++)
    
48.     {a[i]=*(x[0]+i);
    
49.      b[i]=*(x[2]+i);
    
50.      }
    
51.    for(i=0;i<3;i++)
    
52.    free(x[i]);
    
53. }

复制代码

hjfgf.c代码如下：

1. #include "jtf.c"
   
2. double gold(double a[],double b[],double eps,int n,double xx[])
   
3. {int i;
   
4. double f1,f2,*x[2],ff,q,w;
   
5. for(i=0;i<2;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7. for(i=0;i<n;i++)
   
8.   {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
   
9.    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
   
10.   }
    
11.   f1=objf(x[0]);
    
12.   f2=objf(x[1]);
    
13.   do
    
14.    {if(f1>f2)
    
15.      {for(i=0;i<n;i++)
    
16.       {b[i]=*(x[0]+i);
    
17.        *(x[0]+i)=*(x[1]+i);
    
18.        }
    
19.      f1=f2;
    
20.      for(i=0;i<n;i++)
    
21.       *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
    
22.      f2=objf(x[1]);
    
23.      }
    
24.     else
    
25.      { for(i=0;i<n;i++)
    
26.        {a[i]=*(x[1]+i);
    
27.        *(x[1]+i)=*(x[0]+i);}
    
28.      f2=f1;
    
29.     for(i=0;i<n;i++)
    
30.      *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    
31.     f1=objf(x[0]);
    
32.      }
    
33.   q=0;
    
34.   for(i=0;i<n;i++)
    
35.    q=q+(b[i]-a[i])*(b[i]-a[i]);
    
36.   w=sqrt(q);
    
37.   }while(w>eps);
    
38.   for(i=0;i<n;i++)
    
39.    xx[i]=0.5*(a[i]+b[i]);
    
40.   ff=objf(xx);
    
41.   for(i=0;i<2;i++)
    
42.   free(x[i]);
    
43.   return(ff);
    
44. }

复制代码

lzjld502

xiexiele hehe:@)

zhoyi

注释到没有，很难看明白。要是看明白了，也会编写了 。建议上传时把注释也写上。不管怎样，还是谢谢louzhu。

hnsqzrj

怎么不能运行

yangwenli

感谢楼主提供的程序:@D

to楼上:编译在一起就可运行:)

crywolf911

问题是怎么编译在一起啊，我是新手啊~！

modest_xiu

谢谢~~~楼主辛苦