gauleg_dke_yp

    Calculation of the weights for the Gauss-Legendre quadrature. Ribicki's algorithm 

    Input:
    
        - n: number of weigth [1,1]
          (default = 50)
        - epsi: precision [1,1]
          (default = 1e-14)
        - x1: lower limit [1,1]
          (default = -1)
        - x2: upper limit [1,1]
          (default = +1)

    Output:

        - x: Gauss-Legendre quadrature x-values [n,1]
        - w: Gauss-Legendre quadrature weights at x-values [n,1]


by Y.PEYSSON CEA-DRFC 16/05/1991 <yves.peysson@cea.fr>
revised for MatLab 4.1 (30/08/1994)

0001 function [x,w] = gauleg_dke_yp(n,epsi,x1,x2)
0002 %
0003 %    Calculation of the weights for the Gauss-Legendre quadrature. Ribicki's algorithm
0004 %
0005 %    Input:
0006 %
0007 %        - n: number of weigth [1,1]
0008 %          (default = 50)
0009 %        - epsi: precision [1,1]
0010 %          (default = 1e-14)
0011 %        - x1: lower limit [1,1]
0012 %          (default = -1)
0013 %        - x2: upper limit [1,1]
0014 %          (default = +1)
0015 %
0016 %    Output:
0017 %
0018 %        - x: Gauss-Legendre quadrature x-values [n,1]
0019 %        - w: Gauss-Legendre quadrature weights at x-values [n,1]
0020 %
0021 %
0022 %by Y.PEYSSON CEA-DRFC 16/05/1991 <yves.peysson@cea.fr>
0023 %revised for MatLab 4.1 (30/08/1994)
0024 %
0025 if nargin < 0,
0026     infoyp(2,'Wrong number of input arguments for gaulegyp');
0027     return;
0028 end
0029 if nargin == 0,
0030     n = 50;
0031     epsi = 1e-14;
0032     x1 = -1;
0033     x2 = 1;
0034 elseif nargin == 1,
0035     epsi = 1e-14;
0036     x1 = -1;
0037     x2 = 1;
0038 elseif nargin == 2,
0039     x1 = -1;
0040     x2 = 1;
0041 elseif nargin == 3,
0042     x2 = 1;
0043 end
0044 %
0045 m = (n+1)/2;
0046 xm = 0.5*(x2+x1);xl = 0.5*(x2-x1);
0047 for i=1:m,
0048      z = cos(pi*(i-0.25)/(n+0.5));
0049      p1 = 1.0;
0050      p2 = 0.0;
0051      for j=1:n,
0052           p3 = p2;p2 = p1;p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
0053      end
0054      pp = n*(z*p1-p2)/(z*z-1.0);z1 = z;z = z1-p1/pp;
0055      test = abs(z-z1)-epsi;
0056      while test > 0,
0057           p1 = 1.0;p2 = 0.0;
0058           for j=1:n
0059                p3 = p2;p2 = p1;p1 = ((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
0060           end
0061           pp = n*(z*p1-p2)/(z*z-1.0);z1 = z;z = z1-p1/pp;
0062           test = abs(z-z1)-epsi;
0063      end
0064      x(i) = xm-xl*z;x(n+1-i) = xm+xl*z;
0065      w(i) = 2.0*xl/((1.0-z*z)*pp*pp);w(n+1-i) = w(i);
0066 end
0067 %
0068 x = x(:);
0069 w = w(:);