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:
 error Display message and abort function. func2str Construct a string from a function handle. int2str Convert integer to string (Fast version). isa True if object is a given class. plot Linear plot. sprintf Write formatted data to string.
This function is called by:
 b04cdf Brodtkorb (2004) joint (Scf,Hd) CDF of laboratory data. b04jcdf Brodtkorb (2004) joint (Scf,Hd) CDF from Japan Sea. bmr00cdf Brodtkorb et.al. (2000) joint (Scf,Hd) CDF from North Sea. dist2dcdf Joint 2D CDF computed as int F(X1 jhscdf Joint (Scf,Hd) CDF for linear waves with a JONSWAP spectrum. jhsnlcdf Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. jhvcdf Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. jhvnlcdf Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum. mdist2dstat Mean and variance for the MDIST2D distribution. mk87cdf Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF. ohhscdf Joint (Scf,Hd) CDF for linear waves with Ochi-Hubble spectra. ohhsscdf Joint (Scf,Hd) CDF for linear waves in space with Ochi-Hubble spectra. tay81cdf Tayfun (1981) CDF of breaking limited wave heights tay81pdf Tayfun (1981) PDF of breaking limited wave heights tay90cdf Tayfun (1990) CDF of large wave heights tay90pdf Tayfun (1990) PDF of large wave heights thscdf Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. thsnlcdf Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. thsscdf Joint (Scf,Hd) CDF for linear waves in space with Torsethaugen spectra. thvcdf Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra. weib2dcdf Joint 2D Weibull cumulative distribution function weib2dstat Mean and variance for the 2D Weibull distribution

## 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 %
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,
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