gauleg.cpp 1.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950
  1. /****************GAUSS-LEGENDRE QUADRATURE METHOD*****************/
  2. /*****Source - "NUMERICAL RECIPES IN C", www.nr.com*****/
  3. //#include "stdafx.h"
  4. #include <math.h>
  5. #include <stdlib.h>
  6. #include "vector.h"
  7. #include "gauleg.h"
  8. #define EPS 1.0e-12
  9. /*****Abscissas and weights for the Gauss-Legendre quadrature*****/
  10. void gauleg( double x1, double x2, Vector& p, Vector& w )
  11. {
  12. int m, i, j;
  13. int const n = p.dim();
  14. double z1, z, xm, xl, pp, p3, p2, p1;
  15. m = (n+1)/2;
  16. xm = (x2+x1)/2;
  17. xl = (x2-x1)/2;
  18. double atn = 4.0 * atan(1.0);
  19. for(i=1; i<=m; i++)
  20. {
  21. z = cos(atn*(i-0.25)/(n+0.5));
  22. do {
  23. p1 = 1.0;
  24. p2 = 0.0;
  25. for(j=1; j<=n; j++)
  26. {
  27. p3 = p2;
  28. p2 = p1;
  29. p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
  30. }
  31. pp = n*(z*p1-p2)/(z*z-1.0);
  32. z1 = z;
  33. z = z1 - p1/pp;
  34. } while(fabs(z-z1) > EPS);
  35. p[i-1] = xm-xl*z;
  36. p[n-i] = xm+xl*z;
  37. w[i-1] = 2.0*xl/((1.0-z*z)*pp*pp);
  38. w[n-i] = w[i-1];
  39. }
  40. if(double(n)/2>floor(double (n)/2)) p[(n-1)/2] = (x2+x1)/2;
  41. }