Powell法求多元函数极小值（要求目标函数有连续二阶导数）.docVIP

Powell法求多元函数极小值（要求目标函数有连续二阶导数） #include conio.h#include math.h#include stdio.h#include stdlib.h#define XuYimin_Debug_1#define N 3??/*N维空间*/double xq[N]={0,0,0};double fun(double x[N]){??double y;/*目标函数f(X)=(x1-4)^2+(x2-3)^2+1.5*x3^2+(x1-x2+x3-1)^2+sin(x1+x2+x3-7)*sin(x1+x2+x3-7)可根据需要修改，但需要同时修改维数N值*/??y=pow((x[0]-4),2)+pow((x[1]-3),2)+1.5*x[2]*x[2]??? +pow(x[0]-x[1]+x[2]-1,2)+pow(sin(x[0]+x[1]+x[2]-7),2);??return y;}double fd(double x[N],double d[N],double alf){??int i;??double y,t[N];??for(i=0;iN;i++)??? t[i]=x[i]+alf*d[i];??y=fun(t);??return y;}double alf(double x[N],double d[N]){??double h0=0.1,a0=0.0,ep1;??double a1,a2,a3,y1,y2,y3,c1,c2,ap,yp,ax,yx,h;??ep1=pow(0.1,12);??a1=a0,h=h0;a2=h,y1=fd(x,d,a1),y2=fd(x,d,a2);??if(y2y1)??? {??????h=-1*h,a3=a1,y3=y1;??????loop1:a1=a2,y1=y2,a2=a3,y2=y3;??????loop2:a3=a2+h,y3=fd(x,d,a3); if(y3y2) {h=2*h;goto loop1;}??? }??else??goto loop2;??if( fabs(y2)ep1 )??goto loop5;??loop3:??? c1=(y3-y1)/(a3-a1);??? c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);??? if(c2==0) goto loop4;??? ap=0.5*(a1+a3-c1/c2);??? yp=fd(x,d,ap);??while(!fabs((y2-yp)/y2)ep1)??? {??????if((ap-a2)*h0)????????{????????? 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;????????}??????goto loop3;??? }??if(y2yp)??ax=a2,yx=y2;??else?????? ax=ap,yx=yp;??return ax;??loop4:return a2;??loop5:return 0;}double maxt(double del[N+1]){??int i;??double y=del[0];??for (i=1;iN+1;i++)??? if(del[i]y)??????y=del[i];??return y;}void powell(){??int i,j,id,pd[3],xym=1,m;??double al,esp,ALF[N+1],x[2],pa[N+1],X[N][N+1],del[N+1],d[N][N+1],??????F[3],fu[N+1],te[N],tem[N],dk[N],xk[N],delm;??esp=pow(0.1,10);x[0]=xq[0];x[1]=xq[1];??for(i=0;iN;i++)??? for(j=0;jN+1;j++)??????d[i][j]=0;??for(i=0;iN;i++)??d[i][i+1]=1;??loop:??#ifdef XuYimin_Debug_1??? printf(fangxiang juzhen:d=\n);??? for(