MATLAB: Quad 2d giving inconsistant answers

quad2d

Hello,
I'm trying to use quad2d to integrate a function whose value is 0 outside of the (-5,5)(-5,5) box. However, I get different answers if I tell quad to integrate that box specifically or a larger box around it. Any idea what's happening? I'm using R2009a. (Sorry for the inefficent Matlab coding, it's partly because I hacked the code to work with quad2d)
quad2d(@(x,y)(of_Nd_ExpHimmelblau(x,y)),-5,5,-5,5,'AbsTol',10^-14,'RelTol',10^-14)
ans =
38.008891501144596
quad2d(@(x,y)(of_Nd_ExpHimmelblau(x,y)),-5,7,-5,11,'AbsTol',10^-14,'RelTol',10^-14)
ans =
38.008544462736467
function outMat=of_Nd_ExpHimmelblau(x1,y1)
[N,D]=size(x1);
x=zeros(N*D,2);
x(:,1)=x1(:);
x(:,2)=y1(:);
beta=0.01;
[numSamples,numDim]=size(x);
xmin=-5;
xmax=5;
out=zeros(numSamples,1);
for jj=1:numSamples
    if sum(x(jj,:)<xmin)>0 || sum(x(jj,:)>xmax)>0
        out(jj)=0;
    else
        for ii=1:(numDim-1)
            xx=x(jj,ii);
            yy=x(jj,ii+1);
            out(jj)=out(jj)+exp(-((xx^2+yy-11)^2+(xx+yy^2-7)^2)*beta);
        end
    end
end
outMat=zeros(N,D);
outMat(:)=out(:);

Best Answer

Did you turn the warnings off? Your second integration failed:
>> quad2d(@(x,y)(of_Nd_ExpHimmelblau(x,y)),-5,5,-5,5,'AbsTol',10^-14,'RelTol',10^-14)
Warning: RelTol was increased to 100*eps('double') = 2.22045e-14. 
> In quad2d at 173 
ans =
  38.008891501144596
>> quad2d(@(x,y)(of_Nd_ExpHimmelblau(x,y)),-5,7,-5,11,'AbsTol',10^-14,'RelTol',10^-14)
Warning: RelTol was increased to 100*eps('double') = 2.22045e-14. 
> In quad2d at 173 
Warning: Reached the maximum number of function evaluations (2000). The result fails the global error test. 
> In quad2d at 241 
ans =
  38.008544462736474
In general you should use the limits (function handles if need be) to avoid integrating over discontinuities. Setting the integrand to zero outside of the region of interest was a bad habit encouraged of necessity for DBLQUAD, but it is a horrible thing to do, sending accuracy down and run-times up by 10x. If you have sharp edges or discontinuities, especially if you require high accuracy, you really need to break the itegration up into pieces, putting the nasty bits on the boundaries. Of course in this case it is not really necessary since the integrand is zero on the other side of the discontinuity. The reason QUAD2D has more trouble with discontinuities than DBLQUAD, generally speaking, is because it refines the mesh in 2D instead of in one dimension at a time. A curvilinear discontinuity will attract an unnecessarily fine mesh along its length as well as "across" it. -- Mike
Related Question