机械优化设计黄金分割法 外推法

大学机械优化设计部分程序1.外推法2.黄金分割法3.二次插值法4.坐标轮换法5.随机方向法6.四杆机构优化设计1.外推法源程序：#include<stdio.h>#include<math.h>#define R 0.01double fun(double x){ double m;m=x*x-10*x+36;return m;}void main(){double h0=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<y2){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.4576002.黄金分割法源程序：#include<stdio.h>#include<math.h>#define f(x) x*x*x*x-5*x*x*x+4*x*x-6*x+60double hj(double *a,double *b,double e,int *n){double x1,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);}return s;}void main() {double s,a,b,e,m;int n=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值-350.0001a=3.279466,b=3.279793,s=22.659008,m =3.279629,n=213.二次插值法源程序：#include<stdio.h>#include<math.h>int main(void){double a1,a2,a3,ap,y1,y2,y3,yp,c1,c2,m; double j[3];int i,h=1;void finding(double a[3]);finding(j);a1=j[0];a2=j[1];a3=j[2];m=0.001;double f(double x);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)<m)break;else if((ap-a2)*h>0) {if(y2>=yp){a1=a2;y1=y2;a2=ap;y2=yp;}else{a3=ap;y3=yp;} }else if(y2>=yp){a3=a2;y3=y2;a2=ap;y2=yp;} else{a1=ap;y1=yp;} }double x,y;if(y2<=yp){x=a2;y=y2;}else{x=ap;y=yp;}printf("a*=%f\n",x); printf("y*=%f\n",y);return 0;}double f(double x){double y;y=x*x-10*x+36;return y;}void finding(double a[3]) {int h,i;double y[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]<y[1]);}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.000000y*=11.0000004.坐标轮换法源程序：#include <stdio.h>#include <math.h>#include <conio.h>float fun1(float x,float a,float b) {float y;y=x+a*b;return y;}float fun2(float x,float y){float z;z=4*(x-5)*(x-5)+(y-6)*(y-6); return z;}main(){floatd[100][3],x[100][3],xx[3],ax[100][3]; float a1,a2,a3,h,t,y1,y2,y3,e,a,b,l,fi; int i,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(y3<y2);for(;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)/y 2)>=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[i-1][1],d[i][1],ax[k][i]);x[i][2]=fun1(x[i-1][2],d[i][2],ax[k][i]);}}while(sqrt(pow((x[2][1]-x[0][1]),2)+pow( (x[2][2]-x[0][2]),2))>=1e-6);xx[1]=x[2][1];xx[2]=x[2][2];fi=fun2(xx[1],xx[2]);printf("最优解为\nx1*=%f\nx2*=%f\nf*=%f\nk=%d\n",xx[1],xx[2],fi,k);}运行过程及结果：输入初始点坐标89最优解为x1*=5.000000x2*=6.000000f*=0.000000k=25.随机方向法源程序：#include<math.h>#include<stdio.h>#include<stdlib.h> float f(float x,float y){float z;z=(x-2)*(x-2)+(y-1)*(y-1);return z;}float g1(float x,float y){float z;z=x*x-y;return z;}float g2(float x,float y){float z;z=x+y-2;return z;}void main(){int i,j;floatk=8,c=0.000001,a0=-3,b0=3,a1=-3,b1= 3;floatx[10],x0[10],xl[10],e[10],r[10],d[10],h,fl,f0 ,fx;while(g1(x0[0],x0[1])>0||g2(x0[0],x0[1])> 0){x0[0]=a0+(rand()/32767.00)*(b0-a0);x0[1]=a1+(rand()/32767.00)*(b1-a1);}fl=f(x0[0],x0[1]);f0=f(x0[0],x0[1]);while(1){h=0.01;j=1;r[0]=-1+(rand()/32767.00)*(1-(-1));r[1]=-1+(rand()/32767.00)*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]); e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0) {fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}while(j<=k){j++;r[0]=-1+(rand()/32767.00)*(1-(-1));r[1]=-1+(rand()/32767.00)*(1-(-1));e[0]=r[0]/sqrt(r[0]*r[0]+r[1]*r[1]);e[1]=r[1]/sqrt(r[0]*r[0]+r[1]*r[1]);x[0]=x0[0]+h*e[0];x[1]=x0[1]+h*e[1];if(g1(x[0],x[1])<=0&&g2(x[0],x[1])<=0){fx=f(x[0],x[1]);if(fx<fl){fl=fx;for(i=0;i<2;i++){d[i]=e[i];xl[i]=x[i];}}}}x[0]=xl[0];x[1]=xl[1];while(1){h=1.3*h;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)break;fx=f(x[0],x[1]);if(fx<fl) fl=fx;else break;}do{x[0]=x[0]-h*d[0];x[1]=x[1]-h*d[1];h=0.7*h;if(h<c)break;x[0]=x[0]+h*d[0];x[1]=x[1]+h*d[1];if(g1(x[0],x[1])>0||g2(x[0],x[1])>0)continue;fx=f(x[0],x[1]);}while(fx>=fl);if(fabs((f0-fx)/f0)>=c){x0[0]=x[0];x0[1]=x[1];fl=fx;f0=fx;}elsebreak;}printf("输出最优解为\nx1*=%f,x2*=%f, y*=%f\n",x[0],x[1],fx);}运行过程及结果：输出最优解为x1*=0.995421,x2*=1.004521,y*=1.0092 006.四杆机构优化设计源程序：#include<math.h>#include<stdio.h>#include<stdlib.h>#define Pai 3.1415926int g(float l1,float l2){if((-l1<=0)&&(-l2<=0)&&(6-l1-l2<=0)&&(1-l2-4<=0)&&(l2-l1-4<=0)&&(l1*l1+l2*l2-1.414*l1*l2-16<=0)&&(36-l1*l1-l2*l2-1.414*l1*l2<=0))return (1);elsereturn (0);}float fun(float x0[2]){floatf,a[31],b[31],r[31],p[31],q[31],w[31],x1[2];int i;p[0]=acos(((1+x0[0])*(1+x0[0])-x0[1] *x0[1]+25)/(10+10*x0[0]));q[0]=acos(((1+x0[0])*(1+x0[0])-x0[1] *x0[1]-25)/(10*x0[1]));f=0;for(i=1;i<=30;i++){p[i]=p[0]+(Pai/60)*i;r[i]=sqrt(26-10*cos(p[i]));a[i]=acos((r[i]*r[i]+x0[1]*x0[1]-x0[0]* x0[0])/(2*r[i]*x0[1]));b[i]=acos((r[i]*r[i]+24)/(10*r[i]));q[i]=Pai-a[i]-b[i];w[i]=q[0]+(2*(p[i]-p[0])*(p[i]-p[0]))/( 3*Pai);f=f+(Pai/60)*(q[i]-w[i])*(q[i]-w[i])*(p[ i]-p[i-1]);}return f;}void main(){floata,q,f,fl,f0,l[2],z[2],d0[100],d1[100],x[2],xi[ 2],fx,m0,m1,e;int i,j,n,k;printf("输入精度");scanf("%f",&e);do{z[0]=0+5*(rand()/32767.00);z[1]=0+5*(rand()/32767.00);}while(g(z[0],z[1])==0);for(i=0;i<=99;i++){d0[i]=-1+2*(rand()/32767.00);}for(j=0;j<=99;j++){d1[j]=-1+2*(rand()/32767.00);}f0=fun(z);fl=fun(z);ss:a=0.01;for(i=0,j=0;i<=99&&j<=99;i++,j++) {n=1/sqrt((d0[i])*(d0[i])+d1[j]*d1[j]);d0[i]=n*d0[i];d1[j]=n*d1[j];x[0]=z[0]+a*d0[i];x[1]=z[1]+a*d1[j];if(g(x[0],x[1])==1){f=fun(x);if(f<fl){fl=f;m0=d0[i];m1=d1[j];l[0]=x[0];l[1]=x[1];}}}x[0]=l[0];x[1]=l[1];do{a=1.3*a;x[0]=x[0]+a*m0;x[1]=x[1]+a*m1;if(g(z[0],z[1])==0) break;f=fun(x);if(f<fl)fl=f;else break;} while(g(z[0],z[1])==1); do{x[0]=x[0]-a*m0;x[1]=x[1]-a*m1;a=0.7*a;if(a<0.00001)break; x[0]=x[0]+a*m0;x[1]=x[1]+a*m1;if(g(z[0],z[1])==1)f=fun(x); }while(f>=fl);if(fabs((f0-f)/f0)<e){xi[0]=x[0];xi[1]=x[1];fx=f;printf("最优解为\nx1*=%f\nx2*=%f\nfx=%f\n",xi[0],xi[1],fx);}else{f0=f;fl=f;z[0]=x[0];z[1]=x[1];goto ss;}}运行过程及结果：输入精度0.001最优解为x1*=4.161386x2*=2.311257fx=0.000021。