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

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

数值算法/人工智能

开发平台：

Visual C++

NonLinearEquation.inl：源码内容
  1. //NonLinearEquation.inl  非线性方程(组)求解函数(方法)定义
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 何渝, 2002
  4. // 最后修改: 2002.5.31.
  6. #ifndef _NONLINEAREQUATION_INL
  7. #define _NONLINEAREQUATION_INL
  9. //用二分法搜索方程f(x)=0在区间[a,b]内的全部实根
  10. template <class _Ty>
  11. inline size_t 
  12. RootHalves(_Ty a, _Ty b, _Ty step, _Ty eps, valarray<_Ty>& x, size_t m)
  13. {
  14.     int n(0), js;
  15.     _Ty z=a, z1, y, y1, z0, y0;
  16.    
  17. y = FunctionValueRH(z);
  18.     while((z<=b+step/2.0)&&(n!=m))
  19. { 
  20.   if(Abs(y) < eps)
  21.       {
  22. n = n + 1;
  23. x[n-1] = z;
  24. z = z + step / 2.0;
  25. y=FunctionValueRH(z);
  26.       }
  27.       else
  28.       {
  29. z1=z+step;
  30. y1=FunctionValueRH(z1);
  31.         if(Abs(y1)<eps)
  32.         {
  33. n=n+1;
  34. x[n-1]=z1;
  35.             z=z1+step/2.0;
  36. y=FunctionValueRH(z);
  37.         }
  38.         else
  39. if(y*y1>0.0)
  40.             {
  41. y=y1;
  42. z=z1;
  43. }
  44.             else
  45.             {
  46. js=0;
  47.                 while(js==0)
  48.                 {
  49. if(Abs(z1-z)<eps)
  50.                     {
  51. n=n+1;
  52. x[n-1]=(z1+z)/2.0;
  53.                         z=z1+step/2.0;
  54. y=FunctionValueRH(z);
  55.                         js=1;
  56.                     }
  57.                     else
  58.                     {
  59. z0=(z1+z)/2.0;
  60. y0=FunctionValueRH(z0);
  61.                         if(Abs(y0)<eps)
  62.                         {
  63. x[n]=z0; n=n+1; js=1;
  64.                             z=z0+step/2.0;
  65. y=FunctionValueRH(z);
  66.                         }
  67.                         else
  68. if((y*y0) < 0.0)
  69. {
  70. z1 = z0;
  71. y1 = y0;
  72. }
  73. else
  74. {
  75. z = z0;
  76. y = y0;
  77. }
  78.                       }
  79.                   }
  80.               }
  81.           }
  82. }
  83.     return(n); //返回搜索到的实根个数
  84. }
  86. //牛顿(Newton)法求解非线性方程一个实根
  87. template <class _Ty>
  88. inline int 
  89. RootNewton( _Ty& x, _Ty eps, size_t js)
  90. { 
  91.     int k = 1, r = js;
  92.     _Ty  x1;
  93. _Ty x0 = x;
  94. valarray<_Ty> y(2);
  96.     FunctionValueRN(x0, y);
  97.     _Ty d=eps + 1.0;
  98.     while((d > eps || FloatEqual(d, eps))&&(r != 0))
  99.     {
  100. if(FloatEqual(y[1],0)) return(-1);
  101.         x1 = x0 - y[0] / y[1];
  102.         FunctionValueRN(x1,y);
  103.         d=Abs(x1 - x0);
  104. _Ty p=Abs(y[0]);
  105.         if(p > d) d = p;
  106.         x0 = x1;
  107. r -= 1;
  108.     }
  109.     x = x1;
  110.     k = js - r;
  111.     return(k);
  112. }
  114. //埃特金(Aitken)法求解非线性方程一个实根
  115. template <class _Ty>
  116. inline int 
  117. RootAitken( _Ty& x, _Ty eps, size_t js)
  118. {
  119.     _Ty u, v, x0;
  120.     int r(0), flag(0);
  122. x0 = x;
  124.     while((flag == 0) && (r != js))
  125.     {
  126. r++; 
  127.         u = FunctionValueRA(x0);
  128. v = FunctionValueRA(u);
  129.         if(Abs(u-v)<eps)
  130. {
  131. x0 = v;
  132. flag = 1;
  133. }
  134.         else x0=v-(v-u)*(v-u)/(v-2.0*u+x0);
  135.     }
  136.     x = x0;
  137. r = js - r;
  138.     return(r);
  139. }
  141. //连分式(Fraction)法求解非线性方程一个实根
  142. template <class _Ty>
  143. inline int 
  144. RootFraction(_Ty& x, _Ty eps)
  145. {
  146. int r(10), i, j, m, it;
  147.     _Ty x0, q(1.0e+35), h(0), z;
  148. valarray<_Ty> a(10), y(10);
  150. x0 = x; 
  152.     while(r!=0)
  153.     { 
  154. r = r - 1; 
  155. j = 0; 
  156. it = r;
  157.         while(j<=7)
  158.         {
  159. if(j<=2)
  160. z = x0 + 0.1 * j;
  161.             else
  162. z = h;
  163.             y[j] = FunctionValueRF(z);
  164.             h = z;
  165.             if(j==0)
  166. a[0] = z;
  167.             else
  168.             {
  169. m = 0;
  170. i = 0;
  171.                 while((m==0)&&(i<j))
  172.                 {
  173. if(FloatEqual((h-a[i]),0))
  174. m = 1;
  175.                     else
  176. h=(y[j]-y[i])/(h-a[i]);
  177.                     i++;
  178.                 }
  179.                 a[j] = h;
  180.                 if(m!=0) a[j]=q;
  181.                 h=0.0;
  182.                 for(i=j-1; i>=0; i--)
  183.                 {
  184. if(FloatEqual((a[i+1]+h),0))
  185. h = q;
  186.                     else
  187. h=-y[i]/(a[i+1]+h);
  188.                 }
  189.                 h = h + a[0];
  190.             }
  191.             if(Abs(y[j])>=eps) j++;
  192.             else 
  193. {
  194. j = 10;
  195. r = 0;
  196. }
  197.         }
  198.         x0 = h;
  199.     }
  200.     x = h;
  201.     return(10-it);
  202. }
  204. //QR法求代数方程全部根
  205. template <class _Ty, class _Tz>
  206. inline int 
  207. RootQR(valarray<_Ty>& a, valarray<_Tz>& x, _Ty eps, size_t jt)
  208. {
  209. int i, j, n;
  211.     n = x.size(); //多项式方程的次数
  213. matrix<_Ty> q(n, n);
  215.     for(j=0; j<n; j++)
  216. q(0,j) = -a[n-j-1] / a[n];
  218.     for(j=1; j<n; j++)
  219. for(i=0; i<n; i++)
  220. q(j,i) = 0;
  222. for(i=0; i<n-1; i++)
  223.       q(i+1, i) = 1;
  225. i = EigenvalueVectorHessenbergQR(q,x,eps,jt); //求全部特征值的QR方法
  226.     
  227. return(i);
  228. }
  230. //牛顿下山(NewtonHillDown)法求解实系数代数方程全部根(实根和复根)
  231. template <class _Ty, class _Tz>
  232. inline int 
  233. RootNewtonHillDown(valarray<_Ty>& a, valarray<_Tz>& cx)
  234. {
  235. int m, i, jt, k, is, it, n;
  236.     _Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1;
  237. _Ty w, u, p, q, pq, v, u1, v1;
  238.   
  239.     n = cx.size(); //多项式方程的次数
  240.     m = n;
  242.     while((m > 0) && (FloatEqual(a[m], 0.0))) m--;
  244.     if(m <= 0) return(-1);
  245.     
  246. for(i = 0; i <= m; i++) a[i] /= a[m];
  248.     for(i = 0; i <= m / 2; i++)
  249.     {
  250. w = a[i];
  251. a[i] = a[m - i];
  252. a[m - i] = w;
  253. }
  255.     k = m;
  256. is = 0;
  257. w = 1.0;
  258.     jt = 1;
  259.     while(jt == 1)
  260.     {
  261. pq = Abs(a[k]);
  262. while(pq < LONGDOUBLEERROR)
  263.         {
  264. cx[k - 1] = complex<_Ty>(0.0, 0.0);
  265. k = k - 1;
  266.             if(k == 1)
  267.             {
  268. cx[0] = complex<_Ty>(-a[1] * w / a[0], 0.0);
  269.                 return(1);
  270.             }
  271.             pq = Abs(a[k]);
  272.         }
  273. q = log(pq);
  274. q = q / (1.0 * k);
  275. q = exp(q);
  276.         p = q;
  277. w = w * p;
  278.         for(i = 1; i <= k; i++)
  279.         {
  280. a[i] /= q;
  281. q *= p;
  282. }
  283.         x = 0.0001;
  284. x1 = x; 
  285. y = 0.2;
  286. y1 = y; 
  287. dx = 1.0;
  288.         g = 1.0e+37; 
  289.  
  290. zjn:    u = a[0];
  291. v = 0.0;
  292.         for(i = 1; i <= k; i++)
  293.         { 
  294. p = u * x1;
  295. q = v * y1;
  296.             pq = (u + v) * (x1 + y1);
  297.             u = p - q + a[i];
  298. v = pq - p - q;
  299.         }
  300.         g1 = u * u + v * v;
  302.         if(g1 >= g)
  303.         { 
  304. if(is != 0)
  305. { 
  306. it = 1;
  307. g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  308. if(it == 0)
  309. goto zjn;
  310. }
  311. else
  312. {
  313. g60(t, x, y, x1, y1, dx, dy, p, q, k, it);
  314. if(t >= 1.0e-03)  goto zjn;
  315. if(g > 1.0e-18)
  316. {
  317. it = 0;
  318. g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  319. if(it == 0) goto zjn;
  320. }
  322. }
  323. g90(cx , a, x, y, p, q, w, k);
  324.         }
  325.         else
  326.         {
  327. g = g1;
  328. x = x1;
  329. y = y1;
  330. is = 0;
  331.             if(g <= 1.0e-22)
  332. g90(cx, a, x, y, p, q, w, k);
  333.             else
  334.             {
  335. u1 = k * a[0]; 
  336. v1 = 0.0;
  337.                 for(i = 2; i <= k; i++)
  338.                 {
  339. p = u1 * x;
  340. q = v1 * y;
  341. pq = (u1 + v1) * (x + y);
  342.                     u1 = p - q + (k - i + 1) * a[i - 1];
  343.                     v1 = pq - p - q;
  344.                 }
  345.                 p = u1 * u1 + v1 * v1;
  346.                 if(p <= 1.0e-20)
  347.                 {
  348. it = 0;
  349.                     g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  350.                     if(it == 0) goto zjn;
  351.                     g90(cx, a, x, y, p, q, w, k);
  352.                 }
  353.                 else
  354.                 {
  355. dx = (u * u1 + v * v1) / p;
  356.                     dy = (u1 * v - v1 * u) / p;
  357.                     t = 1.0 + 4.0 / k;
  358.                     g60(t, x, y, x1, y1, dx, dy, p, q, k, it);
  359.                     if(t >= 1.0e-03) goto zjn;
  360.                     if(g > 1.0e-18)
  361.                     {
  362. it = 0;
  363.                         g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  364.                         if(it == 0) goto zjn;
  365.                     }
  366.                     g90(cx, a, x, y, p, q, w, k);
  367.                 }
  368.             }
  369.         }
  371.         if(k == 1) jt = 0;
  372.         else jt = 1;
  373.     }
  374.     return(1);
  375. }
  377. //牛顿下山(NewtonHillDown)法求解复系数代数方程全部根(实根和复根)
  378. //重载RootNewtonHillDown()
  379. template <class _Ty>
  380. inline int 
  381. RootNewtonHillDown(complex<_Ty> a[], valarray<complex<_Ty> >& cx)
  382. {
  383. int m, i, jt, k, is, it, n;
  384.     _Ty t, x, y, x1, y1, dx, dy, dd, dc, c, g, g1, mo, sb, xb;
  385. _Ty w, u, p, q, pq, v, u1, v1;
  386.     
  387. n = cx.size(); //多项式方程的次数
  388.     m = n;
  390.     while((m > 0) && (FloatEqual(Abs(a[m]), 0))) m--;
  392.     if(m <= 0) return(-1);
  394.     mo = Abs(a[m]);
  396. for(i = 0; i <= m; i++) a[i] /= mo;
  398.     for(i = 0; i <= m / 2; i++)
  399.     {
  400. swap(a[i], a[m-i]);
  401. }
  403.     k = m;
  404. is = 0;
  405. w = 1.0;
  406.     jt = 1;
  407.     while(jt == 1)
  408.     {
  409. pq = Abs(a[k]);
  410. while(pq < LONGDOUBLEERROR)
  411.         {
  412. cx[k - 1] = complex<_Ty>(0.0, 0.0);
  413. k--;
  414.             if(k == 1)
  415.             {
  416. mo = a[0].real() * a[0].real() + a[0].imag() * a[0].imag();
  417. sb = -w * (a[0].real() * a[1].real() + a[0].imag() * a[1].imag()) / mo;
  418. xb =  w * (a[1].real() * a[0].real() - a[0].imag() * a[1].imag()) / mo;
  419. cx[0] = complex<_Ty>(sb, xb);
  420.                 return(1); 
  421.             }
  422. pq = Abs(a[k]);
  423.         }
  425. q = log(pq);
  426. q = q / (1.0 * k);
  427. q = exp(q);
  428.         p = q;
  429. w = w * p;
  430.         
  431. for(i = 1; i <= k; i++)
  432.         {
  433. a[i] /= q;
  434. q *= p;
  435. }
  436.         
  437. x = 0.0001;
  438. x1 = x; 
  439. y = 0.2;
  440. y1 = y; 
  441. dx = 1.0;
  442.         g = 1.0e+37; 
  443.  
  444. zjn:   
  445. u = a[0].real();
  446. v = a[0].imag();
  447.         
  448. for(i = 1; i <= k; i ++)
  449.         { 
  450. p = u * x1;
  451. q = v * y1;
  452.             pq = (u + v) * (x1 + y1);
  453.             u = p - q + a[i].real();
  454. v = pq - p - q + a[i].imag();
  455.         }
  456.         g1 = u * u + v * v;
  458.         if(g1 >= g)
  459.         { 
  460. if(is != 0)
  461. { 
  462. it = 1;
  463. g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  464. if(it == 0)
  465. goto zjn;
  466. }
  467. else
  468. {
  469. gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
  470. if(t >= 1.0e-03)  goto zjn;
  471. if(g > 1.0e-18)
  472. {
  473. it = 0;
  474. g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  475. if(it == 0) goto zjn;
  476. }
  477. }
  479. gg90(cx, a, x, y, p, q, w, k);
  480.         }
  481.         else
  482.         {
  483. g = g1;
  484. x = x1;
  485. y = y1;
  486. is = 0;
  488.             if(g <= 1.0e-22)
  489. gg90(cx, a, x, y, p, q, w, k);
  490.             else
  491.             {
  492. u1 = k * a[0].real(); 
  493. v1 = a[0].imag();
  494.                 for(i = 2; i <= k; i++)
  495.                 {
  496. p = u1 * x;
  497. q = v1 * y;
  498. pq = (u1 + v1) * (x + y);
  499.                     u1 = p - q + (k - i + 1) * a[i - 1].real();
  500.                     v1 = pq - p - q + (k - i + 1) * a[i - 1].imag();
  501.                 }
  502.                 p = u1 * u1 + v1 * v1;
  503.                 if(p <= 1.0e-20)
  504.                 {
  505. it = 0;
  506.                     g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  507.                     if(it == 0) goto zjn;
  508.                     gg90(cx, a, x, y, p, q, w, k);
  509.                 }
  510.                 else
  511.                 {
  512. dx = (u * u1 + v * v1) / p;
  513.                     dy = (u1 * v - v1 * u) / p;
  514.                     t = 1.0 + 4.0 / k;
  515.                     gg60(t, x, y, x1, y1, dx, dy, p, q, k, it);
  516.                     if(t >= 1.0e-03) goto zjn;
  517.                     if(g > 1.0e-18)
  518.                     {
  519. it = 0;
  520.                         g65(x, y, x1, y1, dx, dy, dd, dc, c, k, is, it);
  521.                         if(it == 0) goto zjn;
  522.                     }
  523.                     gg90(cx, a, x, y, p, q, w, k);
  524.                 }
  525.             }
  526.         }
  527.         if(k == 1) jt = 0;
  528.         else jt = 1;
  529.     }
  530.     return 1;
  531. }
  533. //牛顿下山法(NewtonHillDown)子函数
  534. template <class _Ty>
  535. inline void
  536. g60(_Ty& t, _Ty& x, _Ty& y, _Ty& x1, _Ty& y1, _Ty& dx, _Ty& dy, 
  537. _Ty& p, _Ty& q, int& k, int& it)
  538. { 
  539. it = 1;
  540.     while(it == 1)
  541.     {
  542. t = t / 1.67;
  543. it = 0;
  544.         x1 = x - t * dx;
  545.         y1 = y - t * dy;
  546.         if(k >= 50)
  547. {
  548. p = sqrt(x1 * x1 + y1 * y1);
  549.             q = exp(85.0 / k);
  550.             if(p >= q) it = 1;
  551.         }
  552.     }
  553. }
  555. //牛顿下山法(NewtonHillDown)子函数
  556. template <class _Ty>
  557. inline void
  558. gg60(_Ty& t, _Ty& x, _Ty& y, _Ty& x1, _Ty& y1, _Ty& dx, _Ty& dy, 
  559. _Ty& p, _Ty& q, int& k, int& it)
  560. { 
  561. it = 1;
  562.     while(it == 1)
  563.     {
  564. t = t / 1.67;
  565. it = 0;
  566.         x1 = x - t * dx;
  567.         y1 = y - t * dy;
  568.         if(k >= 30)
  569. {
  570. p = sqrt(x1 * x1 + y1 * y1);
  571.             q = exp(75.0 / k);
  572.             if(p >= q) it = 1;
  573.         }
  574.     }
  575. return;
  576. }
  578. //牛顿下山法(NewtonHillDown)子函数
  579. template <class _Ty, class _Tz>
  580. inline void
  581. g90(valarray<_Tz> &cx, valarray<_Ty> &a, _Ty& x, _Ty& y, _Ty& p, _Ty& q, _Ty& w, int &k)
  582. {
  583. int i;
  584.     if(Abs(y) <= 1.0e-06)
  585.     {
  586. p = -x;
  587. y = 0.0;
  588. q = 0.0;
  589. }
  590.     else
  591.     {
  592. p = -2.0 * x;
  593. q = x * x + y * y;
  594.         cx[k - 1] = complex<_Ty>(x * w, -y * w);
  595.         k = k-1;
  596.     }
  597.     for(i = 1; i <= k; i ++)
  598.     {
  599. a[i] = a[i] - a[i - 1] * p;
  600.         a[i + 1] = a[i + 1] - a[i - 1] * q;
  601.     }
  602.     cx[k - 1] = complex<_Ty>(x * w, y * w);
  603.     k = k - 1;
  604.     if(k == 1)
  605. cx[0] = complex<_Ty>(-a[1] * w / a[0], 0.0);
  606. }
  608. //牛顿下山法(NewtonHillDown)子函数
  609. template <class _Ty, class _Tz>
  610. inline void
  611. gg90(valarray<_Tz> &cx, complex<_Ty> a[], _Ty& x, _Ty& y, _Ty& p, _Ty& q, _Ty& w, int &k)
  612. {
  613. int i;
  614.     double mo, sb, xb;    
  615. for (i = 1; i <= k; i ++)
  616. { 
  617.   xb = a[i].real() + a[i-1].real() * x - a[i - 1].imag() * y;
  618.   sb = a[i].imag() + a[i-1].real() * y + a[i - 1].imag() * x;
  619.   a[i] = complex<_Ty>(xb, sb); 
  620. }
  621.     
  622.     cx[k - 1] = complex<_Ty>(x * w, y * w); 
  623. k --;
  624. if(k == 1)
  625. {
  626. mo = Abs(a[0])*Abs(a[0]);
  627. sb = -w * (a[0].real() * a[1].real() + a[0].imag() * a[1].imag()) / mo;
  628. xb =  w * (a[1].real() * a[0].imag() - a[0].real() * a[1].imag()) / mo;
  629. cx[0] = complex<_Ty>(sb, xb);
  630. }
  631.     return;
  632. }
  634. //牛顿下山法(NewtonHillDown)子函数
  635. template <class _Ty>
  636. inline void
  637. g65(_Ty& x, _Ty& y, _Ty& x1, _Ty& y1, _Ty& dx, _Ty& dy, _Ty& dd, _Ty& dc, _Ty& c, 
  638. int& k, int& is, int& it)
  639. {
  640. if(it == 0)
  641.     {
  642. is = 1;
  643.         dd = sqrt(dx * dx + dy * dy);
  644.         if(dd > 1.0) dd = 1.0;
  645.         dc = 6.28 / (4.5 * k);
  646. c = 0.0;
  647.     }
  648.     while(true)
  649.     {
  650. c += dc;
  651.         dx = dd * cos(c); 
  652. dy = dd * sin(c);
  653.         x1 = x + dx; 
  654. y1 = y + dy;
  655.         if(c <= 6.29)
  656.         {
  657. it = 0;
  658. return;
  659. }
  660. dd /= 1.67;
  661. if(dd <= 1.0e-07)
  662. {
  663. it = 1;
  664. return;
  665. }
  666. c = 0.0;
  667. }
  668. }
  670. //梯度(Gradient)法(最速下降)求解非线性方程组一组实根
  671. template <class _Ty>
  672. inline int 
  673. RootGradient(_Ty eps, valarray<_Ty>& x, size_t js)
  674. { 
  675.     int r, j;
  676.     _Ty f, d, s;
  678. int n = x.size(); //方程的个数，也是未知数的个数
  679. std::valarray<_Ty> y(n); //定义一维数组对象y
  681.     r = js;
  682.     f = FunctionValueRG(x,y);
  683.     while(f>=eps)
  684.     {
  685. r = r - 1;
  686.         if(r==0) return(js);
  688.         d = 0;
  689.         
  690. for(j=0; j<n; j++) d=d+y[j]*y[j];
  692.         if(FloatEqual(d,0)) return(-1);
  694.         s = f / d;
  695.         for(j=0; j<n; j++) x[j]=x[j]-s*y[j];
  697.         f=FunctionValueRG(x, y);
  698.     }
  699.     return(js-r);
  700. }
  702. //拟牛顿(QuasiNewton)法求解非线性方程组一组实根
  703. template <class _Ty>
  704. inline int 
  705. RootQuasiNewton(_Ty eps, _Ty t, _Ty h, valarray<_Ty>& x, int k)
  706. {
  707.     int i, j, l;
  708.     _Ty am, z, beta, d;
  710. int n = x.size(); //方程组中各方程个数，也是未知数个数
  712. matrix<_Ty> a(n, n); //定义二维数组对象a
  713. valarray<_Ty> b(n); //定义一维数组对象b
  714. valarray<_Ty> y(n); //定义一维数组对象y
  715.     
  716. l = k;
  717. am = 1.0 + eps;
  718.     while(am>=eps)
  719.     {
  720. FunctionValueRSN(x, b);
  721.         am = 0.0;
  722.         for(i=0; i<n; i++)
  723.         {
  724. z = Abs(b[i]);
  725.             if(z>am) am = z;
  726.         }
  727.         if(am>=eps)
  728.         {
  729. l--;
  730.             if(l==0)
  731. {
  732. cout << "Fail!" << endl;
  733. return(0);
  734. }
  736.             for(j=0; j<n; j++)
  737.             {
  738. z = x[j];
  739. x[j] += h;
  740.                 FunctionValueRSN(x, y);
  741.                 for(i=0; i<n; i++) 
  742. a(i,j) = y[i];
  743.                 x[j] = z;
  744.             }
  745. if(LE_TotalChoiceGauss(a, b) == 0)
  746. {
  747. cout << "Fail!" << endl;
  748. return(-1); //矩阵a奇异
  749. }
  751.             beta = 1.0;
  753.             for(i=0; i<n; i++) beta -= b[i];
  754.             
  755. if(FloatEqual(beta,0))
  756. {
  757. cout << "Fail!" << endl;
  758. return(-2); //beta=0
  759. }
  761.             d = h / beta;
  762.             
  763. for(i=0; i<n; i++) x[i]=x[i]-d*b[i];
  765.             h = t * h;
  766.         }
  767.     }
  768.     return(k-l); //正常结束，返回迭代次数
  769. }
  771. //非线性方程组最小二乘解的广义逆法
  772. template <class _Ty>
  773. int RootLeastSquareGeneralizedInverse(int m, _Ty eps1, _Ty eps2, valarray<_Ty>& x, int ka)
  774. {
  775.     int i, j, k, l, kk, jt;
  776. bool yn;
  777.     _Ty y[10],b[10],alpha,z,h2,y1,y2,y3,y0,h1;
  779. int n = x.size(); //非线性方程组中未知数个数
  780.     
  781. matrix<_Ty> p(m, n);
  782.     matrix<_Ty> pp(n, m);
  783.     matrix<_Ty> u(m, m);
  784.     matrix<_Ty> v(n, n);
  786. valarray<_Ty> w(ka);
  787. valarray<_Ty> d(m);
  788. valarray<_Ty> dx(n);
  790.     l = 60; 
  791. alpha = 1.0;
  793.     while(l > 0)
  794. {
  795. FunctionValueRLSGI(x, d); //计算非线性方程组左边函数值
  796.         JacobiMatrix(x, p); //计算雅可比矩阵
  797.         
  798. //最小二乘的广义逆法求解线性方程组
  799. jt = LE_LinearLeastSquareGeneralizedInverse(p, d, dx, pp, eps2, u, v);
  801.         if(jt < 0) return(jt);
  802.         
  803. j = 0;
  804. jt = 1;
  805. h2 = 0.0;
  806.         
  807. while(jt == 1)
  808. {
  809. jt = 0;
  810.             
  811. if(j < 3) z = alpha + 0.01 * j;
  812.             else z = h2;
  813.             
  814. for(i = 0; i < n; i++) w[i] = x[i] - z * dx[i];
  815.             
  816. FunctionValueRLSGI(w, d);
  817.             y1 = 0.0;
  818.             
  819. for(i = 0; i < m; i++) y1 = y1 + d[i] * d[i];
  820.             for(i = 0; i < n; i++) w[i] = x[i] - (z + 0.00001) * dx[i];
  821.             
  822. FunctionValueRLSGI(w, d);
  823.             
  824. y2 = 0.0;
  825.             for(i = 0; i < m; i++) y2 = y2 + d[i] * d[i];
  826.             
  827. y0 =(y2 - y1) / 0.00001;
  828.            
  829. if(Abs(y0) > 1.0e-10)
  830. {
  831. h1 = y0;
  832. h2 = z;
  833.                 
  834. if(j == 0)
  835. { 
  836. y[0] = h1;
  837. b[0] = h2;
  838. }
  839.                 else
  840. {
  841. y[j] = h1;
  842. kk = k = 0;
  844.                     while((kk == 0) && (k < j))
  845. {
  846. y3 = h2 - b[k];
  848. yn = FloatEqual(y3, 0);
  849. if(yn) kk = 1;
  850. else h2 = (h1 - y[k]) / y3;
  852. k++;
  853. }
  854.                     
  855. b[j] = h2;
  856.                     
  857. if(kk != 0) b[j] = 1.0e+35;
  858.                     
  859. h2 = 0.0;
  860.                     
  861. for(k = j - 1; k >= 0; k--)
  862. h2 = -y[k] / (b[k + 1] + h2);
  863.                     h2 = h2 + b[0];
  864. }
  865.                 
  866. j++;
  867.                 
  868. if(j <= 7) jt = 1;
  869.                 else z = h2;
  870. }
  871. }
  872.         
  873. alpha = z;
  874. y1 = y2 = 0.0;
  875.        
  876. for(i = 0; i <= n-1; i ++)
  877. {
  878. dx[i] = -alpha * dx[i];
  879.             x[i] = x[i] + dx[i];
  880.             y1 = y1 + Abs(dx[i]);
  881.             y2 = y2 + Abs(x[i]);
  882. }
  883.         
  884. if(y1 < eps1 * y2) return(1);
  886.         l--;
  887. }
  889.     return(0);
  890. }
  892. //蒙特卡洛(MonteCarlo)法求解f(x)=0的一个实根
  893. //f(x)的自变量为与系数都为实数
  894. template <class _Ty>
  895. inline void 
  896. RootMonteCarloReal(_Ty& x, _Ty b, int m, _Ty eps)
  897. {
  898.    // extern double mrnd1();
  899.     //extern double dmtclf();
  900.     //int k;
  901.     _Ty x1, y1;
  902.     _Ty a = b;
  903. size_t k = 1;
  904. double r = 1.0;
  905. _Ty xx = x;
  906. _Ty y = FunctionValueMCR(xx); //计算函数值
  908.     while(a>eps||FloatEqual(a,eps))
  909.     {
  910. x1 = rand_01_One(r); //取随机数
  911. x1 = -a + 2.0 * a * x1;
  912.         x1 = xx + x1; 
  913. y1 = FunctionValueMCR(x1); //计算函数值
  914.         k++;
  915.         if(Abs(y1)>Abs(y)||FloatEqual(Abs(y1),Abs(y)))
  916.         {
  917. if(k>m)
  918. {
  919. k = 1; 
  920. a /= 2.0;
  921. }
  922. }
  923.         else
  924.         {
  925. k = 1;
  926. xx = x1;
  927. y = y1;
  928.             if(Abs(y) < eps)
  929.             {
  930. x=xx;
  931. exit(0);
  932. }
  933.         }
  934.     }
  935.     x = xx;
  936. }
  938. //蒙特卡洛(MonteCarlo)法求解f(x)=0的一个复根
  939. //f(x)的自变量为复数，或自变量与系数都为复数(不能都为实数)
  940. template <class _Tz, class _Ty>
  941. inline void 
  942. RootMonteCarloComplex(_Tz& cxy, _Ty b, int m, _Ty eps)
  943. {
  944.     size_t k = 1;
  945. _Tz xxyy(cxy);
  946.     _Ty a = b;
  947. double r(1);
  948. _Ty z, z1;
  950.     z = FunctionModule(xxyy);
  952.     while(a > eps || FloatEqual(a,eps))
  953.     {
  954. _Ty tempx = -a + 2.0 *a * rand_01_One(r);
  955. _Ty tempy = -a + 2.0 *a * rand_01_One(r);
  957. _Tz x1y1(tempx,tempy);
  958. x1y1 += xxyy;
  960.         z1 = FunctionModule(x1y1);
  961.         
  962. k++;
  964.         if(z1 > z || FloatEqual(z1,z))
  965.         {
  966. if(k > m)
  967. {
  968. k = 1;
  969. a = a / 2.0;
  970. }
  971. }
  972.         else
  973.         {
  974. k = 1;
  975. xxyy = x1y1;
  976. z = z1;
  977.             if(z < eps)
  978. {
  979. cxy = xxyy;
  980. exit(0);
  981. }
  982.           }
  983.       }
  984.     cxy = xxyy;
  985. }
  987. //蒙特卡洛(MonteCarlo)法求解f(x)=0的一组实根
  988. //f(x)的自变量为与系数都为实数
  989. template <class _Ty>
  990. inline void 
  991. RootMonteCarloGroupReal(valarray<_Ty>& x, _Ty b, int m, _Ty eps)
  992. {
  993.     _Ty a=b;
  994. size_t k=1;
  995. double r=1.0;
  997. int n = x.size(); //方程个数，也是未知量的个数
  998. valarray<_Ty> y(n);
  1000. _Ty z = FunctionModule(x);
  1002.     while(a>eps||FloatEqual(a,eps))
  1003.     {
  1004. for(size_t i = 0; i < n; i++)
  1005. y[i] = -a + 2.0 * a * rand_01_One(r) + x[i];
  1007.         _Ty z1 = FunctionModule(y);
  1008.         
  1009. k++;
  1011.         if(z1 > z || FloatEqual(z1,z))
  1012.         {
  1013. if(k > m)
  1014. {
  1015. k = 1;
  1016. a = a / 2.0;
  1017. }
  1018. }
  1019.         else
  1020.         {
  1021. k = 1; 
  1023.             for(i = 0; i < n; i++) x[i] = y[i];
  1025.             z = z1;
  1026.             
  1027. if(z < eps) exit(0);
  1028.         }
  1029.     }
  1030. }
  1032. #endif //_NONLINEAREQUATION_INL