机械优化设计实验c语言程序

二次插值法
#include "stdio.h"
#include "math.h"
#include "conio.h"
void main()
{
float *area(float a1,float p,float a[3]);
float f(float x);
float ar,fr;
float a2,a3;
float f1,f2,f3;
float a1=10,p=0.01,e=0.00001;
float pa[3];
area(a1,p,pa);
a1=pa[0];
a2=pa[1];
a3=pa[2];
f1=f(a1);
f2=f(a2);
f3=f(a3);
do
{
ar=((a3*a3-a2*a2)*f1+(a1*a1-a3*a3)*f2+(a2*a2-a1*a1)*f3); ar=ar/2/((a3-a2)*f1+(a1-a3)*f2+(a2-a1)*f3);
fr=f(ar);
if(ar>a2)
{if(fr>f2)
{a3=ar;f3=fr;}
else if(fr<f2)
{a1=a2;f1=f2;
a2=ar;f2=fr;}
else
{a3=ar;a1=a2;a2=(a1+a3)/2;
f1=f2;f3=fr;f2=f(a2);}
}
else if(ar<a2)
{if(fr>f2)
{a1=ar;f1=fr;}
else if(fr<f2)
{a3=a2;f3=f2;
a2=ar;f2=fr;}
else
{a1=ar;a3=a2;a2=(a1+a3)/2;
f1=fr;f3=f2;f2=f(a2);}
}
if(fabs(a1-a3)<=e) break;
}while(1);
if(f2<fr)
{ar=a2;fr=f2;}
printf("\nx*=%f\nf(x*)=%f",ar,fr);
}
float *area(float a1,float p,float a[3])
{float f(float x);
float a2,f2,a3,f3,temp;
float acc=0.01;
float f1=f(a1);
float p;
while(1)
{a2=a1+p;f2=f(a2);
if(f2>=f1)
{if(fabs(f2-f1)<acc)
p=p/2;
else
p=-p;
}
else break;
}
while(1)
{a3=a2+p;f3=f(a3);
if(f2<=f3) break;
p=2*p;
a1=a2;f1=f2;
a2=a3;f2=f3;
}
if(a1>a3)
{temp=a1;a1=a3;a3=temp;}
a[0]=a1;a[1]=a2;a[2]=a3;
return a;
}
float f(float x)
{float y=pow(x,4)-4*pow(x,3)-6*pow(x,2)-16*x+4;
return y;
}
变尺度法
计算f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2 的无约束极值，初始点x0=[1,1]。

