Numeric evaluation of integrals



Shaun Jackman <sjackman@telus.net> writes:

> I'm trying to evaluate the surface area of a hemisphere within the
> rectangle (-1..1,-1..1). I've managed to get Maxima to give me a
> negative number for the surface area, which is quite confusing. When
> I do something similar in Maple I get a number near 4.22. What have
> I done to confuse Maxima?

> z(x,y) := sqrt(4-x^2-y^2);
> zx(x,y) := diff( z(x,y), x);
> zy(x,y) := diff( z(x,y), y);
> N(x,y) := sqrt( 1 + zx(x,y)^2 + zy(x,y)^2);
> integrate( integrate( N(x,y), x,-1,1), y,-1,1);
> negative
> negative
> ev(%,limit);
> ev(%,integrate);
> negative
> ev(%,numer);
> -8.156085826898927

Here is a variation on this theme which will work, but any slight
change seems to expose yet another known bug or shortcoming in Maxima.
Also, things evolve in a completely different way if you RATSIMP your
integrand first.

dS:(sqrt(4-x^2-y^2),sqrt(diff(%%,x)^2+diff(%%,y)^2+1));
integrate(%,x);
realpart(%);
ev(%,x=1)-ev(%,x=-1);
integrate(%,y,-1,1);
/* Check the result */
%,numer;
4*romberg(romberg(''dS,x,0,1),y,0,1);
load("qq");
4*quanc8(quanc8(''dS,x,0,1),y,0,1);
/* Actually, your result is correct mod %pi */
(-8.156085826898927-%)/%pi,numer;

Actually, in the case of your double integral, internally Maxima does
much the same thing as the snippet above shows.  In particular, it
correctly computes a primitive function with respect to x in terms of
ATAN2, but then expresses it by ATAN instead, and at this point gets
the inner integral wrong by 2*%pi (and hence the double integral wrong
by 4*%pi).  This happens in INTSUBS in src/defint.lisp.

Wolfgang