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
Eigenvalues and eigenvectors. | |
Display message and abort function. | |
Gamma function. | |
computes base points and weight factors for a Gauss- |
compute abscissas and weight factors for Gaussian quadratures |
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 030 % See also gaussq 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, 042 % Academic Press. 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
Comments or corrections to the WAFO group