baoweier.zip_鲍威尔_鲍威尔法_鲍威尔算法

#include <iostream> using namespace std; int ncom; double pcom[51], xicom[51]; const int N=2;//维数 ////////////////定义符号函数////////////////////////// int sgn(double x) { if (x>0) { return 1; } else { if(x<0) { return -1; } } return 0; } /////////////////定义多元函数////////////////////////// double f(double x[]) { return 4*pow(x[1],2)+2*pow(x[2],2)+9*x[1]-32*x[2]+16; } /////////////多元函数沿某指定方向点的函数值////////////// double func(double x) { double xt[51]; for (int j = 1; j<=ncom; j++) { xt[j] = pcom[j] + x * xicom[j]; } return f(xt); } ////////////////////////////////////////搜索区间/////////////////////////////// void find(double& ax, double& bx, double& cx, double& fa, double& fb, double& fc) { double r,q,dum,gold = 1.618034; int glimit = 100; double u,ulim,fu,tiny = 1e-20; fa = func(ax); fb = func(bx); if (fb > fa) { dum = ax; ax = bx; bx = dum; dum = fb; fb = fa; fa = dum; } cx = bx + gold * (bx - ax); fc = func(cx); while (fb >= fc) { r = (bx - ax) * (fb - fc); q = (bx - cx) * (fb - fa); dum = q - r; if (fabs(dum) < tiny) { dum = tiny; } u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum); ulim = bx + glimit * (cx - bx); if ((bx - u) * (u - cx) > 0) { fu = func(u); if (fu < fc) { ax = bx; fa = fb; bx = u; fb = fu; return; } else { if (fu > fb) { cx = u; fc = fu; return; } } u = cx + gold * (cx - bx); fu = func(u); } else { if ((cx - u) * (u - ulim) > 0) { fu = func(u); if (fu < fc) { bx = cx; cx = u; u = cx + gold * (cx - bx); fb = fc; fc = fu; fu = func(u); } } else { if ((u - ulim) * (ulim - cx) >= 0) { u = ulim; fu = func(u); } else { u = cx + gold * (cx - bx); fu = func(u); } } } ax = bx; bx = cx; cx = u; fa = fb; fb = fc; fc = fu; } } ////////////////////求极小点的一维搜索函数///////////////////////////// double brent(double ax, double bx, double cx, double tol, double& xmin) { int done,iter,itmax = 100; double d,fu,r,q,p,xm,tol1,tol2,a,b,cgold = 0.381966; double u,etemp,dum,v,w,x,e,fx,fv1,fw,zeps = 0.0000000001; a = ax; if (cx < ax) { a = cx; } b = ax; if (cx > ax) { b = cx; } v = bx; w = v; x = v; e = 0.0; fx = func(x); fv1 = fx; fw = fx; for (iter = 1; iter<=itmax; iter++) { xm = 0.5 * (a + b); tol1 = tol * fabs(x) + zeps; tol2 = 2.0 * tol1; if (fabs(x - xm) <= tol2 - 0.5 * (b - a)) { break; } done = -1; if (fabs(e) > tol1) { r = (x - w) * (fx - fv1); q = (x - v) * (fx - fw); p = (x - v) * q - (x - w) * r; q = 2.0 * (q - r); if (q > 0.0) { p = -p; } q = fabs(q); etemp = e; e = d; dum = fabs(0.5 * q * etemp); if (fabs(p) < dum && p > q * (a - x) && p < q * (b - x)) { d = p / q; u = x + d; if (u - a < tol2 || b - u < tol2) { d = fabs(tol1) * sgn(xm - x); } done = 0; } } if (done) { if (x >= xm) { e = a - x; } else { e = b - x; } d = cgold * e; } if (fabs(d) >= tol1) { u = x + d; } else { u = x + fabs(tol1) * sgn(d); } fu = func(u); if (fu <= fx) { if (u >= x) { a = x; } else { b = x; } v = w; fv1 = fw; w = x; fw = fx; x = u; fx = fu; } else { if (u < x) { a = u; } else { b = u; } if (fu <= fw || w == x) { v = w; fv1 = fw; w = u; fw = fu; } else { if (fu <= fv1 || v == x || v == w) { v = u; fv1 = fu; } } } } if (iter > itmax) { cout<<" brent exceed maximum iterations."<<endl; } xmin = x; return fx; } //////////////////求多元函数沿某指定方向的极小值点/////////////////////// void min(double p[], double xi[], int n, double& fret) { int j; double tol = 0.0001; ncom = n; for (j = 1; j<=n; j++) { pcom[j] = p[j]; xicom[j] = xi[j]; } double fa,fx,fb,bx,ax = 0.0; double xmin,xx = 1.0; find(ax, xx, bx, fa, fx, fb); fret = brent(ax, xx, bx, tol, xmin); for (j = 1; j<=n; j++) { xi[j] = xmin * xi[j]; p[j] = p[j] + xi[j]; } } /////////////////////置零////////////////////////////////// void erase(double pbar[], double prr[], double pr[]) { for (int i=0; i<=20; i++) { pbar[i]=0; prr[i]=0; pr[i]=0; } } ////////////////////////////////鲍威尔/////////////////////////////////////////////// void powell(double p[], double xi[], int n, double ftol, int& iter, double& fret) { int i,j,ibig,itmax = 200; double pt[21], ptt[21], xit[21]; double t,dum,fp,del,fptt; fret = f(p); for (j = 1; j<=n; j++) { pt[j] = p[j]; } iter = 0; do { do { do { iter = iter + 1; fp = fret; ibig = 0; del = 0.0; for (i = 1; i<=n; i++) { for (j = 1; j<=n; j++) { xit[j] = xi[(j-1)*N+i]; } fptt = fret; min(p, xit, n, fret); if (fabs(fptt - fret) > del) { del = fabs(fptt - fret); ibig = i; } } if (2.0 * fabs(fp - fret) <= ftol * (fabs(fp) + fabs(fret))) { erase(xit, ptt, pt); return; } if (iter == itmax) { cout<<" 超出范围"<<endl; return; } for (j = 1; j<=n; j++) { ptt[j] = 2.0 * p[j] - pt[j]; xit[j] = p[j] - pt[j]; pt[j] = p[j]; } fptt = f(ptt); }while (fptt >= fp); dum = fp - 2 * fret + fptt; t = 2.0 * dum * pow((fp - fret - del) , 2) - del * pow((fp - fptt) , 2); }while (t >= 0.0); min(p, xit, n, fret); for (j = 1; j<=n; j++) { xi[(j-1)*N+ibig] = xit[j]; } }while(1); } void main() { int i,j,iter,ndim; ndim = N;//维数 double fret,ftol = 0.000001;//精度 double p[4], xi[10]; for (i = 1; i<=N; i++) { for (j = 1; j<=N; j++) { xi[(i-1)*N+j] = 0.0; } } xi[1] = 1.0; xi[1*N+2] = 1.0; xi[2*N+3] = 1.0; p[1] = 1.5; p[2] = 1.5; p[3] = 2.5; powell(p, xi, ndim, ftol, iter, fret); cout<<"最小值点为: "<<endl; cout.setf(ios::fixed|ios::left); cout.precision(6); for (i = 1; i<=N; i++) { cout.width(12); cout<<p[i]; } cout<<endl; cout<<"最小函数值为 = "<<fret; cout<<endl; }