Home > wafo > misc > gaussq2d.m

# gaussq2d

## PURPOSE

Numerically evaluates a 2D integral using Gauss quadrature.

## SYNOPSIS

[int, tol1,k] = gaussq2d(fun,xlow,xhigh,ylow,yhigh,tol,p1,p2,p3,p4,p5,p6,p7,p8,p9)

## DESCRIPTION

 GAUSSQ2D Numerically evaluates a 2D integral using Gauss quadrature.

CALL:  [int tol] = gaussq2d('Fun',xlow,xhigh,ylow,yhigh)
or
[int tol] = gaussq2d('Fun',xlow,xhigh,ylow,yhigh,reltol,p1,p2,...)

int = evaluated integral
tol = absolute tolerance abs(int-intold)
Fun = Fun(x,y), string containing the name of the function or
a directly given expression enclosed in parenthesis.
The function may depend on parameters pi.
xlow,xhigh = integration limits for x
ylow,yhigh = integration limits for y (obs. integration limits can
be vectors.)
reltol = relative tolerance (default 1e-3)

Note that if there are discontinuities the region of integration
should be broken up into separate pieces.  And if there are singularities,
a more appropriate integration quadrature should be used
(such as the Gauss-Chebyshev for a specific type of singularity).

Example:%  integration for (x,y) in [0,1]x[-1,3] and [2,4]x[1,2]:

p1=2.; p2=0.5;
gaussq2d('(p2*x.^2.*y+p1)',[0 2],[1 4],[-1 1],[3 2],[],p1,p2)

See also  gaussq, qrule2d

## CROSS-REFERENCE INFORMATION

This function calls:
 comnsize Check if all input arguments are either scalar or of common size. qrule2d compute abscissas and weight factors for Gaussian quadratures error Display message and abort function. exist Check if variables or functions are defined. func2str Construct a string from a function handle. int2str Convert integer to string (Fast version). isa True if object is a given class. permute Permute array dimensions. squeeze Remove singleton dimensions.
This function is called by:
 dist2dprb returns the probability for rectangular regions. weib2dprb returns the probability for rectangular regions.

## SOURCE CODE

