Why Matlab cannot calculate simple double integrals with good accuracy?

Let's say that I want to calculate the area of a circle of radius 1 which is π. I use the following code:

f=@(x,y) double(1-x.^2-y.^2>0);
integral2(f,-2,2,-2,2)

What I get using the long format is: 3.141530523062315.

I put more precision:

integral2(f,-2,2,-2,2,'Method','iterated','AbsTol',1e-10)

Now I get:

3.141535921819190

What is happening??? Why the answer is not with the precision required? And it lasts a long time to compute it. In Mathematica the code:

N[Integrate[HeavisideTheta[1 - x^2 - y^2], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}], 100]

returns the correct result with tons of digits in an instant:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

wolfram alpha

How can I make my Matlab program to get to the same precision?

[Matt J]

Well Mathematica is doing the integral symbolically, isn't it? Of course it's going to be very accurate. And, it's fast because it's a very easy integral to do symbolically.

[Paul]

I thought it would be too. But:

>> f(x,y) = heaviside(1 - x^2 - y^2)
 
f(x, y) =
 
heaviside(- x^2 - y^2 + 1)
 
>> int(int(f(x,y),x,-2,2,'IgnoreAnalyticConstraints',true),y,-2,2,'IgnoreAnalyticConstraints',true)
 
ans =
 
pi/2
 

Did I do something wrong?

Answer 1

[Matt J]

I'm guessing that numerical accuracy might suffer because the integrand is discontinuous, and the curved boundary of the discontinuity is hard to sample accurately with the small rectangular elements used in a Cartesian coordinate integral. In polar coordinates, it seems to work quite well:

format long 

f=@(r,theta) r.*(r<=1);

integral2(f,0,2,0,2*pi,'AbsTol',1e-10)

ans = 
   3.141592653589738