a bug in Integrate
- To: mathgroup at smc.vnet.net
- Subject: [mg122247] a bug in Integrate
- From: dimitris <dimmechan at yahoo.com>
- Date: Sat, 22 Oct 2011 06:07:15 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
The following has a connection with the recent post of mine called Simplification. http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/b6e8fc7f70a0f20f# Let int = 1/((x^2 + x + 1) Sqrt[x^2 - x + 1]); Then Mathematica's result for the definite integral in the range [-1,0] is In[74]:= Integrate[int, {x, -1, 0}] Out[74]= -((1/(4* Sqrt[3]))*(I*(Sqrt[ 1 + I*Sqrt[3]]*(-2*I* ArcTan[(6 + 3*I*Sqrt[3])/(3*I - Sqrt[3] + 3*Sqrt[1 - I*Sqrt[3]])] + 2*ArcTanh[(3*Sqrt[3])/(3*I + Sqrt[3] - Sqrt[3 - 3*I*Sqrt[3]])] - Log[11*I + 4*Sqrt[3] + 10*I*Sqrt[1 - I*Sqrt[3]]] + Log[3*(13*I + 4*Sqrt[3] + 6*I*Sqrt[3 - 3*I*Sqrt[3]])]) + Sqrt[1 - I*Sqrt[3]]*(2*ArcTanh[(3*(2*I + Sqrt[3]))/(3*I + Sqrt[3] - 3*Sqrt[1 + I*Sqrt[3]])] + 2*ArcTanh[(3*Sqrt[3])/(3*I - Sqrt[3] + Sqrt[3 + 3*I*Sqrt[3]])] + Log[(1/16)*(-11*I + 4*Sqrt[3] - 10*I*Sqrt[1 + I*Sqrt[3]])] - Log[(3/16)*(-13*I + 4*Sqrt[3] - 6*I*Sqrt[3 + 3*I*Sqrt[3]])])))) which is incorect. Indeed: In[75]:= N[%]//Chop Out[75]= -1.294232744953466 In[76]:= NIntegrate[int, {x, -1, 0}] Out[76]= 0.9272087241257184 I think the reason is that Mathematica firsts evaluates an antiderivative: In[78]:= intMath = Integrate[int, x] Out[78]= (1/(4*Sqrt[3]))*(-2*Sqrt[1 + I*Sqrt[3]]* ArcTan[(3*(1 - x + x^2)*(I*Sqrt[3] - x + x^2))/ (3*I + Sqrt[3] - 2*Sqrt[3]*x^4 - Sqrt[3 - 3*I*Sqrt[3]]*Sqrt[1 - x + x^2] + 2*x^3*(3*I + Sqrt[3 - 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]) + x*(3*I - 3*Sqrt[3] + Sqrt[3 - 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]) + x^2*(-3*I + Sqrt[3] + Sqrt[3 - 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]))] - 2*Sqrt[1 - I*Sqrt[3]]* ArcTan[(3*(I*Sqrt[3] + x - x^2)*(1 - x + x^2))/(3*I - Sqrt[3] + 2*Sqrt[3]*x^4 + Sqrt[3 + 3*I*Sqrt[3]]*Sqrt[1 - x + x^2] + x^3*(6*I - 2*Sqrt[3 + 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]) + x*(3*I + 3*Sqrt[3] - Sqrt[3 + 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]) - x^2*(3*I + Sqrt[3] + Sqrt[3 + 3*I*Sqrt[3]]*Sqrt[1 - x + x^2]))] + I*((Sqrt[1 - I*Sqrt[3]] - Sqrt[1 + I*Sqrt[3]])* Log[16*(1 + x + x^2)^2] + Sqrt[1 + I*Sqrt[3]]* Log[(1 + x + x^2)*(11*I + 4*Sqrt[3] + (11*I + 4*Sqrt[3])*x^2 + 10*I*Sqrt[1 - I*Sqrt[3]]* Sqrt[1 - x + x^2] - x*(17*I + 4*Sqrt[3] + 8*I*Sqrt[1 - I*Sqrt[3]]*Sqrt[1 - x + x^2]))] - Sqrt[1 - I*Sqrt[3]]* Log[(1 + x + x^2)*(-11*I + 4*Sqrt[3] + (-11*I + 4*Sqrt[3])*x^2 - 10*I*Sqrt[1 + I*Sqrt[3]]* Sqrt[1 - x + x^2] + x*(17*I - 4*Sqrt[3] + 8*I*Sqrt[1 + I*Sqrt[3]]*Sqrt[1 - x + x^2]))])) and then simply does (intMath /. x -> 0) - (intMath /. x -> -1). In[84]:= N[(intMath /. x -> 0) - (intMath /. x -> -1)]//Chop Out[84]= -1.2942327449534659 However the evaluated antiderivative possesses a jump singularity at near x = -0.725 and Mathematica fails unfortunately to take it into acount. In[85]:= Plot[intMath, {x, -3, 3}] Nevertheless, there are and some good news. The antiderivative is correct in the whole complex plane. In[89]:= FullSimplify[D[intMath, x] - int] Out[89]= 0 If the integration range does not include the jump, Mathematica will give correct result. E.g. In[13]:= Integrate[int, {x, 0, Infinity}] (*after several minutes*) Out[13]= (1/(4* Sqrt[3]))*(I*(Sqrt[ 1 + I*Sqrt[3]]*(2*ArcTanh[(3*Sqrt[3])/(3*I + Sqrt[3] - Sqrt[3 - 3*I*Sqrt[3]])] + Log[11*I + 4*Sqrt[3] - 8*I*Sqrt[1 - I*Sqrt[3]]] - Log[11*I + 4*Sqrt[3] + 10*I*Sqrt[1 - I*Sqrt[3]]] - Log[(-2 - I*Sqrt[3] + 2*Sqrt[1 - I*Sqrt[3]])/(-1 + Sqrt[1 - I*Sqrt[3]])] + Log[(-2 + I*Sqrt[3] + 2*Sqrt[1 - I*Sqrt[3]])/(-1 + Sqrt[1 - I*Sqrt[3]])]) + Sqrt[1 - I*Sqrt[3]]*(2* ArcTanh[(3*Sqrt[3])/(3*I - Sqrt[3] + Sqrt[3 + 3*I*Sqrt[3]])] - Log[-11*I + 4*Sqrt[3] + 8*I*Sqrt[1 + I*Sqrt[3]]] + Log[-11*I + 4*Sqrt[3] - 10*I*Sqrt[1 + I*Sqrt[3]]] - Log[(-2 - I*Sqrt[3] + 2*Sqrt[1 + I*Sqrt[3]])/(-1 + Sqrt[1 + I*Sqrt[3]])] + Log[(-2 + I*Sqrt[3] + 2*Sqrt[1 + I*Sqrt[3]])/(-1 + Sqrt[1 + I*Sqrt[3]])]))) In[14]:= N[%] Out[14]= 0.9358813101035703 + 0.*I In[15]:= NIntegrate[int, {x, 0, Infinity}] Out[15]= 0.9358813101035696 At first glance integrand 1/((x^2 + x + 1) Sqrt[x^2 - x + 1]) may look simple. This is far from true. Dimitris