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: This function is called by:

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 % 
028 % See also  gaussq, qrule2d 
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