Re: a bug in Integrate (2nd message)
- To: mathgroup at smc.vnet.net
- Subject: [mg122308] Re: a bug in Integrate (2nd message)
- From: dimitris <dimmechan at yahoo.com>
- Date: Mon, 24 Oct 2011 05:15:09 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201110221007.GAA29655@smc.vnet.net> <j80q6r$ad0$1@smc.vnet.net>
On Oct 23, 1:27 pm, Andrzej Kozlowski <a... at mimuw.edu.pl> wrote:
> I don' think the jump point is quite at the point where you suggest it
> is (and I don't think that your claim that Mathematica "fails to
> evaluate" these limits is correct. You can actually see that the jump
> can't be quite at the point you suggest by evaluating:
>
> N[intMath /. x -> jump - 10^(-2), 10] // Chop
>
> 0.3647401105
>
> Chop[N[intMath /. x -> jump - 10^(-3), 10]]
>
> -1.84923250024426519561423619339918941442`10.014331869778037
>
> where jump is the value you computed. As you can see, the actual jump
> takes place slightly to the left of your "jump".
>
> So where does the jump take place? Well, the function ArcTan has branch
> cuts on the imaginary axis running from I to I Infinity and -I to - I
> Infinity. Hence the jump will occur when the argument of ArcTan becomes
> purely Imaginary. One can see this for one of your ArcTan parts by
> looking at the graphs of the real and imaginary parts:
>
> Plot[{Re[(3 (x^2 - x + 1) (x^2 - x +
> I Sqrt[3]))/(-(2 Sqrt[3]
> x^4) + (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] + Sqr=
t[3] -
> 3 I) x^2 + (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] -
> 3 Sqrt[3] + 3 I) x -
> Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] +
> 2 (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] + 3 I) x^3 + Sqrt[
> 3] + 3 I)],
> Im[(3 (x^2 - x + 1) (x^2 - x +
> I Sqrt[3]))/(-(2 Sqrt[3]
> x^4) + (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] + Sqr=
t[3] -
> 3 I) x^2 + (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] -
> 3 Sqrt[3] + 3 I) x -
> Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] +
> 2 (Sqrt[3 - 3 I Sqrt[3]] Sqrt[x^2 - x + 1] + 3 I) x^3 + Sqrt[
> 3] + 3 I)]}, {x, -1, 0}]
>
> One can see from the graph that this point is around x = -0.7248,
> somewhat to the left of the value -0.7224 you computed. In principle one
> could try to compute it by solving the equation f==0, where f is the
> first of the plotted functions, but I don't think any Mathematica solver
> will return this value in a reasonable time.
>
> Andrzej Kozlowski
>
> On 22 Oct 2011, at 12:07, dimitris wrote:
>
>
>
>
>
>
>
> > int = 1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]);
>
> > As I told in the previous message Mathematica gives wrong result for
> > the following integral
>
> > Integrate[int, {x, -1, 0}]
>
> > the problem is (I guess from what I figured out in the previous
> > message) that after applying the Risch algorithm and gets an
> > antiderivative
>
> > intMath = Integrate[int, x]
>
> > Mathematica uses the fundamental theorem of Calculus incorrectly
> > failing to take into account the jump singularity that possesses the
> > antiderivative in the integration range.
>
> > Let's figure out where is this jump.
>
> > In[43]:= ff = Expand[intMath]
>
> > Out[43]:=(-(1/2))*Sqrt[(1/3)*(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]))] - (1/2)*
> > Sqrt[(1/3)*(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]))] + (1/4)*I*
> > Sqrt[(1/3)*(1 - I*Sqrt[3])]*
> > Log[16*(1 + x + x^2)^2] - (1/4)*I*Sqrt[(1/3)*(1 + I*Sqrt[3])]*
> > Log[16*(1 + x + x^2)^2] +
> > (1/4)*I*Sqrt[(1/3)*(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]))] -
> > (1/4)*I*Sqrt[(1/3)*(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]))]
>
> > In[44]:= Length[ff]
>
> > Out[44]= 6
>
> > In[45]:= ({#1, Plot[Evaluate[Re[ff[[#1]]]], {x, -3, 3}]} & ) /@
> > Range[Length[ff]]
>
> > The first two terms of ff (which contain ArcTan) have jump in the same
> > point.
>
> > In[52]:= ff[[1]]
> > ff[[2]]
>
> > Out[52]= (-(1/2))*Sqrt[(1/3)*(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]))=
]
>
> > Out[53]= (-(1/2))*Sqrt[(1/3)*(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]))]
>
> > In particular the jump is at the point:
>
> > In[100]:= (N[#1, 20] & )[Re[Select[
> > x /. Solve[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]) == 0],=
-1 <
> > Re[#1] < 0 && Abs[Abs[#1] - 0.7] < 0.1 & ][[
> > 1]]]]
>
> > Out[100]= -0.72241927849008646432307690472106948257`20.
>
> > or in symbolic form
>
> > jump = 1/6 + Re[(I*(1 + I*Sqrt[3])*(1 -
> > 5*I*Sqrt[3]))/(12*(-71*I + 9*Sqrt[3] +
> > Sqrt[-5022 - 918*I*Sqrt[3]])^(1/3)) -
> > (1/12)*I*(1 - I*Sqrt[3])*(-71*I + 9*Sqrt[3] +
> > Sqrt[-5022 - 918*I*Sqrt[3]])^(1/3)];
>
> > Unfortunately, Mathematica fails to evaluate the side limits
>
> > Limit[intMath,x->jump,Direction->1]
> > and
> > Limit[intMath,x->jump,Direction->-1]
>
> > The package NLimit fails as well.
>
> > So I think the only solution is to help a bit Mathematica in order to
> > get an antiderivative continuous in the whole real axis.
> > The question is how?
>
> > Dimitris Anagnostou
>
> > P.S. In another forum a user gave the following antiderivative:
>
> > intUser=ArcTan[(Sqrt[2]*(1 + x))/Sqrt[1 - x + x^2]]/Sqrt[2] +
> > ArcTanh[(Sqrt[2/3]*(-1 + x))/Sqrt[1 - x + x^2]]/Sqrt[6];
>
> > which is correct in the whole complex plane
>
> > D[intUser, x] - 1/((x^2 + x + 1) Sqrt[x^2 - x + 1]) // FullSimplify
> > 0
>
> > and possesses no jump in the whole real axis.
>
> > Hence the definite integral in the range [-1,0] is
>
> > In[111]:= FullSimplify[int /. x -> 0 - int /. x -> -1]
>
> > Out[111]= (3/13)*Sqrt[3*(4 + Sqrt[3])]
>
> > In[112]:= N[%] (*check*)
>
> > Out[112]= 0.956959
I finally succeeded in obtaining a continuous antiderivative in the
whole real axis using Euler's substitution:
http://planetmath.org/encyclopedia/EulersSubstitutionsForIntegration.html
The integrand
In[116]:= int = 1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]);
The substitution
In[134]:= sub = x /. Solve[Sqrt[1 - x + x^2] == x + t, x][[1]]
Out[134]= (1 - t^2)/(1 + 2*t)
In[135]:= PowerExpand[Simplify[int*dx /. {x -> sub, dx -> D[sub, t]}]]
Integrate[%, t]
antiDer = % /. t -> Sqrt[1 - x + x^2] - x (*returning to the original
variable*)
Simplify[D[%, x] - int] (*check*)
Out[135]= -((2*(1 + 2*t))/(3 + 6*t + t^2 - 2*t^3 + t^4))
Out[136]= -RootSum[
3 + 6*#1 + #1^2 -
2*#1^3 + #1^4 & , (Log[t - #1] + 2*Log[t - #1]*#1)/(3 + #1 -
3*#1^2 + 2*#1^3) & ]
Out[137]= -RootSum[
3 + 6*#1 + #1^2 -
2*#1^3 + #1^4 & , (Log[-x + Sqrt[1 - x + x^2] - #1] +
2*Log[-x + Sqrt[1 - x + x^2] - #1]*#1)/
(3 + #1 - 3*#1^2 + 2*#1^3) & ]
Out[138]= 0
Hence an antiderivative of 1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]) is
In[139]:= antiDer
Out[139]= -RootSum[3 + 6*#1 + #1^2 - 2*#1^3 + #1^4 & , (Log[-x +
Sqrt[1 - x + x^2] - #1] +
2*Log[-x + Sqrt[1 - x + x^2] - #1]*#1)/(3 + #1 - 3*#1^2 +
2*#1^3) & ]
which is continuous in the whole real axis
In[140]:=Plot[antiDer, {x, -3, 3}]
and can be used for definite integration. E.g.
In[145]:= NIntegrate[1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]), {x, -
Infinity, Infinity}]
Out[145]= 2.2869
In[149]:= Chop[N[Limit[antiDer, x -> Infinity] - Limit[antiDer, x -> -
Infinity]]]
Out[149]= 2.2869
Dimitris
- References:
- a bug in Integrate (2nd message)
- From: dimitris <dimmechan@yahoo.com>
- a bug in Integrate (2nd message)