/*
tt ---- 一维搜索初始步长
ff ---- 差分法求梯度时的步长
ac ---- 终止迭代收敛精度
ad ---- 一维搜索收敛精度
n ----- 设计变量的维数
xk[n] -- 迭代初始点
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>
#define tt 0.01
#define ff 1.0e-6
#define ac 1.0e-6
#define ad 1.0e-6
#define n 2
double ia;
double fny(double *x)
{
double x1=x[0],x2=x[1];
double f;
f=x1*x1+2*x2*x2-4*x1-2*x1*x2;
return f;
}
double * iterate(double *x,double a,double *s) {
double *x1;
int i;
x1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
x1[i]=x[i]+a*s[i];
return x1;
}
double func(double *x,double a,double *s) {
double *x1;
x1=iterate(x,a,s);
f=fny(x1);
return f;
}
void finding(double a[3],double f[3],double *xk,double *s) {
double t=tt;
int i;
double a1,f1;
a[0]=0;f[0]=func(xk,a[0],s);
for(i=0;;i++)
{
a[1]=a[0]+t;
f[1]=func(xk,a[1],s);
if(f[1]<f[0]) break;
if(fabs(f[1]-f[0])>=ad)
{
t=-t;
a[0]=a[1];f[0]=f[1];
}
else
{
if(ia==1) return; //break
t=t/2;ia=1;
}
}
for(i=0;;i++)
{
a[2]=a[1]+t;
f[2]=func(xk,a[2],s);
if(f[2]>f[1]) break;
t=2*t;
a[0]=a[1];f[0]=f[1];
a[1]=a[2];f[1]=f[2];
}
if(a[0]>a[2])
{
a1=a[0];
f1=f[0];
a[0]=a[2];
f[0]=f[2];
f[2]=f1;
}
return;
}
double lagrange(double *xk,double *ft,double *s)
{
int i;
double a[3],f[3];
double b,c,d,aa;
finding(a,f,xk,s);
for(i=0;;i++)
{
if(ia==1) { aa=a[1]; *ft=f[1]; break; }
d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]);
if(fabs(d)==0) break;
c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d;
if(fabs(c)==0) break;
b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]);
aa=-b/(2*c);
*ft=func(xk,aa,s);
if(fabs(aa-a[1])<=ad) {if(*ft>f[1]) aa=a[1];break;}
if(aa>a[1])
{
if(*ft>f[1]) {a[2]=aa;f[2]=*ft;}
else if(*ft<f[1]) {a[0]=a[1];a[1]=aa;f[0]=f[1];f[1]=*ft;}
else if(*ft==f[1])
{
a[2]=aa;a[0]=a[1];
f[2]=*ft;f[0]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
else
{
if(*ft>f[1]) {a[0]=aa;f[0]=*ft;}
else if(*ft<f[1]) {a[2]=a[1];a[1]=aa;f[2]=f[1];f[1]=*ft;}
else if(*ft==f[1])
{a[0]=aa;a[2]=a[1];
f[0]=*ft;f[2]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
}
if(*ft>f[1]) {*ft=f[1];aa=a[1];}
return aa;
}
double *gradient(double *xk)
{
double *g,f1,f2,q;
int i;
g=(double*)malloc(n*sizeof(double));
f1=fny(xk);
for(i=0;i<n;i++)
{q=ff;
xk[i]=xk[i]+q; f2=fny(xk);
g[i]=(f2-f1)/q; xk[i]=xk[i]-q;
}
return g;
}
double * bfgs(double *xk)
{
double u[n],v[n],h[n][n],dx[n],dg[n],s[n]; double aa,ib;
double *ft,*xk1,*g1,*g2,*xx,*x0=xk; double fi;
int i,j,k;
ft=(double *)malloc(sizeof(double));
xk1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{
s[i]=0;
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
g1=gradient(xk);
fi=fny(xk);
x0=xk;
for(k=0;k<n;k++)
{
ib=0;
if(ia==1) { xx=xk; break; }
ib=0;
for(i=0;i<n;i++) s[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
s[i]+= -h[i][j]*g1[j];
aa=lagrange(xk,ft,s);
xk1=iterate(xk,aa,s);
g2=gradient(xk1);
for(i=0;i<n;i++)
if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))
{ib=ib+1;}
if(ib==0) { xx=xk1; break; }
fi=*ft;
if(k==n-1)
{ int j;
xk=xk1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
g1=g2; k=-1;
}
else
{
int j;
double a1=0,a2=0;
for(i=0;i<n;i++)
{
dg[i]=g2[i]-g1[i];
dx[i]=xk1[i]-xk[i];
for(i=0;i<n;i++)
{
int j;
u[i]=0;v[i]=0;
for(j=0;j<n;j++)
{
u[i]=u[i]+dg[j]*h[j][i];
v[i]=v[i]+dg[j]*h[i][j];
}
}
for(j=0;j<n;j++)
{
a1+=dx[j]*dg[j];
a2+=v[j]*dg[j];
}
if(fabs(a1)!=0)
{
a2=1+a2/a1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1;
}
xk=xk1; g1=g2;
}
}
if(*ft>fi) { *ft=fi; xx=xk;}
xk=x0;
return xx;
}
void main ()
{
int k;
double *xx,f;
double xk[n]={1,1};
xx=bfgs(xk);
f=fny(xx);
printf("\n\nThe Optimal Design Result Is:\n");
for(k=0;k<n;k++)
{printf("\n\tx[%d]*=%f",k+1,xx[k]);}
printf("\n\tf*=%f",f);
getch();
}
鲍威尔法程序
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
double objf(double x[])
{double ff;
ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
return(ff);
}
void jtf(double x0[],double h0,double s[],int n,double a[],double b[]) {int i;
double *x[3],h,f1,f2,f3;
for(i=0;i<3;i++)
x[i]=(double *)malloc(n*sizeof(double));
h=h0;
for(i=0;i<n;i++)
*(x[0]+i)=x0[i];
f1=objf(x[0]);
for(i=0;i<n;i++)
*(x[1]+i)=*(x[0]+i)+h*s[i];
f2=objf(x[1]);
if(f2>=f1)
{h=-h0;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[0]+i);
f3=f1;
for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
}
f1=f2;
f2=f3;
}
for(;;)
{h=2*h;
for(i=0;i<n;i++)
*(x[2]+i)=*(x[1]+i)+h*s[i];
f3=objf(x[2]);
if(f2<f3) break;
else
{ for(i=0;i<n;i++)
{*(x[0]+i)=*(x[1]+i);
*(x[1]+i)=*(x[2]+i);
}
f1=f2;
f2=f3;
}
}
if(h<0)
for(i=0;i<n;i++)
{a[i]=*(x[2]+i);
b[i]=*(x[0]+i);
}
else
for(i=0;i<n;i++)
{a[i]=*(x[0]+i);
b[i]=*(x[2]+i);
}
for(i=0;i<3;i++)
free(x[i]);
}
double gold(double a[],double b[],double eps,int n,double xx[]) {int i;
double f1,f2,*x[2],ff,q,w;
for(i=0;i<2;i++)
x[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
}
f1=objf(x[0]);
f2=objf(x[1]);
do
{if(f1>f2)
{for(i=0;i<n;i++)
{b[i]=*(x[0]+i);
*(x[0]+i)=*(x[1]+i);
}
f1=f2;
for(i=0;i<n;i++)
*(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
f2=objf(x[1]);
}
else
{ for(i=0;i<n;i++)
{a[i]=*(x[1]+i);
*(x[1]+i)=*(x[0]+i);}
f2=f1;
for(i=0;i<n;i++)
*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
f1=objf(x[0]);
}
q=0;
for(i=0;i<n;i++)
q=q+(b[i]-a[i])*(b[i]-a[i]);
w=sqrt(q);
}while(w>eps);
for(i=0;i<n;i++)
xx[i]=0.5*(a[i]+b[i]);
ff=objf(xx);
for(i=0;i<2;i++)
free(x[i]);
return(ff);
}
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 powell(double p[],double h0,double eps,double epsg,int n,double x[]) {int i,j,m;
double *xx[4],*ss,*s;
double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
ss=(double *)malloc(n*(n+1)*sizeof(double));
s=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{for(j=0;j<=n;j++)
*(ss+i*(n+1)+j)=0;
*(ss+i*(n+1)+i)=1;
}
for(i=0;i<4;i++)
xx[i]=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
*(xx[0]+i)=p[i];
for(;;)
{for(i=0;i<n;i++)
{*(xx[1]+i)=*(xx[0]+i);
x[i]=*(xx[1]+i);
}
f0=f1=objf(x);
dlt=-1;
for(j=0;j<n;j++)
{for(i=0;i<n;i++)
{*(xx[0]+i)=x[i];
*(s+i)=*(ss+i*(n+1)+j);
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
df=f0-f;
if(df>dlt)
{dlt=df;
m=j;
}
}
sdx=0;
for(i=0;i<n;i++)
sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
if(sdx<eps)
{free(ss);
free(s);
for(i=0;i<4;i++)
free(xx[i]);
return(f);
}
for(i=0;i<n;i++)
*(xx[2]+i)=x[i];
f2=f;
for(i=0;i<n;i++)
{*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
x[i]=*(xx[3]+i);
}
fx=objf(x);
f3=fx;
q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
d=0.5*dlt*(f1-f3)*(f1-f3);
if((f3<f1)||(q<d))
{if(f2<=f3)
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[2]+i);
else
for(i=0;i<n;i++)
*(xx[0]+i)=*(xx[3]+i);
}
else
{for(i=0;i<n;i++)
{*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
*(s+i)=*(ss+(i+1)*(n+1));
}
f=oneoptim(xx[0],s,h0,epsg,n,x);
for(i=0;i<n;i++)
*(xx[0]+i)=x[i];
for(j=m+1;j<=n;j++)
for(i=0;i<n;i++)
*(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
}
}
}
void main()
{double p[]={1,2};
double ff,x[2];
ff=powell(p,0.3,0.001,0.0001,2,x);
printf("x[0]=%f,x[1]=%f,ff=%f\n",x[0],x[1],ff); getchar();
}
变尺度法
计算f(x1,x2)=x1^2+2*x2^2-4*x1-2*x1*x2 的无约束极值，初始点x0=[1,1]。

/*
tt ---- 一维搜索初始步长
ff ---- 差分法求梯度时的步长
ac ---- 终止迭代收敛精度
ad ---- 一维搜索收敛精度
n ----- 设计变量的维数
xk[n] -- 迭代初始点
*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<conio.h>
#define tt 0.01
#define ff 1.0e-6
#define ac 1.0e-6
#define ad 1.0e-6
#define n 2
double ia;
double fny(double *x)
{
double x1=x[0],x2=x[1];
double f;
f=x1*x1+2*x2*x2-4*x1-2*x1*x2;
return f;
}
double * iterate(double *x,double a,double *s)
{
double *x1;
int i;
x1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
x1[i]=x[i]+a*s[i];
return x1;
double func(double *x,double a,double *s)
{
double *x1;
double f;
x1=iterate(x,a,s);
f=fny(x1);
return f;
}
void finding(double a[3],double f[3],double *xk,double *s) {
double t=tt;
int i;
double a1,f1;
a[0]=0;f[0]=func(xk,a[0],s);
for(i=0;;i++)
{
a[1]=a[0]+t;
f[1]=func(xk,a[1],s);
if(f[1]<f[0]) break;
if(fabs(f[1]-f[0])>=ad)
{
t=-t;
a[0]=a[1];f[0]=f[1];
}
else
{
if(ia==1) return; //break
t=t/2;ia=1;
}
}
for(i=0;;i++)
{
a[2]=a[1]+t;
f[2]=func(xk,a[2],s);
if(f[2]>f[1]) break;
t=2*t;
a[0]=a[1];f[0]=f[1];
a[1]=a[2];f[1]=f[2];
}
if(a[0]>a[2])
a1=a[0];
f1=f[0];
a[0]=a[2];
f[0]=f[2];
a[2]=a1;
f[2]=f1;
}
return;
}
double lagrange(double *xk,double *ft,double *s)
{
int i;
double a[3],f[3];
double b,c,d,aa;
finding(a,f,xk,s);
for(i=0;;i++)
{
if(ia==1) { aa=a[1]; *ft=f[1]; break; }
d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]);
if(fabs(d)==0) break;
c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d;
if(fabs(c)==0) break;
b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]);
aa=-b/(2*c);
*ft=func(xk,aa,s);
if(fabs(aa-a[1])<=ad) {if(*ft>f[1]) aa=a[1];break;}
if(aa>a[1])
{
if(*ft>f[1]) {a[2]=aa;f[2]=*ft;}
else if(*ft<f[1]) {a[0]=a[1];a[1]=aa;f[0]=f[1];f[1]=*ft;}
else if(*ft==f[1])
{
a[2]=aa;a[0]=a[1];
f[2]=*ft;f[0]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
else
{
if(*ft>f[1]) {a[0]=aa;f[0]=*ft;}
else if(*ft<f[1]) {a[2]=a[1];a[1]=aa;f[2]=f[1];f[1]=*ft;} else if(*ft==f[1])
{a[0]=aa;a[2]=a[1];
f[0]=*ft;f[2]=f[1];
a[1]=(a[0]+a[2])/2;
f[1]=func(xk,a[1],s);
}
}
}
if(*ft>f[1]) {*ft=f[1];aa=a[1];}
return aa;
}
double *gradient(double *xk)
{
double *g,f1,f2,q;
int i;
g=(double*)malloc(n*sizeof(double));
f1=fny(xk);
for(i=0;i<n;i++)
{q=ff;
xk[i]=xk[i]+q; f2=fny(xk);
g[i]=(f2-f1)/q; xk[i]=xk[i]-q;
}
return g;
}
double * bfgs(double *xk)
{
double u[n],v[n],h[n][n],dx[n],dg[n],s[n];
double aa,ib;
double *ft,*xk1,*g1,*g2,*xx,*x0=xk;
double fi;
int i,j,k;
ft=(double *)malloc(sizeof(double));
xk1=(double *)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{
s[i]=0;
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
}
g1=gradient(xk);
fi=fny(xk);
x0=xk;
for(k=0;k<n;k++)
{
ib=0;
if(ia==1) { xx=xk; break; }
ib=0;
for(i=0;i<n;i++) s[i]=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
s[i]+= -h[i][j]*g1[j];
aa=lagrange(xk,ft,s);
xk1=iterate(xk,aa,s);
g2=gradient(xk1);
for(i=0;i<n;i++)
if((fabs(g2[i])>=ac)&&(fabs(g2[i]-g1[i])>=ac))
{ib=ib+1;}
if(ib==0) { xx=xk1; break; }
fi=*ft;
if(k==n-1)
{ int j;
xk=xk1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
h[i][j]=0;
if(j==i) h[i][j]=1;
}
g1=g2; k=-1;
}
else
{
int j;
double a1=0,a2=0;
for(i=0;i<n;i++)
{
dg[i]=g2[i]-g1[i];
dx[i]=xk1[i]-xk[i];
}
for(i=0;i<n;i++)
{
int j;
u[i]=0;v[i]=0;
for(j=0;j<n;j++)
{
u[i]=u[i]+dg[j]*h[j][i];
v[i]=v[i]+dg[j]*h[i][j];
}
}
for(j=0;j<n;j++)
{
a1+=dx[j]*dg[j];
a2+=v[j]*dg[j];
}
if(fabs(a1)!=0)
{
a2=1+a2/a1;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
h[i][j]+=(a2*dx[i]*dx[j]-v[i]*dx[j]-dx[i]*u[j])/a1;
}
xk=xk1; g1=g2;
}
}
if(*ft>fi) { *ft=fi; xx=xk;}
xk=x0;
return xx;
}
void main ()
{
int k;
double *xx,f;
double xk[n]={1,1};
xx=bfgs(xk);
f=fny(xx);
printf("\n\nThe Optimal Design Result Is:\n");
for(k=0;k<n;k++)
{printf("\n\tx[%d]*=%f",k+1,xx[k]);}
printf("\n\tf*=%f",f);
getch();
}
共轭梯度法的C+=语言实现
#include <iostream.h>
#include "Matrix.h"
#include<LIMITS.H>
#define MAX_M 2048
#include<math.h>
#define beta (sqrt(5)-1)/2
#define eclipse 0.01
CMatrix hesse()
{
CMatrix temp(2,2);
temp.A[0][0] = 2;
temp.A[0][1] = 0;
temp.A[1][0] = 0;
temp.A[1][1] = 8;
return temp;
}
double fun(double t,CMatrix &X1,CMatrix &X2)
{
double x1 = X1.A[0][0] + t*X2.A[0][0];
double x2 = X1.A[1][0] + t*X2.A[1][0];
return x1*x1 + 4*x2*x2;
}
double Max_Num(double a,double b)
{
if(a>b)
{
return a;
}
else
{
return b;
}
}
double Min_Num(double a,double b)
{
if(a<b)
{
return a;
}
else
{
return b;
}
}
void StepAdding_Search(double &a,double &b,CMatrix &mt1,CMatrix &mt2) {
double t[MAX_M]={0};
double h[MAX_M]={0};
double f[MAX_M]={0};
double result=0;
double p=2.1;
t[0]=0;
h[0]=1;
f[0]=fun(t[0],mt1,mt2);
for(int k=0;k<MAX_M-1;k++)
{
t[k+1]=t[k]+h[k];
f[k+1]=fun(t[k+1],mt1,mt2);
if(f[k+1]<f[k])
{
h[k+1]=p*h[k];
result=t[k];
t[k]=t[k+1];
f[k]=f[k+1];
}
else if(k == 0)
{
h[k]=-h[k];
result=t[k+1];
}
else
{
a=Min_Num(result,t[k+1]);
b=Max_Num(result,t[k+1]);
}
}
}
double Fibonacci(double &a,double &b,CMatrix mt1,CMatrix mt2) {
double t2,t1,result_1,result_2;
t2 = a + beta*(b-a);
result_2 = fun(t2,mt1,mt1);
while(true)
{
t1 = a+b-t2;
result_1 = fun(t1,mt1,mt2);
if(fabs(t1-t2)<eclipse)
{
return (t1+t2)/2;
}
else if(result_1 < result_2)
{
b = t1;
t1 = t2;
result_1 = result_2;
t2 = a+beta*(b-a);
result_2 = fun(t2,mt1,mt2);
}
else if(result_1 == result_2)
{
a = t1;
b = t2;
result_2 = fun(a,mt1,mt2);
}
else
{
a = t2;
t2 = t1;
result_2 = result_1;
}
}
}
double distance(CMatrix &X1,CMatrix &X2)
{
double x1 = X1.A[0][0] - X2.A[0][0];
double x2 = X1.A[1][0] - X2.A[1][0];
return x1*x1 + x2*x2;
}
CMatrix diff_fun(CMatrix &mt)
{
CMatrix temp(2,1);
temp.A[0][0] = 2*mt.A[0][0];
temp.A[1][0] = 8*mt.A[1][0];
return temp;
}
double fanshu(CMatrix mt)
{
return sqrt(mt.A[0][0]*mt.A[0][0] + mt.A[1][0]*mt.A[1][0]);
}
void main()
{
int i =0;
int n = 2000;
CMatrix temp_X1;
CMatrix X1(2,1);
X1.A[0][0] = 1;
X1.A[1][0] = 1;
CMatrix m_temp = hesse();
CMatrix n_temp = diff_fun(X1);
CMatrix X2 = (m_temp.Inverse()* n_temp).Matrix_Multiply(-1);
double a,b;
StepAdding_Search(a,b,X1,X2);
double m = Fibonacci(a,b,X1,X2);
temp_X1 = X1;
CMatrix X3 = X1 + X2.Matrix_Multiply(m);
while(distance(X1,X3) > eclipse)
{
X1 = X3;
if(i+1 == n)
{
i = 0;
X2 = m_temp.Inverse()* diff_fun(X1);
}
else
{
m = (fanshu(temp_X1)*fanshu(temp_X1))/(fanshu(X1)*fanshu(X1));
X2 = diff_fun(X1).Matrix_Multiply(-1) + X2.Matrix_Multiply(m);
}
X1 = X3;
X2 = m_temp.Inverse()* diff_fun(X1);
StepAdding_Search(a,b,X1,X2);
m = Fibonacci(a,b,X1,X2);
X3 = X1 + X2.Matrix_Multiply(m);
i++;
}
cout<<"求解的结果是:"<<endl;
cout<<"("<<X1.A[0][0]<<","<<X1.A[1][0]<<")"<<endl;
cout<<"最优解是:"<<endl;
cout<<fun(m,X1,X2)<<endl;
}
用最速下降法求极值。

#include <stdio.h>
#include <math.h>
#define EP 0.001
float H,s[2];
float x[2]={1,1};
float f(float *q)
{
float p;
p=pow(q[0],2)+2*pow(q[1],2)-4*q[0]-2*q[0]*q[1];
return(p);
}
float h()
{
float hh;
hh=(x[0]*s[0]+2*x[1]*s[1]-2*s[0]-s[1]*x[0]-s[0]*x[1])/(pow(s[0],2)+2*pow(s[1],2 )-2*s[0]*s[1]);
/* you wai bu suan shu zui jia bu chang.*/
return(hh);
}
main()
{
float y,R;
int i=0,n=0;
do
{
n=n+1;
s[0]=2*x[0]-4-2*x[1];
s[1]=4*x[1]-2*x[0];
R=sqrt(pow(s[0],2)+pow(s[1],2));
printf("\n R is %f.",R);
H=h();
for(i=0;i<=1;i++)
x[i]=x[i]-H*s[i];
}
while(R>EP);
y=f(x);
printf("\n The Min is %f.",y);
for(i=0;i<=1;i++)
printf("\n The X%d is %f.",i,x[i]); printf("\n The Number is %d.",n);
}。