Home > wafo > misc > gaussq.m

gaussq

PURPOSE ^

Numerically evaluates a integral using a Gauss quadrature.

SYNOPSIS ^

[int, tol1,ix]= gaussq(fun,xlow,xhigh,tol,trace,varargin)

DESCRIPTION ^

 GAUSSQ Numerically evaluates a integral using a Gauss quadrature. 
        The Quadrature integrates a (2m-1)th order  polynomial exactly  
        and the  integral is of the form 
               b                               
              Int (p(x)* Fun(x)) dx  
               a                  
  
  CALL: 
  [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],p1,p2,....) 
  [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,p1,p2,....) 
  [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,beta,p1,p2,....) 
  
        int = evaluated integral 
        tol = absolute tolerance abs(int-intold) 
        Fun = string containing the name of the function or a directly given  
              expression enclosed in parenthesis. The function may depend on 
              parameters alpha,beta, pi. 
 xlow,xhigh = integration limits 
     reltol = relative tolerance default=1e-3 
       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) 
  
    trace = for non-zero TRACE traces the function evaluations  
               with a point plot of the integrand. 
       gn = number of base points and weight points to start the  
             integration with (default=2) 
 p1,p2,...= coefficients to be passed directly to function Fun:    
                   G = Fun(x,p1,p2,...). 
  
  Note that int is the common size of xlow, xhigh and p1,p2,.... 
  Example:% a) integration of x.^2 from 0 to 2 and from 1 to 4 
          % b) integration of x^2*exp(-x) from zero to infinity  
  
     gaussq('(x.^2)',[0 1],[2 4])             % a) 
     gaussq('(1)',0,inf,[1e-3 8],[],2)        % b) 
     gaussq('(x.^2)',0,inf,[1e-3 8],[],0)     % b) 
  
  See also  qrule, gaussq2d

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [int, tol1,ix]= gaussq(fun,xlow,xhigh,tol,trace,varargin) 
002 %GAUSSQ Numerically evaluates a integral using a Gauss quadrature. 
003 %       The Quadrature integrates a (2m-1)th order  polynomial exactly  
004 %       and the  integral is of the form 
005 %              b                               
006 %             Int (p(x)* Fun(x)) dx  
007 %              a                  
008 % 
009 % CALL: 
010 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],p1,p2,....) 
011 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,p1,p2,....) 
012 % [int, tol] = gaussq('Fun',xlow,xhigh,[reltol wfun],[trace,gn],alpha,beta,p1,p2,....) 
013 % 
014 %       int = evaluated integral 
015 %       tol = absolute tolerance abs(int-intold) 
016 %       Fun = string containing the name of the function or a directly given  
017 %             expression enclosed in parenthesis. The function may depend on 
018 %             parameters alpha,beta, pi. 
019 %xlow,xhigh = integration limits 
020 %    reltol = relative tolerance default=1e-3 
021 %      wfun = weight function 
022 %         1  p(x)=1                       a =-1,   b = 1 Legendre (default) 
023 %         2  p(x)=1/sqrt((x-a)*(b-x)),    a =-1,   b = 1 Chebyshev of the 
024 %                                                             first kind 
025 %         3  p(x)=sqrt((x-a)*(b-x)),      a =-1,   b = 1 Chebyshev of the  
026 %                                                            second kind 
027 %         4  p(x)=sqrt((x-a)/(b-x)),      a = 0,   b = 1 
028 %         5  p(x)=1/sqrt(b-x),            a = 0,   b = 1 
029 %         6  p(x)=sqrt(b-x),              a = 0,   b = 1 
030 %         7  p(x)=(x-a)^alpha*(b-x)^beta  a =-1,   b = 1 Jacobi  
031 %                                     alpha, beta >-1 (default alpha=beta=0) 
032 %         8  p(x)=x^alpha*exp(-x)         a = 0,   b = inf generalized Laguerre 
033 %         9  p(x)=exp(-x^2)               a =-inf, b = inf Hermite 
034 %        10  p(x)=1                       a =-1,   b = 1 Legendre (slower than 1) 
035 % 
036 %   trace = for non-zero TRACE traces the function evaluations  
037 %              with a point plot of the integrand. 
038 %      gn = number of base points and weight points to start the  
039 %            integration with (default=2) 
040 %p1,p2,...= coefficients to be passed directly to function Fun:    
041 %                  G = Fun(x,p1,p2,...). 
042 % 
043 % Note that int is the common size of xlow, xhigh and p1,p2,.... 
044 % Example:% a) integration of x.^2 from 0 to 2 and from 1 to 4 
045 %         % b) integration of x^2*exp(-x) from zero to infinity  
046 % 
047 %    gaussq('(x.^2)',[0 1],[2 4])             % a) 
048 %    gaussq('(1)',0,inf,[1e-3 8],[],2)        % b) 
049 %    gaussq('(x.^2)',0,inf,[1e-3 8],[],0)     % b) 
050 % 
051 % See also  qrule, gaussq2d 
052  
053  
054 % References  
055 % 
056 % [1]  Golub, G. H. and Welsch, J. H. (1969) 
057 % 'Calculation of Gaussian Quadrature Rules' 
058 %  Mathematics of Computation, vol 23,page 221-230, 
059 % 
060 % [2] Davis and Rabinowitz (1975) 'Methods of Numerical Integration', page 365, 
061 %     Academic Press. 
062 % 
063 % [3]. Stroud and Secrest (1966), 'gaussian quadrature formulas',  
064 %      prentice-hall, Englewood cliffs, n.j. 
065 %  
066 % [4] Abromowitz and Stegun (1954) 'Handbook of mathematical functions' 
067  
068  
069  
070  
071 % tested on: Matlab 5.3 
072 % history: 
073 % Revised pab 22Nov2004 
074 % -Added the possibility of using a function handle.   
075 % Revised pab 09.09.2002 
076 % -added the possibility of using a inline function. 
077 % revised pab 27.03.2000 
078 %  - fixed a bug for p1,p2,... changed to varargin in input 
079 % revised pab 19.09.1999 
080 %   documentation 
081 %  by Per A. Brodtkorb 30.03.99, 17.02.99 : 
082 %   wfun 1: from grule.m in NIT toolbox, see ref [2]  
083 %   wfun 2-6: see ref [4] 
084 %   wfun 7-10:  Adapted from Netlib routine gaussq.f see ref [1,3] 
085 %  -accept multiple integrationlimits, int is the common size of xlow and xhigh 
086 %  -optimized by only computing the integrals which did not converge. 
087 %  -enabled the integration of directly given functions enclosed in  
088 %     parenthesis. Example: integration from 0 to 2 and from 2 to 4 for x is done by: 
089 %                        gaussq('(x.^2)',[0 2],[2 4]) 
090  
091 global ALPHA1 BETA1 ALPHA2 
092 wfun=1; 
093 if nargin<4| isempty(tol), 
094   tol=1e-3; 
095 elseif length(tol)==2, 
096   wfun=tol(2); 
097    tol=tol(1); 
098 end 
099  
100 P0=varargin; 
101 NP=length(P0); 
102  
103 istart      = 0;  
104 alpha1      = 0; 
105 beta1       = 0; 
106 FindWeights = 1; 
107  
108 x = []; 
109 y = []; 
110  
111 switch wfun 
112   case 7, 
113     istart=2; 
114     if ((NP>=1) & (~isempty(P0{1}))), alpha1 = P0{1};    end 
115     if ((NP>=2) & (~isempty(P0{1}))), beta1  = P0{2};    end 
116     if isempty(ALPHA1)|isempty(BETA1), 
117     elseif ALPHA1==alpha1 & BETA1==beta1, 
118       FindWeights=0; 
119     end 
120     ALPHA1=alpha1;BETA1=beta1; 
121     %remember what type of weights are saved as global constants 
122   case 8, 
123     istart=1; 
124     if ((NP>=1) & (~isempty(P0{1}))), alpha1 = P0{1};    end 
125     if isempty(ALPHA2), 
126     elseif ALPHA2==alpha1, 
127       FindWeights=0; 
128     end 
129     ALPHA2=alpha1; 
130     %remember what type of weights are saved as global constants 
131   otherwise,FindWeights=0; 
132 end 
133  
134 P0(1:istart)=[]; 
135  
136 gn=2; 
137 if nargin <5|isempty(trace) , 
138   trace = 0;  
139 elseif length(trace)==2, 
140   gn    = trace(2); 
141   trace = trace(1); 
142 end 
143  
144  
145 if prod(size(xlow))==1,% make sure the integration limits have correct size 
146   xlow=xlow(ones(size(xhigh)));; 
147 elseif prod(size(xhigh))==1, 
148   xhigh=xhigh(ones(size(xlow)));; 
149 elseif any( size(xhigh)~=size(xlow) ) 
150   error('The input must have equal size!') 
151 end 
152 [N M]=size(xlow);%remember the size of input 
153  
154 num_parameters=NP-istart; 
155 if num_parameters>0, 
156   ptxt = sprintf('p%d,',1:num_parameters); 
157   ptxt(end)=[]; % remove ',' 
158   eval(sprintf('[%s]=deal(P0{:});',ptxt)); 
159 end 
160 %Old call 
161 %for ix=1:num_parameters, 
162 %  eval(['p' int2str(ix) '=P0{ix};']); % variable # i 
163 %end 
164  
165 if (isa(fun,'char') &  any(fun=='(')), %  & any(fun=='x'), 
166   exec_string=['y=',fun ';']; %the call function is already setup 
167 else 
168   %setup string to call the function 
169   if isa(fun,'function_handle') 
170     fun = func2str(fun); 
171   end 
172   exec_string=['y=feval(fun,x']; 
173   for ix=1:num_parameters, 
174     xvar=['p' int2str(ix)]; % variable # i 
175     if eval(['isnumeric(' ,xvar,')  & length(',xvar,'(:)) ~=1' ]) , 
176     if N*M==1,      
177       eval(['[N M]=size(', xvar, ');']);  
178     elseif  eval(['N*M~=prod(size(', xvar, '))']), 
179       error('The input must have equal size!') 
180     end     
181       eval([xvar, '=' ,xvar ,'(:);']); %make sure it is a column  
182       exec_string=[exec_string,',' xvar '(k,ones(1,gn) )']; %enable integration with multiple arguments 
183     else 
184       exec_string=[exec_string,',' xvar]; 
185     end 
186   end 
187   exec_string=[exec_string,');']; 
188 end 
189  
190  
191 nk   = N*M;% # of integrals we have to compute 
192 k    = (1:nk)'; 
193 int  = zeros(nk,1); 
194 tol1 = int; 
195  
196 wtxt    = sprintf('%d_%d',gn,wfun); % number of weights and 
197                                     % weightfunction used 
198 cbtxt = sprintf('cb%s',wtxt); %base points string 
199 cwtxt = sprintf('cw%s',wtxt); %weights string 
200 eval(sprintf('global %s %s ;',cbtxt,cwtxt)); 
201 if isempty(eval(['cb',wtxt]))|FindWeights ,   
202   % calculate the weights if necessary 
203   eval(sprintf('[%s,%s]=qrule(gn,wfun,alpha1,beta1);',cbtxt,cwtxt)); 
204 end 
205  
206 %setup mapping parameters and execution strings 
207 xlow=xlow(:); 
208 jacob=(xhigh(:)-xlow(:))/2; 
209  
210 switch wfun,% this is clumsy and can written more tidy 
211   case {1 ,10}, 
212     calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 
213     int_string=['(ones(nk,1),:).*y,2).*jacob(k);']; 
214   case 2,   
215     calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 
216     int_string=['(ones(nk,1),:).*y,2);']; 
217   case 3,  
218     calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 
219     int_string=['(ones(nk,1),:).*y,2).*jacob(k).^2;']; 
220   case 4,   
221     calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 
222     int_string=['(ones(nk,1),:).*y,2).*jacob(k)*2;']; 
223   case 5,   
224     calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 
225     int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2);']; 
226   case 6,   
227     calcx_string=['(ones(nk,1),:)).*jacob(k,ones(1, gn ))*2+xlow(k,ones(1,gn ));']; 
228     int_string=['(ones(nk,1),:).*y,2).*sqrt(jacob(k)*2).^3;']; 
229   case 7,  
230     calcx_string=['(ones(nk,1),:)+1).*jacob(k,ones(1, gn ))+xlow(k,ones(1,gn ));']; 
231     int_string=['(ones(nk,1),:).*y,2).*jacob(k).^(alpha1+beta1+1);']; 
232   case {8,9} 
233      calcx_string=['(ones(nk,1),:));']; 
234      int_string=['(ones(nk,1),:).*y,2);']; 
235 otherwise error('unknown option') 
236 end 
237  
238  
239 eval(['x=(',cbtxt, calcx_string]); % calculate the x values             
240 eval(exec_string); % calculate function values  y=fun(x)                         
241 eval(['int(k)=sum(',cwtxt, int_string]);       % do the integration 
242 int_old=int; 
243  
244 if trace==1, 
245   x_trace=x(:); 
246   y_trace=y(:); 
247 end 
248  
249 % Break out of the iteration loop for three reasons: 
250 %  1) the last update is very small (compared to int  and  compared to tol) 
251 %  2) There are more than 11 iterations. This should NEVER happen.  
252  
253 converge='n'; 
254 for i=1:10, 
255   gn=gn*2;% double the # of weights 
256   wtxt = sprintf('%d_%d',gn,wfun); % # of weights and weight function used 
257   %eval(['global cb',wtxt,' cw',wtxt]); 
258   cbtxt = sprintf('cb%s',wtxt); %base points string 
259   cwtxt = sprintf('cw%s',wtxt); %weights string 
260   eval(sprintf('global %s %s ;',cbtxt,cwtxt)); 
261   if isempty(eval(cbtxt))|FindWeights ,   
262     % calculate the weights if necessary 
263      eval(sprintf('[%s,%s]=qrule(gn,wfun,alpha1,beta1);',cbtxt,cwtxt)); 
264   end 
265   
266   eval(['x=(',cbtxt, calcx_string]); % calculate the x values   
267   eval(exec_string);                 % calculate function values  y=fun(x) 
268   eval(['int(k)=sum(',cwtxt, int_string]);% do the integration 
269    
270   if trace==1, 
271     x_trace=[x_trace;x(:)]; 
272     y_trace=[y_trace;y(:)]; 
273   end 
274  
275   tol1(k) = abs(int_old(k)-int(k)); %absolute tolerance 
276   k       = find(tol1 > abs(tol*int)); %| tol1 > abs(tol));%indices to integrals which did not converge 
277     
278   if any(k),% compute integrals again 
279       nk=length(k);%# of integrals we have to compute again 
280   else 
281     converge='y'; 
282     break; 
283   end 
284   int_old=int; 
285 end 
286  
287 int=reshape(int,N,M); % make sure int is the same size as the integration  limits 
288 if nargout>1, 
289   tol1=reshape(tol1,N,M); 
290 end 
291  
292 if converge=='n' 
293   if nk>1 
294     if (nk==N*M), 
295       disp('All integrals did not converge--singularities likely') 
296     else 
297       disp(sprintf('%d integrals did not converge--singularities likely', ... 
298            nk)) 
299     end 
300   else 
301     disp('Integral did not converge--singularity likely') 
302   end 
303 end 
304  
305 if trace==1, 
306   plot(x_trace,y_trace,'+') 
307 end 
308  
309  
310

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