机械优化设计黄金分割法外推法.docx
《机械优化设计黄金分割法外推法.docx》由会员分享，可在线阅读，更多相关《机械优化设计黄金分割法外推法.docx（20页珍藏版）》请在冰豆网上搜索。
机械优化设计黄金分割法外推法
郑州大学
机械优化设计部分程序
1.外推法
2.黄金分割法
3.二次插值法
4.坐标轮换法
5.随机方向法
6.四杆机构优化设计
1.外推法
源程序：
#include
#include
#defineR0.01
doublefun（doublex）
{doublem;
m=x*x-10*x+36;
returnm;
}
voidmain（）
{
doubleh0=R,y1,y2,y3,x1,x2,x3,h;
x1=0;h=h0;x2=h;
y1=fun（x1）;y2=fun（x2）;
if（y2>y1）
{h=-h;
x3=x1;
y3=y1;
x1=x2;
y1=y2;
x2=x3;
y2=y3;
}
x3=x2+h;y3=fun（x3）;
while（y3{h*=2.0;
x1=x2;
y1=y2;
x2=x3;
y2=y3;
x3=x2+h;
y3=fun（x3）;
}
printf（"fun（%f）=%f,fun（%f）=%f,fun（%f）=%f\n",x1,y1,x2,y2,x3,y3）;
}
运行过程及结果：
fun（2.560000）=16.953600,
fun（5.120000）=11.014400,
fun（10.240000）=38.457600
2.黄金分割法
源程序：
#include
#include
#definef（x）x*x*x*x-5*x*x*x+4*x*x-6*x+60
doublehj（double*a,double*b,doublee,int*n）
{
doublex1,x2,s;
if（fabs（（*b-*a）/（*b））<=e）
s=f（（*b+*a）/2）;
else
{
x1=*b-0.618*（*b-*a）;
x2=*a+0.618*（*b-*a）;
if（f（x1）>f（x2））
*a=x1;
else
*b=x2;
*n=*n+1;
s=hj（a,b,e,n）;
}
returns;
}
voidmain（）
{
doubles,a,b,e,m;
intn=0;
printf（"输入a，b值和精度e值\n"）;
scanf（"%lf%lf%lf",&a,&b,&e）;
s=hj（&a,&b,e,&n）;
m=（a+b）/2;
printf（"a=%lf,b=%lf,s=%lf,m=%lf,n=%d\n",a,b,s,m,n）;
}
运行过程及结果：
输入a，b值和精度e值
-3
5
0.0001
a=3.279466,b=3.279793,s=22.659008,m=3.279629,n=21
3.二次插值法
源程序：
#include
#include
intmain（void）
{
doublea1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m;
doublej[3];
inti,h=1;
voidfinding（doublea[3]）;
finding（j）;
a1=j[0];
a2=j[1];
a3=j[2];
m=0.001;
doublef（doublex）;
y1=f（a1）;
y2=f（a2）;
y3=f（a3）;
for（i=1;1>=1;i++）
{
c1=（y3-y1）/（a3-a1）;
c2=（（y2-y1）/（a2-a1）-c1）/（a2-a3）;
ap=0.5*（a1+a3-c1/c2）;
yp=f（ap）;
if（fabs（（y2-yp）/y2）break;
elseif（（ap-a2）*h>0）
{
if（y2>=yp）{
a1=a2;y1=y2;
a2=ap;y2=yp;}
else{
a3=ap;y3=yp;}
}
elseif（y2>=yp）{
a3=a2;y3=y2;
a2=ap;y2=yp;}
else{a1=ap;y1=yp;}
}
doublex,y;
if（y2<=yp）{
x=a2;y=y2;}
else{
x=ap;y=yp;}
printf（"a*=%f\n",x）;
printf（"y*=%f\n",y）;
return0;
}
doublef（doublex）
{
doubley;
y=x*x-10*x+36;
returny;
}
voidfinding（doublea[3]）
{
inth,i;
doubley[3];
a[0]=0;
h=1;
a[1]=h;
y[0]=f（a[0]）;y[1]=f（a[1]）;
if（y[1]>y[0]）
{
h=-h;
a[2]=a[0];y[2]=y[0];
do{
a[0]=a[1];a[1]=a[2];
y[0]=y[1];y[1]=y[2];
a[2]=a[1]+h;y[2]=f（a[2]）;
h=2*h;
}while（y[2]}
else{
for（i=1;i>=1;i++）{
a[2]=a[1]+h;y[2]=f（a[2]）;
if（y[2]>=y[1]）
break;
h=2*h;
a[0]=a[1];y[0]=y[1];
a[1]=a[2];y[1]=y[2];}
}
return;
}
运行过程及结果：
a*=5.000000
y*=11.000000
4.坐标轮换法
源程序：
#include
#include
#include
floatfun1（floatx,floata,floatb）
{floaty;
y=x+a*b;
returny;
}
floatfun2（floatx,floaty）
{floatz;
z=4*（x-5）*（x-5）+（y-6）*（y-6）;
returnz;
}
main（）
{floatd[100][3],x[100][3],xx[3],ax[100][3];
floata1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi;
inti,k;
printf（"输入初始点坐标\n"）;
scanf（"%f%f",&x[0][1],&x[0][2]）;
e=0.000001;
l=0.618;
x[2][1]=x[0][1];
x[2][2]=x[0][2];
k=0;
k--;
do
{x[0][1]=x[2][1];
x[0][2]=x[2][2];
k++;
for（i=1;i<=2;i++）
{
if（i==1）
{d[i][1]=1;
d[i][2]=0;
}
else
{d[i][1]=0;
d[i][2]=1;
}
h=0.1;
a1=0;
a2=h;
x[i][1]=fun1（x[i-1][1],d[i][1],a1）;
x[i][2]=fun1（x[i-1][2],d[i][2],a1）;
y1=fun2（x[i][1],x[i][2]）;
x[i][1]=fun1（x[i-1][1],d[i][1],a2）;
x[i][2]=fun1（x[i-1][2],d[i][2],a2）;
y2=fun2（x[i][1],x[i][2]）;
if（y2>y1）
{h=-h;
a3=a1;
y3=y1;
a1=a2;
a2=a3;
y1=y2;
y2=y3;
}
a3=a2+h;
x[i][1]=fun1（x[i-1][1],d[i][1],a3）;
x[i][2]=fun1（x[i-1][2],d[i][2],a3）;
y3=fun2（x[i][1],x[i][2]）;
do
{a1=a2;
y1=y2;
a2=a3;
y2=y3;
a3=a2+h;
x[i][1]=fun1（x[i-1][1],d[i][1],a3）;
x[i][2]=fun1（x[i-1][2],d[i][2],a3）;
y3=fun2（x[i][1],x[i][2]）;
}
while（y3for（;a1>a3;）
{t=a3;
a3=a1;
a1=t;
t=y1;
y3=y1;
y1=t;
}
a=a1;
b=a3;
a1=b-l*（b-a）;
a2=a+l*（b-a）;
x[i][1]=fun1（x[i-1][1],d[i][1],a1）;
x[i][2]=fun1（x[i-1][2],d[i][2],a1）;
y1=fun2（x[i][1],x[i][2]）;
x[i][1]=fun1（x[i-1][1],d[i][1],a2）;
x[i][2]=fun1（x[i-1][2],d[i][2],a2）;
y2=fun2（x[i][1],x[i][2]）;
if（b<1e-3）
{
for（;fabs（b-a）>e;）
{
if（y1>=y2）
{a=a1;
a1=a2;
y1=y2;
a2=a+l*（b-a）;
x[i][1]=fun1（x[i-1][1],d[i][1],a2）;
x[i][2]=fun1（x[i-1][2],d[i][2],a2）;
y2=fun2（x[i][1],x[i][2]）;
}
else
{b=a2;
a2=a1;
y2=y1;
a1=b-l*（b-a）;
x[i][1]=fun1（x[i-1][1],d[i][1],a1）;
x[i][2]=fun1（x[i-1][2],d[i][2],a1）;
y1=fun2（x[i][1],x[i][2]）;
}
}
}
else
{
for（;fabs（（b-a）/b）>=e||fabs（（y2-y1）/y2）>=e;）
{
if（y1>=y2）
{a=a1;
a1=a2;
y1=y2;
a2=a+l*（b-a）;
x[i][1]=fun1（x[i-1][1],d[i][1],a2）;
x[i][2]=fun1（x[i-1][2],d[i][2],a2）;
y2=fun2（x[i][1],x[i][2]）;
}
else
{b=a2;
a2=a1;
y2=y1;
a1=b-l*（b-a）;
x[i][1]=fun1（x[i-1][1],d[i][1],a1）;
x[i][2]=fun1（x[i-1][2],d[i][2],a1）;
y1=fun2（x[i][1],x[i][2]）;
}
}
}
ax[k][i]=0.5*（a+b）;
x[i][1]=fun1（x[
升级会员 