gauss_jac.m

Integraton routines in matlab

function [x, w] = gauss_jac(N, alpha, beta)
%
%   compute the  nodes and weights for the 
%   Gauss-quadrature.
%
%   The nodes are the roots of the orthogonal
%   Gauss-Jacoby polynomial
%

  
  gamma2     = zeros(N,1);
  gamma2(2)  = 4*(alpha+1)*(beta+1) / ((alpha+beta+3)*(alpha+beta+2)^2);
  nn         = (3:N)
  gamma2(nn) = 4*(nn-1).*(nn+alpha-1).*(nn+beta-1).*(nn+alpha+beta-1) ...
                ./ ( (2*nn+alpha+beta-1).*(2*nn+alpha+beta-2).^2 ...
                     .*(2*nn+alpha+beta-3) );
         
  gamma = sqrt(gamma2)
  
  if( alpha == beta | alpha == -beta ) 
     delta = zeros(N,1);
  else
     nn     = (1:N)
     delta  = (alpha+beta)*(alpha-beta) ...
              ./((2*nn+alpha+beta).*(2*nn+alpha+beta-2));
  end
         
  T = diag(delta,0) + diag(gamma(2:N),1) + diag(gamma(2:N),-1);
  [V,D] = eig(T);
  x = diag(D); [x,i] = sort(x);
  
  
  w = zeros(N);


end