开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
Extremum.inl
资源名称：C++常用算法之06矩阵特征值与特征向量的计算.rar　[点击查看]
上传用户：fxromeo
上传日期：2010-04-08
资源大小：89k
文件大小：20k
源码类别：

数值算法/人工智能

开发平台：

Visual C++

Extremum.inl：源码内容
  1. // Extremum.inl 计算极值函数(方法)定义头文件
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 何渝, 2002
  4. // 最后修改: 2002.5.31.
  6. #ifndef _EXTREMUM_INL //避免多次编译
  7. #define _EXTREMUM_INL
  9. //一维极值连分式法
  10. template <class _Ty>
  11. void ExtremumFraction1D(valarray<_Ty>& x, _Ty eps, 
  12. int k, valarray<int>& js)
  13. {
  14.     _Ty xx,h1,dx,y[10],b[10],h2(0);
  15. valarray<_Ty> z(2);
  16.     int jt = 1;
  17. js[0] = 0;
  18.     while(jt == 1)
  19.     {
  20. int j=0;
  21.         while(j<=7)
  22.         {
  23. if(j <=2) 
  24. xx=x[0]+0.01*j;
  25.             else xx=h2;
  26.             FunctionValueEF1D(xx,z);
  27.             if(Abs(z[1]) < eps)
  28.             {
  29. jt=0; 
  30. j=10;
  31. }
  32.             else
  33.             {
  34. h1=z[1]; 
  35. h2=xx;
  36.                 if(j==0)
  37. {
  38. y[0]=h1;
  39. b[0]=h2;
  40. }
  41.                 else
  42.                 {
  43. y[j]=h1;
  44. int m=0;
  45. int i=0;
  46.                     while((m==0)&&(i<=j-1))
  47. {
  48. if(FloatEqual((h2-b[i]),0)) 
  49. m=1;
  50. else
  51. h2=(h1-y[i])/(h2-b[i]);
  52.                         i=i+1;
  53. }
  54.                     b[j]=h2;
  55.                     if(m!=0)
  56. b[j]=1.0e+35;
  57.                     h2=0.0;
  58.                     for(i=j-1; i>=0; i--)
  59. h2=-y[i]/(b[i+1]+h2);
  60.                     h2=h2+b[0];
  61.                 }
  62.                 j=j+1;
  63.             }
  64.         }
  65.         x[0]=h2;
  66.         js[0]=js[0]+1;
  67.         if(js[0]==k)  jt=0;
  68.     }
  69.     xx = x[0];
  70.     FunctionValueEF1D(xx,z);
  71. x[1] = z[0];
  72.     if(Abs(x[0]) < 1.0 || FloatEqual(Abs(x[0]), 1)) dx=1.0e-05;
  73.     else dx = Abs(x[0] * 1.0e-05);
  74.     xx = x[0] - dx;
  75.     FunctionValueEF1D(xx,z);
  76. h1 = z[0];
  77.     xx = x[0] + dx;
  78.     FunctionValueEF1D(xx,z);
  79. h2 = z[0];
  80.     js[1] = -1;
  81.     if((h1 + h2 - 2.0 * x[1]) < 0.0 || FloatEqual((h1 + h2 - 2.0 * x[1]),0)) js[1]=1;
  82. }
  84. //n维极值连分式法
  85. template <class _Ty>
  86. void ExtremumFractionND(valarray<_Ty>& x, _Ty eps, 
  87. int k, valarray<int>& js)
  88. {
  89.     int i,j,m,l,jt(1),il;
  90.     _Ty y[10],b[10],p,z,t,h1,h2(0),f,dx;
  92. int n = x.size()-1; //自变量个数
  94.     js[0] = 0;
  95. //jt = 1;
  96. //h2 = 0;
  97.     while(jt == 1)
  98.     {
  99. t=0.0; 
  100. js[0]=js[0]+1;
  101.         for(i=1; i<=n; i++)
  102.         {
  103. f=FunctionValueEFND(x,i);
  104.             t=t+Abs(f);
  105.         }
  106.         if(t<eps) jt=0;
  107.         else
  108.         {
  109. for(i = 0; i < n; i++)
  110.               { 
  111. il=5;
  112.                 while(il!=0)
  113.                 {
  114. j=0;
  115. t=x[i];
  116. il=il-1;
  117.                     while(j<=7)
  118.                     {
  119. if(j<=2) z=t+j*0.01;
  120.                         else z=h2;
  121.                         x[i]=z;
  122.                         f=FunctionValueEFND(x,i+1);
  123.                         if(FloatEqual(f,0))
  124.                         {
  125. j=10;
  126. il=0;
  127. }
  128.                         else
  129.                         {
  130. h1=f;
  131. h2=z;
  132.                             if(j==0)
  133.                             {
  134. y[0]=h1;
  135. b[0]=h2;
  136. }
  137.                             else
  138.                             {
  139. y[j]=h1;
  140. m=0;
  141. l=0;
  142.                                 while((m==0)&&(l<=j-1))
  143.                                 {
  144. p=h2-b[l];
  145.                                     if(FloatEqual(p,0)) m=1;
  146.                                     else h2=(h1-y[l])/p;
  147.                                     l=l+1;
  148.                                 }
  149.                                 b[j]=h2;
  150.                                 if(m!=0) b[j]=1.0e+35;
  151.                                 h2=0.0;
  152.                                 for(l=j-1; l>=0; l--)
  153.                                 h2=-y[l]/(b[l+1]+h2);
  154.                                 h2=h2+b[0];
  155.                             }
  156.                             j=j+1;
  157.                         }
  158.                     }
  159.                     x[i]=h2; 
  160.                 }
  161.                 x[i]=z;
  162.             }
  163.             if(js[0]==k) jt=0;
  164.         }
  165.     }
  166.     js[1]=1;
  167.     f=FunctionValueEFND(x,0); 
  168. x[n]=f;
  169.     dx=0.00001;
  170. t=x[0];
  171.     x[0]=t+dx;
  172. h1=FunctionValueEFND(x,0);
  173.     x[0]=t-dx;
  174. h2=FunctionValueEFND(x,0);
  175.     x[0]=t;
  176.     t=h1+h2-2.0*f;
  177.     if(t>0.0) js[1]=-1;
  178.     j=1; jt=1;
  179.     while(jt==1)
  180.     {
  181. j=j+1;
  182. dx=0.00001;
  183. jt=0;
  184.         t=x[j-1];
  185.         x[j-1]=t+dx;
  186. h2=FunctionValueEFND(x,0);
  187.         x[j-1]=t-dx;
  188. h1=FunctionValueEFND(x,0);
  189.         x[j-1]=t;
  190. t=h1+h2-2.0*f;
  191.         if((t*js[1]<0.0)&&(j<n)) jt=1;
  192.     }
  193.     if(t*js[1]>0.0) js[1]=0;
  194. }
  196. //线性规划
  197. template <class _Ty>
  198. int ExtremumLinePrograming(matrix<_Ty>& a, valarray<_Ty>& b, 
  199. valarray<_Ty>& c, valarray<_Ty>& x)
  200. {
  201.     _Ty y;
  203. int m = b.size(); //不等式约束条件的个数
  204. int n = c.size() - m; //变量个数
  206. valarray<int> js(m);
  207. matrix<_Ty> p(m, m);
  208. matrix<_Ty> d(m, (m+n));
  209. for(int i = 0; i < m; i++) js[i] = n + i;
  210.     int mn = m + n;
  211. _Ty s(0);
  212.     while(1==1)
  213.     {
  214. for(i = 0; i < m; i++)
  215.           for(int j = 0; j < m; j++)
  216.   p(i,j) = a(i, js[j]);
  217.         int l=MatrixInversionGS(p);
  218.         if(l==0)
  219.         {
  220. x[n]=s;
  221.             return(l);
  222.         }
  224. d = p * a; //矩阵相乘
  225.         
  226.         for(i=0; i<mn; i++) x[i]=0.0;
  227.         for(i=0; i<m; i++)
  228.         {
  229. s=0.0;
  230.             for(int j=0; j<m; j++)
  231.               s=s + p(i,j) * b[j];
  232.             x[js[i]]=s;
  233.         }
  234.         int k = -1;
  235. _Ty dd=1.0e-35;
  236.         for(int j = 0; j < mn; j++)
  237.         {
  238. _Ty z(0);
  239.             for(i=0; i<m; i++)
  240.               z = z + c[js[i]] * d(i, j);
  241.             z = z - c[j];
  242.             if(z > dd)
  243. {
  244. dd = z;
  245. k = j;
  246. }
  247.         }
  248.         if(k == -1)
  249.           {
  250. s=0.0;
  251.             for(j=0; j<n; j++)
  252. s = s + c[j] * x[j];
  253. x[n] = s;
  254.             return(1);
  255.         }
  256.         j = -1;
  257.         dd = 1.0e+20;
  258.         for(i=0; i<m; i++)
  259. if(d(i, k) >= 1.0e-20 )
  260.             { 
  261.   y=x[js[i]] / d(i, k);
  262.               if(y < dd)
  263.   {
  264.   dd = y; 
  265.   j = i;
  266.   }
  267. }
  268.         if(j==-1) 
  269.         {
  270. x[n]=s;
  271.             return(0);
  272.         }
  273.         js[j]=k;
  274.     }
  275. }
  277. //n维极值单形调优法
  278. template <class _Ty>
  279. int ExtremumSimplexND(_Ty d, _Ty u, _Ty v, valarray<_Ty>& x,  
  280. _Ty eps, int k, matrix<_Ty>& xx, valarray<_Ty>& f)
  281. {
  282. int r,g,l;
  283. _Ty fe,ft,ff;
  285. int n = x.size() - 1; //变量个数
  287. valarray<_Ty> xt(n);
  288. valarray<_Ty> xf(n);
  289. valarray<_Ty> xe(n);
  291.     int kk=0;
  292. _Ty nn=1.0*n;
  293.     _Ty fr=sqrt(nn+1.0);
  294.     _Ty fl=d*(fr-1.0)/(1.414*nn);
  295.     _Ty fg=d*(fr+nn-1.0)/(1.414*nn);
  296.     for(int i=0; i<n; i++)
  297.     for(int j=0; j<=n; j++)
  298.        xx(i,j) = 0;
  299.     for(i=1; i<=n; i++)
  300.     for(int j=0; j<n; j++)
  301.        xx(j,i) = fl;
  302.     for(i=1; i<=n; i++)
  303.        xx((i-1),i) = fg;
  304.     for(i=0; i<=n; i++)
  305.     { 
  306. for(int j=0; j<n; j++) xt[j]=xx(j,i);
  307.         f[i]=FunctionValueESND(xt,n);
  308.     }
  309.     ft=1.0+eps;
  310.     while((kk<k)&&(ft>eps))
  311.     {
  312. kk=kk+1;
  313.         fr=f[0]; fl=f[0]; r=0; l=0;
  314.         for(i=1; i<=n; i++)
  315.         {
  316. if(f[i]>fr)
  317. {
  318. r=i;
  319. fr=f[i];
  320. }
  321.             if(f[i]<fl)
  322. {
  323. l=i;
  324. fl=f[i];
  325. }
  326.         }
  327.         g=0;
  328. fg=f[0];
  329.         int j=0;
  330.         if(r==0)
  331. {
  332. g=1;
  333. fg=f[1];
  334. j=1;
  335. }
  336.         for(i=j+1; i<=n; i++)
  337. if((i!=r)&&(f[i]>fg))
  338.             {
  339. g=i;
  340. fg=f[i];
  341. }
  342.         for(j=0; j<n; j++)
  343.         {
  344. xf[j]=0.0;
  345.             for(i=0; i<=n; i++)
  346. if(i!=r) xf[j]=xf[j]+xx(j,i)/nn;
  347.             xt[j]=2.0*xf[j]-xx(j,r);
  348.         }
  349.         ft=FunctionValueESND(xt,n);
  350.         if(ft<f[l])
  351.         {
  352. for(j=0; j<n; j++) xf[j]=(1.0+u)*xt[j]-u*xf[j];
  353.             ff=FunctionValueESND(xf,n);
  354.             if(ff<f[l])
  355.             {
  356. for(j=0; j<n; j++) xx(j,r)=xf[j];
  357.                 f[r]=ff;
  358.             }
  359.             else
  360. {
  361. for(j=0; j<n; j++) xx(j,r)=xt[j];
  362.                 f[r]=ft;
  363.             }
  364.         }
  365.         else
  366.         {
  367. if(ft<=f[g])
  368.             {
  369. for(j=0; j<=n-1; j++) xx(j,r)=xt[j];
  370.                 f[r]=ft;
  371.             }
  372.             else 
  373.             {
  374. if(ft<=f[r])
  375.                 {
  376. for(j=0; j<=n-1; j++) xx(j,r)=xt[j];
  377.                     f[r]=ft;
  378.                 }
  379.                 for(j=0; j<n; j++) xf[j]=v*xx(j,r)+(1.0-v)*xf[j];
  380.                 ff=FunctionValueESND(xf,n);
  381.                 if(ff>f[r])
  382.                   for(i=0; i<=n; i++)
  383.                   {
  384.   for(j=0; j<=n-1; j++)
  385.   {
  386.   xx(j,i)=(xx(j,i)+xx(j,l))/2.0;
  387.                           x[j]=xx(j,i);
  388.   xe[j]=x[j];
  389.                       }
  390.                       fe=FunctionValueESND(xe,n);
  391.   f[i]=fe;
  392.                   }
  393.                 else
  394.                 {
  395. for(j=0; j<=n-1; j++) xx(j,r)=xf[j];
  396.                     f[r]=ff;
  397.                 }
  398.             }
  399.         }
  400.         ff=0.0;
  401. ft=0.0;
  402.         for(i=0; i<=n; i++)
  403.         {
  404. ff=ff+f[i]/(1.0+nn);
  405.             ft=ft+f[i]*f[i];
  406.         }
  407.         ft=(ft-(1.0+n)*ff*ff)/nn;
  408.     }
  409.     for(int j=0; j<n; j++)
  410.     {
  411. x[j] = 0.0;
  412.         for(i=0; i<=n; i++) x[j] = x[j] + xx(j,i) / (1.0+nn);
  413.         xe[j] = x[j];
  414.     }
  415.     fe = FunctionValueESND(xe,n);
  416. x[n] = fe;
  417.     return(kk);
  418. }
  420. //n维极值复形调优法
  421. template <class _Ty>
  422. int ExtremumComplexND(int m, valarray<_Ty>& a, valarray<_Ty>& b,
  423. _Ty alpha, _Ty eps, valarray<_Ty>& x, matrix<_Ty>& xx, int k)
  424. {
  425.     int r, g, kt, jt, kk;
  426.     _Ty fj, fr, fg, z;
  427. double rr = 0; //随机数种子
  428. int n = a.size(); //变量个数
  430. valarray<_Ty> c(m);
  431. valarray<_Ty> d(m);
  432. valarray<_Ty> w(m);
  433. valarray<_Ty> xt(m);
  434. valarray<_Ty> xf(m);
  436.     for(int i=0; i<n; i++) xx(i,0) = x[i];
  437.     xx(n,0) = FunctionValueTarget(x,n);
  438.     for(int j=1; j<2*n; j++)
  439.     {
  440. for(i=0; i<n; i++)
  441. {
  442. xx(i,j) = a[i] + (b[i] - a[i]) * rand_01_One(rr);
  443.             x[i] = xx(i,j);
  444.         }
  445.         int it = 1;
  446.         while(it==1)
  447.         {
  448. it = 0;
  449. r = 0;
  450. g = 0;
  451.             while((r<n)&&(g==0))
  452.             {
  453. if((a[r]<=x[r])&&(b[r]>=x[r])) r = r + 1;
  454.                 else g = 1;
  455.             }
  456.             if(g==0)
  457.             {
  458. ConditionValue(n,m,x,c,d,w);
  459.                 r = 0;
  460.                 while((r<m)&&(g==0))
  461.                 {
  462. if((c[r]<=w[r])&&(d[r]>=w[r])) r=r+1;
  463.                     else g=1;
  464.                 }
  465.             }
  466.             if(g!=0)
  467.             {
  468. for(r=0; r<n; r++)
  469.                 {
  470. z = 0.0;
  471.                     for(g=0; g<j; g++) z = z + xx(r,g) / (1.0*j);
  472.                     xx(r,j) = (xx(r,j) + z) / 2.0;
  473.                     x[r]=xx(r,j);
  474.                 }
  475.                 it = 1;
  476.             }
  477.             else xx(n,j) = FunctionValueTarget(x,n);
  478.         }
  479.     }
  480.     kk = 1;
  481. int it = 1;
  482.     while(it==1)
  483.     {
  484. it = 0;
  485.         fr = xx(n,0);
  486. r = 0;
  487.         for(i=1; i<2*n; i++)
  488.           if(xx(n,i)>fr)
  489.           {
  490.   r = i;
  491.   fr = xx(n,i);
  492.   }
  493.         g = 0;
  494. j = 0;
  495. fg = xx(n,0);
  496.         if(r==0)
  497.         {
  498. g = 1;
  499. j = 1;
  500. fg = xx(n,1);
  501. }
  502.         for(i=j+1; i<2*n; i++)
  503.           if(i!=r)
  504.             if(xx(n,i)>fg)
  505.             {
  506. g = i;
  507. fg = xx(n,i);
  508. }
  509.         for(i=0; i<n; i++)
  510.         {
  511. xf[i] = 0.0;
  512.             for(j=0; j<2*n; j++)
  513.               if(j!=r) xf[i]=xf[i]+xx(i,j)/(2.0*n-1.0);
  514.             xt[i]=(1.0+alpha)*xf[i]-alpha*xx(i,r);
  515.         }
  516.         jt = 1;
  517.         while(jt==1)
  518.         {
  519. jt = 0;
  520.             z = FunctionValueTarget(xt,n);
  521.             while(z>fg)
  522.             {
  523. for(i=0; i<n; i++) xt[i]=(xt[i]+xf[i])/2.0;
  524.                 z = FunctionValueTarget(xt,n);
  525.             }
  526.             j=0;
  527.             for(i=0; i<n; i++)
  528.             { 
  529. if(a[i]>xt[i])
  530.                 {
  531. xt[i]=xt[i]+0.000001;
  532. j=1;
  533. }
  534. if(b[i]<xt[i])
  535.                 {
  536. xt[i]=xt[i]-0.000001;
  537. j=1;
  538. }
  539.             }
  540.             if(j!=0) jt=1;
  541.             else
  542.             {
  543. ConditionValue(n,m,xt,c,d,w);
  544.                 j=0;
  545. kt=1;
  546.                 while((kt==1)&&(j<m))
  547.                 {
  548. if((c[j]<=w[j])&&(d[j]>=w[j])) j=j+1;
  549.                     else kt=0;
  550.                 }
  551.                 if(j<m)
  552.                 { 
  553. for(i=0; i<n; i++) xt[i]=(xt[i]+xf[i])/2.0;
  554.                     jt=1;
  555.                 }
  556.             }
  557.         }
  558.         for(i=0; i<n; i++) xx(i,r)=xt[i];
  559.         xx(n,r) = z;
  560.         fr = 0.0; 
  561. fg = 0.0;
  562.         for(j=0; j<2*n; j++)
  563.         {
  564. fj = xx(n,j);
  565.             fr = fr + fj / (2.0*n);
  566.             fg = fg + fj * fj;
  567.         }
  568.         fr = (fg - 2.0 * n * fr * fr) / (2.0 * n - 1.0);
  569.         if(fr>=eps)
  570.         {
  571. kk = kk + 1;
  572.             if(kk<k) it = 1;
  573.         }
  574.     }
  575.     for(i=0; i<n; i++)
  576.     {
  577. x[i] = 0.0;
  578.         for(j=0; j<2*n; j++) x[i]=x[i]+xx(i,j)/(2.0*n);
  579.     }
  580.     z = FunctionValueTarget(x,n);
  581. x[n] = z;
  582.     return(kk);
  583. }
  585. //确定1维函数极小值点所在区间
  586. template <class _Ty>
  587. void MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, _Ty& fa, _Ty& fb, _Ty& fc)
  588. {
  589.     _Ty r,q,dum,gold = 1.0 + GoldNo;
  590.     _Ty u,ulim,fu,tiny = 1e-20;
  591. int glimit = 100;
  593.     fa = FunctionValue(ax);
  594.     fb = FunctionValue(bx);
  595.     if (fb > fa)
  596. {
  597.         dum = ax;
  598.         ax = bx;
  599.         bx = dum;
  600.         dum = fb;
  601.         fb = fa;
  602.         fa = dum;
  603.     }
  604.     cx = bx + gold * (bx - ax);
  605.     fc = FunctionValue(cx);
  606. while (fb >= fc)
  607. {
  608.         r = (bx - ax) * (fb - fc);
  609.         q = (bx - cx) * (fb - fa);
  610.         dum = q - r;
  611.         if (Abs(dum) < tiny)
  612. {
  613. dum = tiny;
  614. }
  615.         u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
  616.         ulim = bx + glimit * (cx - bx);
  617.         if ((bx - u) * (u - cx) > 0)
  618. {
  619.             fu = FunctionValue(u);
  620.             if (fu < fc)
  621. {
  622.                 ax = bx;
  623.                 fa = fb;
  624.                 bx = u;
  625.                 fb = fu;
  626.                 return;
  627. }
  628.             else
  629. {
  630. if (fu > fb)
  631. {
  632. cx = u;
  633. fc = fu;
  634. return;
  635. }
  636.             }
  637.             u = cx + gold * (cx - bx);
  638.             fu = FunctionValue(u);
  639. }
  640.         else
  641. {
  642. if ((cx - u) * (u - ulim) > 0)
  643. {
  644. fu = FunctionValue(u);
  645. if (fu < fc)
  646. {
  647. bx = cx;
  648. cx = u;
  649. u = cx + gold * (cx - bx);
  650. fb = fc;
  651. fc = fu;
  652. fu = FunctionValue(u);
  653. }
  654.             }
  655. else
  656. {
  657. if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
  658. {
  659. u = ulim;
  660. fu = FunctionValue(u);
  661. }
  662. else
  663. {
  664. u = cx + gold * (cx - bx);
  665. fu = FunctionValue(u);
  666. }
  667. }
  668. }
  669. ax = bx;
  670. bx = cx;
  671. cx = u;
  672. fa = fb;
  673. fb = fc;
  674. fc = fu;
  675. }
  676. }
  678. //确定一维函数极小值点所在区间(重载)
  679. template <class _Ty>
  680. int MinimizerInterval(_Ty& ax, _Ty& bx, _Ty& cx, int& sign, 
  681. _Ty b, _Ty& fa, _Ty& fb, _Ty& fc)
  682. {
  683.     _Ty r,q,dum,gold = GoldNo + 1.0;
  684.     _Ty u,ulim,fu,tiny = 1e-20;
  685. int glimit = 100;
  687.     fa = FunctionValue(ax);
  688.     fb = FunctionValue(bx);
  689.     if (fb > fa)
  690. {
  691. sign = -1;
  692. fa *= sign;
  693. fb *= sign;
  694.     }
  695.     cx = bx + gold * (bx - ax);
  696. fc = sign * FunctionValue(cx);
  697. while (fb >= fc)
  698. {
  699.         r = (bx - ax) * (fb - fc);
  700.         q = (bx - cx) * (fb - fa);
  701.         dum = q - r;
  702.         if (Abs(dum) < tiny)
  703. {
  704. dum = tiny;
  705. }
  706.         u = bx - ((bx - cx) * q - (bx - ax) * r) / (2 * dum);
  707.         ulim = bx + glimit * (cx - bx);
  708. if(u > b) //超出f(x)右边界
  709. {
  710. cout<<"  The MinimizerInterval exceed right border. Program is over!"<<endl;
  711. return 0;
  712. }
  713.         if ((bx - u) * (u - cx) > 0)
  714. {
  715.             fu = sign * FunctionValue(u);
  716.             if (fu < fc)
  717. {
  718.                 ax = bx;
  719.                 fa = fb;
  720.                 bx = u;
  721.                 fb = fu;
  722.                 return 1;
  723. }
  724.             else
  725. {
  726. if (fu > fb)
  727. {
  728. cx = u;
  729. fc = fu;
  730. return 1;
  731. }
  732.             }
  733.             u = cx + gold * (cx - bx);
  734.             fu = sign * FunctionValue(u);
  735. }
  736.         else
  737. {
  738. if ((cx - u) * (u - ulim) > 0)
  739. {
  740. fu = sign * FunctionValue(u);
  741. if (fu < fc)
  742. {
  743. bx = cx;
  744. cx = u;
  745. u = cx + gold * (cx - bx);
  746. fb = fc;
  747. fc = fu;
  748. fu = sign * FunctionValue(u);
  749. }
  750.             }
  751. else
  752. {
  753. if ((u - ulim) * (ulim - cx) > 0 || FloatEqual((u - ulim) * (ulim - cx),0))
  754. {
  755. u = ulim;
  756. fu = sign * FunctionValue(u);
  757. }
  758. else
  759. {
  760. u = cx + gold * (cx - bx);
  761. fu = sign * FunctionValue(u);
  762. }
  763. }
  764. }
  765. ax = bx;
  766. bx = cx;
  767. cx = u;
  768. fa = fb;
  769. fb = fc;
  770. fc = fu;
  771. }
  772. return 1;
  773. }
  775. //1维函数极小值的黄金分割法
  776. template <class _Ty>
  777. _Ty ExtremumGold1D(_Ty ax, _Ty bx, _Ty cx, int sign, _Ty tol, _Ty& xmin)
  778. {
  779.     _Ty x0,x1,x2,x3,f0,f1,f2,f3,r = GoldNo;
  780.     _Ty c = 1.0 - GoldNo;
  781.     x0 = ax;
  782.     x3 = cx;
  783.     if (Abs(cx - bx) > Abs(bx - ax))
  784. {
  785. x1 = bx;
  786.         x2 = bx + c * (cx - bx);
  787. }
  788.     else
  789. {
  790.         x2 = bx;
  791.         x1 = bx - c * (bx - ax);
  792.     }
  793.     f1 = sign * FunctionValue(x1);
  794.     f2 = sign * FunctionValue(x2);
  795.     while (Abs(x3 - x0) > tol * (Abs(x1) + Abs(x2)))
  796. {
  797.         if (f2 < f1)
  798. {
  799.             x0 = x1;
  800.             x1 = x2;
  801.             x2 = r * x1 + c * x3;
  802.             f0 = f1;
  803.             f1 = f2;
  804.             f2 = sign * FunctionValue(x2);
  805. }
  806.         else
  807. {
  808.             x3 = x2;
  809.             x2 = x1;
  810.             x1 = r * x2 + c * x0;
  811.             f3 = f2;
  812.             f2 = f1;
  813.             f1 = sign * FunctionValue(x1);
  814.         }
  815.     }
  816.     if (f1 < f2)
  817. {
  818.         xmin = x1;
  819. if (sign == -1)
  820. f1 = sign * f1;
  821.         return f1;
  822. }
  823.     else
  824. {
  825.         xmin = x2;
  826. if (sign == -1)
  827. f2 = sign * f2;
  828.         return f2;
  829.     }
  830. }
  832. //1维函数极小值点的不用导数布伦特(Brent)法
  833. template <class _Ty>
  834. _Ty ExtremumBrentNonDerivative1D(_Ty ax, _Ty bx, _Ty cx, int sign, 
  835. _Ty tol, _Ty& xmin)
  836. {
  837.     int done,iter,itmax = 100;
  838.     _Ty d,fu,r,q,p,xm,tol1,tol2,a,b,cgold = 1.0 - GoldNo;
  839.     _Ty u,etemp,dum,v,w,x,e,fx,fv1,fw,zeps = DOUBLEERROR;
  840.     a = ax;
  841.     if (cx < ax)
  842. {
  843. a = cx;
  844. }
  845.     b = ax;
  846.     if (cx > ax)
  847. {
  848. b = cx;
  849. }
  850.     v = bx;
  851.     w = v;
  852.     x = v;
  853.     e = 0.0;
  854.     fx = sign * FunctionValue(x);
  855.     fv1 = fx;
  856.     fw = fx;
  857.     for (iter = 1; iter<=itmax; iter++)
  858. {
  859.         xm = 0.5 * (a + b);
  860.         tol1 = tol * Abs(x) + zeps;
  861.         tol2 = 2.0 * tol1;
  862.         if (Abs(x - xm) <= tol2 - 0.5 * (b - a))
  863. {
  864. break;
  865. }
  866.         done = -1;
  867.         if (Abs(e) > tol1)
  868. {
  869.             r = (x - w) * (fx - fv1);
  870.             q = (x - v) * (fx - fw);
  871.             p = (x - v) * q - (x - w) * r;
  872.             q = 2.0 * (q - r);
  873.             if (q > 0.0)
  874. {
  875. p = -p;
  876. }
  877.             q = Abs(q);
  878.             etemp = e;
  879.             e = d;
  880.             dum = Abs(0.5 * q * etemp);
  881.             if (Abs(p) < dum && p > q * (a - x) && p < q * (b - x))
  882. {
  883.                 d = p / q;
  884.                 u = x + d;
  885.                 if (u - a < tol2 || b - u < tol2)
  886. {
  887.                     d = Abs(tol1) * Sgn(xm - x);
  888.                 }
  889.                 done = 0;
  890.             }
  891.         }
  892.         if (done)
  893. {
  894.             if (x > xm || FloatEqual(x,xm))
  895. {
  896.                 e = a - x;
  897. }
  898.             else
  899. {
  900.                 e = b - x;
  901.             }
  902.             d = cgold * e;
  903.         }
  904.         if (Abs(d) > tol1 || FloatEqual(d,tol1))
  905. {
  906.             u = x + d;
  907. }
  908.         else
  909. {
  910.             u = x + Abs(tol1) * Sgn(d);
  911.         }
  912.         fu = sign * FunctionValue(u);
  913.         if (fu < fx || FloatEqual(fu,fx))
  914. {
  915.             if (u > x || FloatEqual(u,x))
  916. {
  917.                 a = x;
  918. }
  919.             else
  920. {
  921.                 b = x;
  922.             }
  923.             v = w;
  924.             fv1 = fw;
  925.             w = x;
  926.             fw = fx;
  927.             x = u;
  928.             fx = fu;
  929. }
  930.         else
  931. {
  932.             if (u < x)
  933. {
  934.                 a = u;
  935. }
  936.             else
  937. {
  938.                 b = u;
  939.             }
  940.             if (fu < fw || FloatEqual(fu,fw) || FloatEqual(w,x))
  941. {
  942.                 v = w;
  943.                 fv1 = fw;
  944.                 w = u;
  945.                 fw = fu;
  946. }
  947.             else
  948. {
  949. if (fu < fv1 || FloatEqual(fu,fv1) ||  FloatEqual(v,x) ||  FloatEqual(v,w))
  950. {
  951. v = u;
  952. fv1 = fu;
  953. }
  954.             }
  955.         }
  956.     }
  957.     if (iter > itmax)
  958. {
  959. cout<<" ExtremumBrentNonDerivative1D exceed maximum iterations."<<endl;
  960. }
  961.     xmin = x;
  962.     if(sign == -1) 
  963. {
  964. fx *= sign;
  965. }
  966. return fx;
  967. }
  969. #endif // _EXTREMUM_INL