Home > wafo > misc > qrule.m

# qrule

## PURPOSE

compute abscissas and weight factors for Gaussian quadratures

## SYNOPSIS

[bp,wf]=qrule(n,wfun,alpha,beta)

## DESCRIPTION

 QRULE compute abscissas and weight factors for Gaussian quadratures

CALL:  [bp,wf]=qrule(n,wfun,alpha,beta)

bp = base points (abscissas)
wf = weight factors
n  = number of base points (abscissas) (integrates a (2n-1)th order
polynomial exactly)
wfun = weight function%
1  p(x)=1                       a =-1,   b = 1 Legendre (default)
2  p(x)=1/sqrt((x-a)*(b-x)),    a =-1,   b = 1 Chebyshev of the
first kind
3  p(x)=sqrt((x-a)*(b-x)),      a =-1,   b = 1 Chebyshev of the
second kind
4  p(x)=sqrt((x-a)/(b-x)),      a = 0,   b = 1
5  p(x)=1/sqrt(b-x),            a = 0,   b = 1
6  p(x)=sqrt(b-x),              a = 0,   b = 1
7  p(x)=(x-a)^alpha*(b-x)^beta  a =-1,   b = 1 Jacobi
alpha, beta >-1 (default alpha=beta=0)
8  p(x)=x^alpha*exp(-x)         a = 0,   b = inf generalized Laguerre
9  p(x)=exp(-x^2)               a =-inf, b = inf Hermite
10  p(x)=1                       a =-1,   b = 1 Legendre (slower than 1)

The Gaussian Quadrature integrates a (2n-1)th order
polynomial exactly and the integral is of the form
b                         n
Int ( p(x)* F(x) ) dx  =  Sum ( wf_j* F( bp_j ) )
a                        j=1
See also  gaussq

## CROSS-REFERENCE INFORMATION

This function calls:
 eig Eigenvalues and eigenvectors. error Display message and abort function. gamma Gamma function. grule computes base points and weight factors for a Gauss-
This function is called by:
 qrule2d compute abscissas and weight factors for Gaussian quadratures

## SOURCE CODE

