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
Check if all input arguments are either scalar or of common size. | |
compute abscissas and weight factors for Gaussian quadratures | |
Display message and abort function. | |
Check if variables or functions are defined. | |
Construct a string from a function handle. | |
Convert integer to string (Fast version). | |
True if object is a given class. | |
Permute array dimensions. | |
Remove singleton dimensions. |
returns the probability for rectangular regions. | |
returns the probability for rectangular regions. |
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
Comments or corrections to the WAFO group