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

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

数值算法/人工智能

开发平台：

Visual C++

Integral.inl：源码内容
  1. // Integral.inl 数值积分头文件
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 2002
  4. // 最后修改: 2002.5.31.
  6. #ifndef _INTEGRAL_INL //避免多次编译
  7. #define _INTEGRAL_INL
  9. //变步长梯形法求积
  10. template <class _Ty>
  11. _Ty IntegralTrapezia(_Ty a, _Ty b, _Ty eps)
  12. {
  13. _Ty t;
  14.     _Ty fa = FunctionValueIT(a);
  15. _Ty fb = FunctionValueIT(b);
  16.     int n(1);
  17. _Ty h = b - a;
  18.     _Ty t1 = h * (fa + fb) / 2.0;
  19.     _Ty p = eps + 1.0;
  20.     while(p > eps || FloatEqual(p, eps))
  21.     {
  22. _Ty s(0);
  23.         for(int k=0; k<n; k++)
  24.         {
  25. _Ty x = a + (k + 0.5) * h;
  26.             s = s + FunctionValueIT(x);
  27.           }
  28.         t = (t1 + h * s) / 2.0;
  29.         p = Abs(t1 - t);
  30.         t1 = t;
  31. n = n + n;
  32. h = h / 2.0;
  33.     }
  34.     return(t);
  35. }
  37. //变步长辛卜生法求积
  38. template <class _Ty>
  39. _Ty IntegralSimpson1D(_Ty a, _Ty b, _Ty eps)
  40. {
  41.     _Ty s2;
  42.     int n(1);
  43. _Ty h = b - a; //步长
  44.     _Ty t1 = h * (FunctionValueIS1D(a) + FunctionValueIS1D(b)) / 2.0;
  45.     _Ty s1 = t1;
  46.     _Ty ep = eps + 1.0;
  47.     while(ep > eps || FloatEqual(ep, eps))
  48.     {
  49. _Ty p(0);
  50.         for(int k=0; k<n; k++)
  51.         {
  52. _Ty x = a + (k + 0.5) * h;
  53.             p = p + FunctionValueIS1D(x);
  54.         }
  55.         _Ty t2 = (t1 + h * p) / 2.0;
  56.         s2 = (4.0 * t2 - t1) / 3.0;
  57.         ep = Abs(s2 - s1);
  58.         t1 = t2;
  59. s1 = s2;
  60. n = n + n;
  61. h = h / 2.0;
  62.     }
  63.     return(s2);
  64. }
  66. //自适应梯形法求积
  67. template <class _Ty>
  68. _Ty IntegralTrapeziaSelfAdapt(_Ty a, _Ty b, _Ty eps, _Ty d)
  69. {
  70.     valarray<_Ty> t(2);
  71.     _Ty h = b - a;
  72. t[0]=0.0;
  73.     _Ty f0 = FunctionValueITSA(a);
  74. _Ty f1 = FunctionValueITSA(b);
  75.     _Ty t0 = h * (f0 + f1) / 2.0;
  76.     ITSA(a,b,h,f0,f1,t0,eps,d,t);
  77.     _Ty z = t[0];
  78. return(z);
  79. }
  81. //自适应梯形法求积辅助函数
  82. template <class _Ty>
  83. int ITSA(_Ty x0, _Ty x1, _Ty h, _Ty f0, _Ty f1, _Ty t0, 
  84. _Ty eps, _Ty d, valarray<_Ty>& t)
  85. { 
  86.     _Ty x = x0 + h / 2.0;
  87. _Ty f = FunctionValueITSA(x);
  88.     _Ty t1=h*(f0+f)/4.0;
  89. _Ty t2=h*(f+f1)/4.0;
  90.     _Ty p=Abs(t0-(t1+t2));
  91.     if((p<eps)||(h/2.0<d))
  92. {
  93. t[0]=t[0]+(t1+t2);
  94. return(0);
  95. }
  96.     else
  97.     {
  98. _Ty g=h/2.0; 
  99. _Ty eps1=eps/1.4;
  100.         ITSA(x0,x,g,f0,f,t1,eps1,d,t);
  101.         ITSA(x,x1,g,f,f1,t2,eps1,d,t);
  102.         return(1);
  103.     }
  104. }
  106. //龙贝格法求积
  107. template <class _Ty>
  108. _Ty IntegralRomberg(_Ty a, _Ty b, _Ty eps)
  109. {
  110.     _Ty y[10], q;
  111.     _Ty h = b - a; //步长
  112.     y[0] = h * (FunctionValueR(a) + FunctionValueR(b)) / 2.0;
  113.     int m(1), n(1);
  114. _Ty ep=eps+1.0;
  115.     while((ep>eps || FloatEqual(ep,eps)) && (m <= 9))
  116.     {
  117. _Ty p(0);
  118.         for(int i=0;i<n;i++)
  119.         {
  120. _Ty x = a + (i + 0.5) * h;
  121.             p = p + FunctionValueR(x);
  122.         }
  123.         p = (y[0] + h * p) / 2.0;
  124.         _Ty s(1);
  125.         for(int k=1; k<=m; k++)
  126.         {
  127. s = 4.0 * s;
  128.             q = (s * p - y[k-1]) / (s - 1.0);
  129.             y[k-1] = p;
  130. p = q;
  131.         }
  132.         ep = Abs(q - y[m-1]);
  133.         m = m + 1;
  134. y[m-1] = q; 
  135. n = n + n;
  136. h = h / 2.0;
  137.     }
  138.     return(q);
  139. }
  141. //一维连分式法求积
  142. template <class _Ty>
  143. _Ty IntegralFraction1D(_Ty a, _Ty b, _Ty eps)
  144. {
  145.     valarray<_Ty> h(10), bb(10);
  146.     int m(1), n(1), k, l, j;
  147.     _Ty g, hh = b - a;
  148. h[0] = hh;
  149.     _Ty t1 = hh * (FunctionValueF1D(a) + FunctionValueF1D(b)) / 2.0;
  150.     _Ty s1 = t1;
  151. bb[0]=s1;
  152. _Ty ep = 1.0 + eps;
  153.     while((ep>eps || FloatEqual(ep,eps)) && (m <= 9))
  154.     { 
  155. _Ty s(0);
  156.         for(k=0; k<n; k++)
  157.         {
  158. _Ty x = a + (k + 0.5) * hh;
  159.             s = s + FunctionValueF1D(x);
  160.         }
  161.         _Ty t2 = (t1 + hh * s) / 2.0;
  162.         m++;
  163.         h[m-1] = h[m-2] / 2.0;
  164.         g = t2;
  165.         l = 0;
  166. j = 2;
  167.         while((l==0) && (j<=m))
  168.         {
  169. s = g - bb[j-2];
  170.             if(FloatEqual(s,0)) l = 1;
  171.             else g = (h[m-1] - h[j-2]) / s;
  172.             j++;
  173.         }
  174.         bb[m-1] = g;
  175.         if(l!=0) bb[m-1] = 1.0e+35;
  176.         g = bb[m-1];
  177.         for(j=m; j>=2; j--) g = bb[j-2] - h[j-2] / g;
  178.         ep = Abs(g-s1);
  179.         s1 = g;
  180. t1 = t2;
  181. hh = hh / 2.0;
  182. n = n + n;
  183.     }
  184.     return(g);
  185. }
  187. //高振荡函数法求积
  188. template <class _Ty>
  189. void IntegralSurge(_Ty a, _Ty b, int m, valarray<_Ty>& fa, 
  190. valarray<_Ty>& fb, valarray<_Ty>& s)
  191. {
  192.     _Ty sa[4],sb[4],ca[4],cb[4];
  194. int n = fa.size(); //给定积分区间两端点上的导数最高阶数+1
  195. int mm=1;
  197.     _Ty sma=sin(m*a);
  198. _Ty smb=sin(m*b);
  199.     _Ty cma=cos(m*a);
  200. _Ty cmb=cos(m*b);
  201.     
  202. sa[0] = ca[3] = sma;
  203. sa[1] = ca[0] = cma;
  204. sa[2] = ca[1] = -sma; 
  205. sa[3] = ca[2] = -cma;
  206.     sb[0] = cb[3] = smb;
  207. sb[1] = cb[0] = cmb;
  208. sb[2] = cb[1] = -smb;
  209. sb[3] = cb[2] = -cmb;
  211.     s[0] = s[1] = 0.0; 
  213.     for(int k=0; k<n; k++)
  214.     {
  215. int j = k;
  216.         while(j>=4) j = j - 4;
  217.         mm = mm * m;
  218.         s[0] = s[0]+(fb[k]*sb[j]-fa[k]*sa[j])/(1.0*mm);
  219.         s[1] = s[1]+(fb[k]*cb[j]-fa[k]*ca[j])/(1.0*mm);
  220.     }
  221.     
  222. s[1] = -s[1];
  223. }
  225. //勒让德-高斯法求积
  226. template <class _Ty>
  227. _Ty IntegralLegendreGauss(_Ty a, _Ty b, _Ty eps)
  228. {
  229.     _Ty t[5] = 
  230. {
  231. -0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459
  232. };
  233.     
  234. _Ty c[5] = 
  235. {
  236. 0.2369268851, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268851
  237. };
  239.     int m(1);
  240. _Ty g;
  241.     _Ty h = b - a;
  242. _Ty s = Abs(0.001 * h);
  243.     _Ty p = 1.0e+35;
  244. _Ty ep = eps + 1.0;
  246.     while((ep > eps || FloatEqual(ep,eps)) && (Abs(h) > s))
  247.     {
  248. g = 0;
  249.         for(int i=1; i<=m; i++)
  250.         {
  251. _Ty aa = a + (i - 1.0) * h;
  252. _Ty bb = a + i * h;
  253.             _Ty w(0);
  254.             for(int j=0; j<5; j++)
  255.             {
  256. _Ty x = ((bb - aa) * t[j] + (bb + aa)) / 2.0;
  257.                 w = w + FunctionValueLegG(x) * c[j];
  258.             }
  259.             g += w;
  260.         }
  261.         g = g * h / 2.0;
  262.         ep = Abs(g - p) / (1.0 + Abs(g));
  263.         p = g;
  264. m++;
  265. h = (b - a) / m;
  266.     }
  267.     return(g);
  268. }
  270. //拉盖尔-高斯法求积
  271. template <class _Ty>
  272. void IntegralLaguerreGauss(_Ty& dValue)
  273. {
  274. _Ty x;
  276. dValue = 0.0;
  278.     long double t[5] = 
  279. { 
  280. 0.26355990, 1.41340290, 3.59642600, 7.08580990, 12.64080000
  281. };
  282.     
  283. long double c[5] = 
  284. {
  285. 0.6790941054, 1.638487956, 2.769426772, 4.315944000, 7.104896230
  286. };
  287.     
  288. for(int i=0; i<5; i++)
  289.     {
  290. x = t[i];
  291. dValue = dValue + c[i] * FunctionValueLagG(x); 
  292. }
  293. }
  295. //埃尔米特-高斯法求积
  296. template <class _Ty >
  297. void IntegralHermiteGauss(_Ty& dValue)
  298. {
  299.     _Ty t[5] = 
  300. {
  301. -2.02018200, -0.95857190, 0.0, 0.95857190, 2.02018200
  302. };
  303.     
  304. _Ty c[5] = 
  305. {
  306. 1.181469599, 0.9865791417,0.9453089237, 0.9865791417,1.181469599
  307. };
  308.     
  309. dValue = 0;
  310.     
  311. for(int i=0; i<5; i++)
  312.     {
  313. _Ty x = t[i];
  314. dValue = dValue + c[i] * FunctionValueHG(x);
  315. }
  316. }
  318. //切比雪夫法求积
  319. template <class _Ty >
  320. _Ty IntegralChebyshev(_Ty a, _Ty b, _Ty eps)
  321. {
  322. _Ty t[5]={-0.8324975, -0.3745414, 0.0, 0.3745414, 0.8324975};
  323. int m(1);
  324. _Ty g;
  325.     _Ty h = b - a;
  326. _Ty d = Abs(0.001 * h);
  327.     _Ty p = 1.0e+35;
  328. _Ty ep = eps + 1.0;
  329.     while((ep > eps || FloatEqual(ep,eps)) && (Abs(h) > d))
  330.     {
  331. g = 0.0;
  332.         for(int i=1; i<=m; i++)
  333.         {
  334. _Ty aa = a + (i - 1.0) * h;
  335. _Ty bb = a + i * h;
  336.             _Ty s(0);
  337. for(int j=0; j<5; j++)
  338.             {
  339. _Ty x = ((bb - aa) * t[j] + (bb + aa)) / 2.0;
  340.                 s = s + FunctionValueCb(x);
  341.             }
  342.             g += s;
  343.         }
  344.         g = g * h / 5.0;
  345.         ep = Abs(g - p) / (1.0 + Abs(g));
  346.         p = g;
  347. m++;
  348. h = (b - a) / m;
  349.     }
  350.     return(g);
  351. }
  353. //一维蒙特卡洛法求积
  354. template <class _Ty >
  355. _Ty IntegralMonteCarlo1D(_Ty a, _Ty b)
  356. {
  357.     _Ty r(1), s(0), d(10000);
  359.     for(int m=0; m<10000; m++)
  360.     {
  361. _Ty x = a + (b - a) * rand_01_One(r); //取随机数
  363.         s = s + FunctionValueMC1D(x) / d; //调用被积函数值
  364.     }
  365.     
  366. s = s * (b - a);
  367.     
  368. return(s);
  369. }
  371. //二重变步长辛卜生法求积
  372. template <class _Ty>
  373. _Ty IntegralSimpson2D(_Ty a, _Ty b, _Ty eps)
  374. {
  375. int n(1);
  376. _Ty h = 0.5 * (b - a); //步长
  377. _Ty d = Abs((b - a) * 1.0e-06);
  378.     _Ty s1 = IntegralSimp2(a, eps); //调用IntegralSimp2函数
  379. _Ty s2 = IntegralSimp2(b, eps);
  380.     _Ty t1 = h * (s1 + s2);
  381.     _Ty s0(1.0e+35), s;
  382. _Ty ep = eps + 1.0;
  383.     while((ep>eps||FloatEqual(ep,eps))&&((Abs(h)>d)||(n<16)))
  384. {
  385. _Ty x = a - h;
  386. _Ty t2 = 0.5 * t1;
  387.         for(int j=1; j<=n; j++)
  388.         {
  389. x = x + 2.0 * h;
  390.             t2 = t2 + h * IntegralSimp2(x,eps);
  391.         }
  392.         s = (4.0 * t2 - t1) / 3.0;
  393.         ep = Abs(s - s0) / (1.0 + Abs(s));
  394.         n = n + n;
  395. s0 = s;
  396. t1 = t2; 
  397. h = h * 0.5;
  398.     }
  399.     return(s);
  400. }
  402. //二重变步长辛卜生法求积辅助函数
  403. template <class _Ty>
  404. _Ty IntegralSimp2(_Ty x, _Ty eps)
  405. {
  406.     valarray<_Ty> y(2);
  407. int n(1);
  408.     FunctionBoundaryIS2D(x, y);
  409.     _Ty h = 0.5 *(y[1] - y[0]);
  410.     _Ty d = Abs(h * 2.0e-06);
  411.     _Ty t1 = h * (FunctionValueIS2D(x, y[0]) + FunctionValueIS2D(x, y[1]));
  412.     _Ty ep = 1.0 + eps;
  413. _Ty g0(1.0e+35), g;
  414.     while((ep>eps||FloatEqual(ep,eps))&&((Abs(h)>d)||(n<16)))
  415.     {
  416. _Ty yy = y[0] - h;
  417.         _Ty t2 = 0.5 * t1;
  418.         for(int i=1; i<=n; i++)
  419.         {
  420. yy = yy + 2.0 * h;
  421.             t2 = t2 + h * FunctionValueIS2D(x, yy);
  422.         }
  423.         g = (4.0 * t2 - t1) / 3.0;
  424.         ep = Abs(g - g0) / (1.0 + Abs(g));
  425.         n = n + n;
  426. g0 = g;
  427. t1 = t2;
  428. h = 0.5 * h;
  429.     }
  430.     return(g);
  431. }
  433. //多重高斯法求积
  434. template <class _Ty >
  435. _Ty IntegralGaussMD(valarray<_Ty>& js)
  436. {
  437.     valarray<_Ty> y(2);
  438. _Ty p, s;
  439.     _Ty t[5] = 
  440. {
  441. -0.9061798459,-0.5384693101,0.0, 0.5384693101,0.9061798459
  442. };
  443.     
  444. _Ty c[5] = 
  445. {
  446. 0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851
  447. };
  449. int n = js.size(); //积分重数
  451. valarray<int> is(2*(n+1));
  452.     valarray<_Ty> x(n);
  453.     valarray<_Ty> a(2*(n+1));
  454.     valarray<_Ty> b(n+1);
  455. int m(1), l(1);
  456.     
  457. a[n] = a[2*n+1] = 1.0;
  459.     while(l==1)
  460.     {
  461. for(int j=m; j<=n; j++)
  462.         {
  463. FunctionBoundaryGMD(j-1,x,y);
  464.             a[j-1]=0.5*(y[1]-y[0])/js[j-1];
  465.             b[j-1]=a[j-1]+y[0];
  466.             x[j-1]=a[j-1]*t[0]+b[j-1];
  467.             a[n+j]=0.0;
  468.             is[j-1]=1; is[n+j]=1;
  469.         }
  470.         j = n;
  471. int q = 1;
  472.         while(q==1)
  473.         {
  474. int k = is[j-1];
  475.             if(j==n) p = FunctionValueGMD(x);
  476.             else p = 1.0;
  477.             a[n+j]=a[n+j+1]*a[j]*p*c[k-1]+a[n+j];
  478.             is[j-1] = is[j-1] + 1;
  479.             if(is[j-1] > 5)
  480.               if(is[n+j] >= js[j-1])
  481.               {
  482.   j = j - 1;
  483.   q = 1;
  484.                   if(j==0)
  485.                   {
  486.   s = a[n+1] * a[0]; 
  487.                       return(s);
  488.                   }
  489.               }
  490.               else
  491.               {
  492.   is[n+j]=is[n+j]+1;
  493.                   b[j-1]=b[j-1]+a[j-1]*2.0;
  494.                   is[j-1]=1; k=is[j-1];
  495.                   x[j-1]=a[j-1]*t[k-1]+b[j-1];
  496.                   if(j==n) q = 1;
  497.                   else q = 0;
  498.               }
  499.             else
  500.             {
  501. k = is[j-1];
  502.                 x[j-1] = a[j-1]*t[k-1]+b[j-1];
  503.                 if(j==n) q = 1;
  504.                 else q = 0;
  505.             }
  506.         }
  507.         m = j + 1;
  508.     }
  509. }
  511. //二重连分式法求积
  512. template <class _Ty >
  513. _Ty IntegralFraction2D(_Ty a, _Ty b, _Ty eps)
  514. {
  515.     _Ty bb[10],h[10],x,s;
  516. int m(1), n(1);
  517.     _Ty hh = b - a;
  518. h[0] = hh;
  519.     _Ty s1 = IntegralFraction2D2(a, eps);
  520. _Ty s2 = IntegralFraction2D2(b, eps);
  521.     _Ty t1 = hh * (s1 + s2) / 2.0;
  522.     _Ty s0 = t1;
  523. bb[0] = t1;
  524. _Ty ep = 1.0 + eps;
  525.     while((ep > eps || FloatEqual(ep,eps)) && ( m <= 9))
  526.     {
  527. _Ty t2 = 0.5 * t1;
  528.         for(int k=0; k<n; k++)
  529.         {
  530. x = a + (k + 0.5) * hh;
  531.             s1 = IntegralFraction2D2(x,eps);
  532.             t2 = t2 + 0.5 * s1 * hh;
  533.         }
  534.         m = m + 1;
  535.         h[m-1] = h[m-2] / 2.0;
  536.         _Ty g = t2;
  537. int l(0), j(2);
  538.         while((l==0)&&(j<=m))
  539.          {
  540. s = g - bb[j-2];
  541.             if(FloatEqual(s,0)) l = 1;
  542.             else g = (h[m-1] - h[j-2]) / s;
  543.             j = j + 1;
  544.         }
  545.         bb[m-1] = g;
  546.         if(l!=0) bb[m-1] = 1.0e+35;
  547.         s = bb[m-1];
  548.         for(j=m; j>=2; j--) s = bb[j-2] -h [j-2] / s;
  549.         ep = Abs(s-s0) / (1.0 + Abs(s));
  550.         n = n + n;
  551. t1 = t2; 
  552. s0 = s;
  553. hh = hh / 2.0;
  554.     }
  555.     return(s);
  556. }
  558. //二重连分式法求积辅助函数
  559. template <class _Ty>
  560. _Ty IntegralFraction2D2(_Ty x, _Ty eps)
  561. {
  562.     _Ty b[10],h[10],yy,s;
  563. valarray<_Ty> y(2);
  564.     int m(1), n(1);
  565.     FunctionBoundaryF2D(x, y);
  566.     _Ty hh = y[1] - y[0];
  567. h[0] = hh;
  568.     _Ty t1=0.5*hh*(FunctionValueF2D(x,y[0])+FunctionValueF2D(x,y[1]));
  569.     _Ty s0 = t1;
  570. b[0] = t1;
  571. _Ty ep = 1.0 + eps;
  572.     while((ep > eps || FloatEqual(ep, eps)) && (m <= 9))
  573.     {
  574. _Ty t2 = 0.5 * t1;
  575.         for(int k=0; k<n; k++)
  576.         {
  577. yy = y[0] +(k + 0.5) * hh;
  578.             t2 = t2 + 0.5 * hh * FunctionValueF2D(x, yy);
  579. }
  580.         m = m + 1;
  581.         h[m-1] = h[m-2] / 2.0;
  582.         _Ty g = t2;
  583. int l(0), j(2);
  584.         while((l==0)&&(j<=m))
  585.         {
  586. s = g - b[j-2];
  587.             if(FloatEqual(s,0)) l = 1;
  588.             else g = (h[m-1] - h[j-2]) / s;
  589.             j = j + 1;
  590.         }
  591.         b[m-1] = g;
  592.         if(l!=0) b[m-1] = 1.0e+35;
  593.         s = b[m-1];
  594.         for(j=m; j>=2; j--) s = b[j-2] - h[j-2] / s;
  595.         ep = Abs(s - s0) / (1.0 + Abs(s));
  596.         n = n + n;
  597. t1 = t2;
  598. s0 = s;
  599. hh = 0.5 * hh;
  600.     }
  601.     return(s);
  602. }
  604. //多重蒙特卡洛法求积
  605. template <class _Ty >
  606. _Ty IntegralMonteCarlo2D(valarray<_Ty>& a, valarray<_Ty>& b)
  607. {
  608. int i;
  609. int n = a.size(); //积分的重数
  611. valarray<_Ty> x(n);
  612.     _Ty r(1), d(10000), s(0);
  614.     for(int m=0; m<10000; m++)
  615.     {
  616. for(i=0; i<n; i++)
  617. x[i] = a[i] + (b[i] - a[i]) * rand_01_One(r);
  619.         s = s + FunctionValueMC2D(n, x) / d;
  620.     }
  621.     
  622. for(i=0; i<n; i++) s = s * (b[i] - a[i]);
  623.     
  624. return(s);
  625. }
  627. #endif // _INTEGRAL_INL