DFP变尺度法优化子程序
本程序包含5个C文件mdfp.c
dfpopt.c
funct.c
jtf.c
hjfgf.c
本程序由heizi友情提供 mdfp.c代码如下：
#include "dfpopt.c"
main()
{double x[]={1,2};
double ff;
ff=dfpopt(x,0.2,0.0001,0.000001,2);
printf("\nx=%f,x=%f,ff=%f",x,x,ff);
getch();
}
dfpopt.c代码如下：
#include "hjfgf.c"
#include "stdlib.h"
void gradient(double x[],double g[],int n)
{int i;
double af,f1,f2,dltx=0.000001;
for(i=0;i<n;i++)
g=0;
f1=objf(x);
for(i=0;i<n;i++)
{af=*(x+i);
*(x+i)=af+dltx;
f2=objf(x);
g=(f2-f1)/dltx;
*(x+i)=af;
}
}
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
a=(double *)malloc(n*sizeof (double ));
b=(double *)malloc(n*sizeof (double ));
jtf(x0,h0,s,n,a,b);
ff=gold(a,b,epsg,n,x);
free(a);
free(b);
return(ff);
}
double dfpopt(double xx[],double h0,double eps,double epsg,int n)
{int i,j,k;
double ae,zcc;
double *s,*x,*ay,*df,*zd,*zc,*zh;
s=(double *)malloc(n*sizeof (double ));
x=(double *)malloc(n*sizeof (double ));
for(i=0;i<2;i++)
{ay=(double*)malloc(n*sizeof(double));
df=(double*)malloc(n*sizeof(double));
zd=(double*)malloc(n*sizeof(double));
zc=(double*)malloc(n*sizeof(double));
zh=(double*)malloc(n*n*sizeof(double));
}
for(i=0;i<n;i++)
*(ay+i)=xx;
L1:
k=0;
for (i=0;i<n;i++)
for (j=0;j<n;j++)
{*(zh+i*n+j)=0;
if(i==j)
*(zh+i*n+j)=1;
}
for(i=0;i<n;i++)
{*(x+i)=*(ay+i);
*(ay+i)=*(ay+i);
}
gradient(x,df,n);
for(i=0;i<n;i++)
{*(s+i)=0;
*(df+i)=*(df+i);
for(j=0;j<n;j++)
*(s+i)=*(s+i)-(*(zh+i*n+j))*(*(df+j));
}
L2:
oneoptim(x,s,h0,epsg,n,ay);
for(i=0;i<n;i++)
*(x+i)=*(ay+i);
gradient(x,df,n);
ae=0;
for(i=0;i<n;i++)
ae=ae+(*(df+i))*(*(df+i));
if(ae<=eps)
{for(i=0;i<n;i++)
*(xx+i)=*(x+i);
free(s);
free(x);
for(i=0;i<2;i++)
{free(ay);
free(df);
free(zd);
free(zc);
free(zh);
}
return(objf(xx));
}
if(k==n) goto L1;
zcc=0;
for(i=0;i<n;i++)
{*(zd+i)=*(ay+i)-(*(ay+i));
*(zd+i)=*(df+i)-(*(df+i));
*(df+i)=*(df+i);
zcc=zcc+(*(zd+i))*(*(zd+i));
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
*(zh+i*n+j)=*(zd+i)*(*(zd+j))/zcc;
for(i=0;i<n;i++)
{ *(zc+i)=0;
for(j=0;j<n;j++)
*(zc+i)=*(zc+i)+(*(zd+j))*(*(zh+j*n+i));
}
zcc=0;
for(i=0;i<n;i++)
zcc=zcc+(*(zc+i))*(*(zd+i));
for(i=0;i<n;i++)
{*(zc+i)=0;
*(zc+i)=0;
for (j=0;j<n;j++)
{*(zc+i)=*(zc+i)+(*(zh+i*n+j))*(*(zd+j));
*(zc+i)=*(zc+i)+(*(zh+j))*(*(zh+j*n+i));
}
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
*(zh+i*n+j)=*(zh+i*n+j)+(*(zh+i*n+j))-(*(zc+i))*(*(zc+j))/zcc;
for(i=0;i<n;i++)
{*(s+i)=0;
for(j=0;j<n;j++)
*(s+i)=*(s+i)-(*(zh+i*n+j))*(*(df+j));
}
k=k+1;
goto L2;
} funct.c代码如下：
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
double objf(double x[])
{double ff;
ff=x*x+x*x-x*x-10*x-4*x+60;
return(ff);
}
jtf.c代码如下：
#include "funct.c"
void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
{int i;
double *x,h,f1,f2,f3;
for(i=0;i<3;i++)
x=(double *)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x+i)=x0;
f1=objf(x);
for(i=0;i<n;i++)
*(x+i)=*(x+i)+h*s;
f2=objf(x);
if(f2>=f1)
{h=-h0;
for(i=0;i<n;i++)
*(x+i)=*(x+i);
f3=f1;
for(i=0;i<n;i++)
{*(x+i)=*(x+i);
*(x+i)=*(x+i);
}
f1=f2;
f2=f3;
}
for(;;)
{h=2*h;
for(i=0;i<n;i++)
*(x+i)=*(x+i)+h*s;
f3=objf(x);
if(f2<f3) break;
else
{ for(i=0;i<n;i++)
{*(x+i)=*(x+i);
*(x+i)=*(x+i);
}
f1=f2;
f2=f3;
}
}
if(h<0)
for(i=0;i<n;i++)
{a=*(x+i);
b=*(x+i);
}
else
for(i=0;i<n;i++)
{a=*(x+i);
b=*(x+i);
}
for(i=0;i<3;i++)
free(x);
}
hjfgf.c代码如下：
#include "jtf.c"
double gold(double a[],double b[],double eps,int n,double xx[])
{int i;
double f1,f2,*x,ff,q,w;
for(i=0;i<2;i++)
x=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*(x+i)=a+0.618*(b-a);
*(x+i)=a+0.382*(b-a);
}
f1=objf(x);
f2=objf(x);
do
{if(f1>f2)
{for(i=0;i<n;i++)
{b=*(x+i);
*(x+i)=*(x+i);
}
f1=f2;
for(i=0;i<n;i++)
*(x+i)=a+0.382*(b-a);
f2=objf(x);
}
else
{ for(i=0;i<n;i++)
{a=*(x+i);
*(x+i)=*(x+i);}
f2=f1;
for(i=0;i<n;i++)
*(x+i)=a+0.618*(b-a);
f1=objf(x);
}
q=0;
for(i=0;i<n;i++)
q=q+(b-a)*(b-a);
w=sqrt(q);
}while(w>eps);
for(i=0;i<n;i++)
xx=0.5*(a+b);
ff=objf(xx);
for(i=0;i<2;i++)
free(x);
return(ff);
} xiexie xiexielehehe:@D thank you 谢谢分享
页:
[1]