# Integral not being computed correctly

Hello I'm trying to compute this triple integral

assume(H>0);
integrate(1,z,sqrt(x^2+y^2),H+sqrt(3)*y/3);
integrate(%,y,-sqrt(H^2*(3-sqrt(3))^2/4-x^2),sqrt(H^2*(3-sqrt(3))^2/4-x^2));
integrate(%,x,-H*(3-sqrt(3))/2,H*(3-sqrt(3))/2);


And there is a variable H. So I want the final result to be a function of H. But I get a different answer from sage and Mathematica. I hope the integrals are the same. What am I doing wrong? The sage output produces complex solutions because of the argument of inverse sine being larger than 1 while the Mathematica output is nice real solution. Which one is correct?

This is the sage output

This is the Mathematica output

edit retag close merge delete

to be precise, it's not Sage that computes it there, it's Maxima. (these %i and %o prompts are from Maxima).

( 2014-12-12 05:47:40 -0600 )edit

Sort by » oldest newest most voted

The complex pieces cancel, I believe.

sage: f=integrate(integrate(integrate(1,z,sqrt(x^2+y^2),H+sqrt(3)*y/3),y,-sqrt(H^2*(3-sqrt(3))^2/4-x^2),sqrt(H^2*(3-sqrt(3))^2/4-x^2)),x,-H*(3-sqrt(3))/2,H*(3-sqrt(3))/2)
sage: f.subs(H=1)
1/4*sqrt(2)*(((3*sqrt(3)) - 6)*sqrt(-3*sqrt(3) + 6)*arcsin(1/2*sqrt(2)*sqrt(-3*sqrt(3) + 6)*(5*sqrt(3) - 9)/(12*sqrt(3) - 21)) + I*((3*sqrt(3)) - 6)*sqrt(sqrt(3) - 2)*arcsin(1/2*sqrt(2)*sqrt(-3*sqrt(3) + 6)*(5*sqrt(3) - 9)/(12*sqrt(3) - 21)) - ((3*sqrt(3))*(2*sqrt(2)) - 12*sqrt(2))*arcsin(1/2*sqrt(2)*sqrt(-3*sqrt(3) + 6)*(5*sqrt(3) - 9)/(12*sqrt(3) - 21)))
sage: f.subs(H=1).n()
1.09351366847695 - 1.41499750148729e-17*I
sage: f.subs(H=2).n()
8.74810934781556 - 1.13199800118983e-16*I


As you can see, the problem is that Sage gets something that isn't simplifying enough to get rid of vanishingly small imaginary parts when numerically evaluated, though

sage: plot(lambda x: f.subs(H=x).n().real(),(x,0,5))


probably gives what you want.

However, in general simplification is a computationally hard problem. You might have to 'unsimplify' before a good simplification is possible.

sage: f.expand().canonicalize_radical()
-3/4*(sqrt(3)*sqrt(2) - sqrt(2))*H^3*sqrt(sqrt(3) - 2)*arcsinh((3*sqrt(3) - 5)*sqrt(sqrt(3) - 2)/(4*sqrt(3)*sqrt(2) - 7*sqrt(2))) - 1/4*H^3*(-12*I*sqrt(3) + 24*I)*arcsinh((3*sqrt(3) - 5)*sqrt(sqrt(3) - 2)/(4*sqrt(3)*sqrt(2) - 7*sqrt(2)))


Unfortunately, due to the 'canonical' nature of this simplification, this gives quite different answers, presumably wrong, which the documentation for canonicalize_radical points out.

This all doesn't answer the question of what the right answer is, but I can't see your picture in any case, and hopefully it does help some.

more

I have corrected the image hosting... And they both have H^3 term but the constants in front of it are different...

( 2014-12-12 04:04:09 -0600 )edit

how do you know that Mathematica's answer is correct?

( 2014-12-12 06:32:41 -0600 )edit

you can certainly compare the numeric outputs for constants: in maxima you do

 assume(H>0);
integrate(1,z,sqrt(x^2+y^2),H+sqrt(3)*y/3);
integrate(%,y,-sqrt(H^2*(3-sqrt(3))^2/4-x^2),sqrt(H^2*(3-sqrt(3))^2/4-x^2));
r:integrate(%,x,-H*(3-sqrt(3))/2,H*(3-sqrt(3))/2);
expand(ev(subst([H=1],r),numer));


and you get 0.729009112317963 - 13.831305954985441e-9 %i, while the Mathematica's constant term is

 ev(-(-3+sqrt(3))*%pi/4,numer);


is 0.9958449670166816. So it's really different answers, and probably Maxima is wrong, as the imaginary term looks suspiciously large.

The only way to find out for sure is to compute this by some other means...

more