001 function [bp,wf]=qrule(n,wfun,alpha,beta)
002 %QRULE compute abscissas and weight factors for Gaussian quadratures
003 %
004 %CALL:  [bp,wf]=qrule(n,wfun,alpha,beta)
005 %
006 %  bp = base points (abscissas)
007 %  wf = weight factors
008 %  n  = number of base points (abscissas) (integrates a (2n-1)th order
009 %       polynomial exactly)
010 %wfun = weight function%
011 %     1  p(x)=1                       a =-1,   b = 1 Legendre (default)
012 %     2  p(x)=1/sqrt((x-a)*(b-x)),    a =-1,   b = 1 Chebyshev of the
013 %                                                             first kind
014 %     3  p(x)=sqrt((x-a)*(b-x)),      a =-1,   b = 1 Chebyshev of the
015 %                                                            second kind
016 %     4  p(x)=sqrt((x-a)/(b-x)),      a = 0,   b = 1
017 %     5  p(x)=1/sqrt(b-x),            a = 0,   b = 1
018 %     6  p(x)=sqrt(b-x),              a = 0,   b = 1
019 %     7  p(x)=(x-a)^alpha*(b-x)^beta  a =-1,   b = 1 Jacobi
020 %                                     alpha, beta >-1 (default alpha=beta=0)
021 %     8  p(x)=x^alpha*exp(-x)         a = 0,   b = inf generalized Laguerre
022 %     9  p(x)=exp(-x^2)               a =-inf, b = inf Hermite
023 %    10  p(x)=1                       a =-1,   b = 1 Legendre (slower than 1)
024 %
025 %  The Gaussian Quadrature integrates a (2n-1)th order
026 %  polynomial exactly and the integral is of the form
027 %           b                         n
028 %          Int ( p(x)* F(x) ) dx  =  Sum ( wf_j* F( bp_j ) )
029 %           a                        j=1
031
032 % Reference
033 %   wfun 1: copied from grule.m in NIT toolbox, see ref [2]
034 %   wfun 2-6: see ref [4]
035 %   wfun 7-10:  Adapted from Netlib routine gaussq.f see ref [1,3]
036 %
037 % [1]  Golub, G. H. and Welsch, J. H. (1969)
038 % 'Calculation of Gaussian Quadrature Rules'
039 %  Mathematics of Computation, vol 23,page 221-230,
040 %
041 % [2] Davis and Rabinowitz (1975) 'Methods of Numerical Integration', page 365,
043 %
044 % [3]. Stroud and Secrest (1966), 'gaussian quadrature formulas',
045 %      prentice-hall, Englewood cliffs, n.j.
046 %
047 % [4] Abromowitz and Stegun (1954) ''
048
049 %  By Bryce Gardner, Purdue University, Spring 1993.
050 % Modified by Per A. Brodtkorb 19.02.99 pab@marin.ntnu.no
051 % to compute other quadratures  than the default
052 if nargin<4|isempty(beta),
053  beta=0;
054 end
055
056 if nargin<3|isempty(alpha),
057   alpha=0;
058 end
059 if alpha<=-1 | beta <=-1,
060   error('alpha and beta must be greater than -1')
061 end
062
063 if nargin<2|isempty(wfun),
064   wfun=1;
065 end
066
067
068 switch wfun, %
069   case 1,
070     %  This routine computes Gauss base points and weight factors
071     %  using the algorithm given by Davis and Rabinowitz in 'Methods
072     %  of Numerical Integration', page 365, Academic Press, 1975.
073     bp=zeros(n,1); wf=bp; iter=2; m=fix((n+1)/2); e1=n*(n+1);
074     mm=4*m-1; t=(pi/(4*n+2))*(3:4:mm); nn=(1-(1-1/n)/(8*n*n));
075     xo=nn*cos(t);
076     for j=1:iter
077       pkm1=1; pk=xo;
078       for k=2:n
079     t1=xo.*pk; pkp1=t1-pkm1-(t1-pkm1)/k+t1;
080     pkm1=pk; pk=pkp1;
081       end
082       den=1.-xo.*xo; d1=n*(pkm1-xo.*pk); dpn=d1./den;
083       d2pn=(2.*xo.*dpn-e1.*pk)./den;
084       d3pn=(4*xo.*d2pn+(2-e1).*dpn)./den;
085       d4pn=(6*xo.*d3pn+(6-e1).*d2pn)./den;
086       u=pk./dpn; v=d2pn./dpn;
087       h=-u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
088       p=pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
089       dp=dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
090       h=h-p./dp; xo=xo+h;
091     end
092     bp=-xo-h;
093     fx=d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(...
094     d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
095     wf=2*(1-bp.^2)./(fx.*fx);
096     if (m+m) > n, bp(m)=0; end
097     if ~((m+m) == n), m=m-1; end
098     jj=1:m; n1j=(n+1-jj); bp(n1j)=-bp(jj); wf(n1j)=wf(jj);
099     % end
100
101  case 2, % p(x)=1/sqrt((x-a)*(b-x)), a=-1 and b=1 (default)
102   j=[1:n];
103   wf = ones(1,n) * pi / n;
104   bp=cos( (2*j-1)*pi / (2*n) );
105
106  case 3, %p(x)=sqrt((x-a)*(b-x)),   a=-1   and b=1
107   j=[1:n];
108   wf = pi/ (n+1) *sin( j*pi / (n+1) ).^2;
109   bp=cos( j*pi / (n+1) );
110
111  case 4, %p(x)=sqrt((x-a)/(b-x)),   a=0   and b=1
112     j=[1:n];
113     bp=cos( (2*j-1)*pi /2/ (2*n+1) ).^2;
114     wf=2*pi.*bp/(2*n+1) ;
115
116  case 5, % %p(x)=1/sqrt(b-x),   a=0   and b=1
117    [bp wf]=grule(2*n);
118   wf(bp<0)=[];
119   wf=wf*2;
120    bp(bp<0)=[];
121   bp=1-bp.^2;
122
123  case 6, % %p(x)=sqrt(b-x),   a=0   and b=1
124    [bp wf]=grule(2*n+1);
125   wf(bp<=0)=[];
126    bp(bp<=0)=[];
127   wf=2*bp.^2.*wf;
128   bp=1-bp.^2;
129
130  case {7,8,9,10} ,%
131   %7 p(x)=(x-a)^alpha*(b-x)^beta a=-1 b=1 Jacobi
132   %8 p(x)=x^alpha*exp(-x) a=0,   b=inf generalized Laguerre
133   %9 p(x)=exp(-x^2)       a=-inf, b=inf Hermite
134   %10 p(x)=1               a=-1 b=1        Legendre slower than 1
135   % this procedure uses the coefficients a(j), b(j) of the
136   %      recurrence relation
137   %
138   %           b p (x) = (x - a ) p   (x) - b   p   (x)
139   %            j j            j   j-1       j-1 j-2
140   %
141   %      for the various classical (normalized) orthogonal polynomials,
142   %      and the zero-th moment
143   %
144   %           muzero = integral w(x) dx
145   %
146   %      of the given polynomial's weight function w(x).  since the
147   %      polynomials are orthonormalized, the tridiagonal matrix is
148   %      guaranteed to be symmetric.
149   %
150   %
151   %         the input parameter alpha is used only for laguerre and
152   %      jacobi polynomials, and the parameter beta is used only for
153   %      jacobi polynomials.  the laguerre and jacobi polynomials
154   %      require the gamma function.
155
156   a=zeros(n,1);
157   b=zeros(n-1,1);
158   switch wfun
159     case 7,  %jacobi
160       ab = alpha + beta;
161       abi = 2 + ab;
162       muzero = 2^(ab + 1) * gamma(alpha + 1) * gamma(beta + 1) / gamma(abi);
163       a(1) = (beta - alpha)/abi;
164       b(1) = sqrt(4*(1 + alpha)*(1 + beta)/((abi + 1)*abi^2));
165       a2b2 = beta^2 - alpha^2;
166
167       i = (2:n-1)';
168       abi = 2*i + ab;
169       a(i) = a2b2./((abi - 2).*abi);
170       a(n) =a2b2./((2*n - 2+ab).*(2*n+ab));
171       b(i) = sqrt (4*i.*(i + alpha).*(i + beta)*(i + ab)./((abi.^2 - 1).*abi.^2));
172
173     case 8, % Laguerre
174       muzero=gamma(alpha+1);
175       i = (1:n-1)';
176       a(i) = 2 .* i - 1 + alpha;
177       a(n)=2*n-1+alpha;
178       b = sqrt( i .* (i + alpha) );
179     case 9, %Hermite
180       i = (1:(n-1))';
181       muzero = sqrt(pi);
182       %a=zeros(m,1);
183       b=sqrt(i/2);
184     case 10,  % legendre NB! much slower than wfun=1
185       muzero = 2;
186       i = (1:n-1)';
187       abi = i;
188       b(i) = abi./sqrt(4*abi.^2 - 1);
189
190   end
191
192   %[v d] = eig( full(spdiags([b a b],-1:1,n,n )));
193   [v d ] = eig( diag(a) + diag(b,1) + diag(b,-1) );
194   wf = v(1,:);
195   if 1,
196     [bp i] = sort( diag(d) );
197     wf = wf(i);
198   else % save some valuable time by not sorting
199     bp = diag(d) ;
200   end
201   bp=bp';
202
203   wf = muzero.* wf.^2;
204
205 otherwise, error('unknown weight function')
206 end
207
208 % end
209

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003