001 function [int, tol1,k] = gaussq2d(fun,xlow,xhigh,ylow,yhigh,tol,p1,p2,p3,p4,p5,p6,p7,p8,p9)
002 %GAUSSQ2D Numerically evaluates a 2D integral using Gauss quadrature.
003 %
004 % CALL:  [int tol] = gaussq2d('Fun',xlow,xhigh,ylow,yhigh)
005 %      or
006 %        [int tol] = gaussq2d('Fun',xlow,xhigh,ylow,yhigh,reltol,p1,p2,...)
007 %
008 %       int = evaluated integral
009 %       tol = absolute tolerance abs(int-intold)
010 %       Fun = Fun(x,y), string containing the name of the function or
011 %             a directly given expression enclosed in parenthesis.
012 %             The function may depend on parameters pi.
013 %xlow,xhigh = integration limits for x
014 %ylow,yhigh = integration limits for y (obs. integration limits can
015 %                                             be vectors.)
016 %    reltol = relative tolerance (default 1e-3)
017 %
018 % Note that if there are discontinuities the region of integration
019 % should be broken up into separate pieces.  And if there are singularities,
020 % a more appropriate integration quadrature should be used
021 % (such as the Gauss-Chebyshev for a specific type of singularity).
022 %
023 %Example:%  integration for (x,y) in [0,1]x[-1,3] and [2,4]x[1,2]:
024 %
025 %   p1=2.; p2=0.5;
026 %   gaussq2d('(p2*x.^2.*y+p1)',[0 2],[1 4],[-1 1],[3 2],[],p1,p2)
027 %
029
030 % tested on: matlab 5.2
031 % history:
032 % revised pab Nov2004
033 % -wrong input to comnsize fixed
034 % -enabled function handle as input
035 % revised pab 12Nov2003
036 % replaced call to distchk with comnsize
037 %modified by Per A. Brodtkorb 17.11.98 :
038 % -accept multiple integrations limits
039 % -optimized by saving the weights as global constants and only
040 %   computing the integrals which did not converge
041 % -enabled the integration of directly given functions enclosed in
042 %   parenthesis
043 % -adopted from NIT toolbox changed name from quad2dg to gaussq2d
044
045 global bpx2 bpy2 wfxy2;
046 fv = [];
047 if isempty(bpx2)| isempty(bpy2)| isempty(wfxy2)
048   [bpx2,bpy2,wfxy2] = qrule2d(2,2);
049 end
050
051 if exist('tol')~=1,
052   tol=1e-3;
053 elseif isempty(tol),
054   tol=1e-3;
055 end
056 [errorcode ,xlow,xhigh,ylow,yhigh]=comnsize(xlow,xhigh,ylow,yhigh);
057 if errorcode > 0
058   error('Requires non-scalar arguments to match in size.');
059 end
060
061 [N M]=size(xlow);
062 %setup mapping parameters &  make sure they all are column vectors
063 xlow=xlow(:); xhigh=xhigh(:);ylow=ylow(:);yhigh=yhigh(:);
064 nk=N*M;
065
066 %Map to x
067 qx=(xhigh-xlow)/2;
068 px=(xhigh+xlow)/2;
069 x=permute(qx(:,ones(1,2), ones(1,2) ),[ 2 3 1]).*bpx2(:,:,ones(1,nk)) ...
070     +permute(px(:,ones(1,2), ones(1,2) ),[2 3 1]);
071 %Map to y
072 qy=(yhigh-ylow)/2;
073 py=(yhigh+ylow)/2;
074 y=permute(qy(:,ones(1,2), ones(1,2) ),[ 2 3 1]).*bpy2(:,:,ones(1,nk)) ...
075     +permute(py(:,ones(1,2), ones(1,2) ),[2 3 1]);
076
077 if (isa(fun,'char') & any(fun=='(')),
078   exec_string=['fv=', fun, ';']; %the call function is already setup
079 else
080   %set up the call function
081   if isa(fun,'function_handle')
082     fun = func2str(fun);
083   end
084   exec_string=['fv=', fun, ' (x,y'];
085   num_parameters=nargin-6;
086   for i=1:num_parameters,
087     exec_string=[exec_string,',p',int2str(i)];
088   end
089   exec_string=[exec_string,');'];
090 end
091 eval(exec_string);
092 int_old = squeeze(sum(sum(wfxy2(:,:,ones(1,nk)).*fv))).*qx.*qy;
093 int=zeros(size(int_old));tol1=int;
094 k=1:nk;
095
096 % Break out of the iteration loop for three reasons:
097 %  1) the last update is very small (compared to int)
098 %  2) the last update is very small (compared to tol)
099 %  3) There are more than 8 iterations. This should NEVER happen.
100
101 converge='n';
102 for i=1:8,
103   gnum=int2str(2^(i+1));
104   eval(['global bpx',gnum,' wfxy',gnum,'  bpy',gnum,';']);
105   if isempty(eval(['bpx',gnum]))| isempty(eval(['bpy',gnum])) ,
106     eval(['[bpx',gnum,',bpy',gnum,',wfxy',gnum,']=qrule2d(',gnum,',', gnum,');']);
107   end
108   eval(['x=permute(qx(k,ones(1,',gnum,'), ones(1,',gnum,') ),[ 2 3 1]).*bpx',gnum,'(:,:,ones(1,nk)) +permute(px(k,ones(1,',gnum,'), ones(1,',gnum,') ),[2 3 1]);']);
109
110   eval(['y=permute(qy(k,ones(1,',gnum,'), ones(1,',gnum,') ),[ 2 3 1]).*bpy',gnum,'(:,:,ones(1,nk)) +permute(py(k,ones(1,',gnum,'), ones(1,',gnum,') ),[2 3 1]);']);
111     eval(exec_string);
112     eval(['int(k) = squeeze(sum(sum((wfxy',gnum,'(:,:,ones(1,nk)).*fv)))).*qx(k).*qy(k);']);
113     tol1(k)=abs(int_old(k)-int(k)); % absolute tolerance
114     k=find(tol1 > abs(tol*int)|tol1 > abs(tol));%indices to integrals which did not converge
115
116     if any(k), % compute integrals again
117       nk=length(k);%# of integrals we have to compute again
118     else
119       converge='y';
120       break;
121     end
122
123      int_old(k)=int(k);
124 end
125 int=reshape(int,[N M]); % reshape int to the same size as input arguments
126 if converge=='n',
127   disp('Integral did not converge--singularity likely')
128 end
129
130
131
132

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