       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

```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))/(-(2 Sqrt
>           x^4) + (Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] + Sqr=
t -
>          3 I) x^2 + (Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] -
>          3 Sqrt + 3 I) x -
>       Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] +
>       2 (Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] + 3 I) x^3 + Sqrt[
>       3] + 3 I)],
>   Im[(3 (x^2 - x + 1) (x^2 - x +
>         I Sqrt))/(-(2 Sqrt
>           x^4) + (Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] + Sqr=
t -
>          3 I) x^2 + (Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] -
>          3 Sqrt + 3 I) x -
>       Sqrt[3 - 3 I Sqrt] Sqrt[x^2 - x + 1] +
>       2 (Sqrt[3 - 3 I Sqrt] 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:= ff = Expand[intMath]
>
> > Out:=(-(1/2))*Sqrt[(1/3)*(1 + I*Sqrt)]*
> >  ArcTan[(3*(1 - x + x^2)*(I*Sqrt - x + x^2))/(3*I + Sqrt -
> >      2*Sqrt*x^4 -
> >            Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2] +
> >      2*x^3*(3*I + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]) +
> >      x*(3*I - 3*Sqrt + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]=
)
> > +
> >      x^2*(-3*I + Sqrt +
> >         Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]))] - (1/2)*
> >  Sqrt[(1/3)*(1 - I*Sqrt)]*
> >     ArcTan[(3*(I*Sqrt + x - x^2)*(1 - x + x^2))/(3*I - Sqrt +
> >      2*Sqrt*x^4 + Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2] +
> >            x^3*(6*I - 2*Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + =
x^2]) +
> >      x*(3*I + 3*Sqrt - Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2]=
)
> > -
> >      x^2*(3*I + Sqrt +
> >         Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2]))] + (1/4)*I*
> >  Sqrt[(1/3)*(1 - I*Sqrt)]*
> >     Log[16*(1 + x + x^2)^2] - (1/4)*I*Sqrt[(1/3)*(1 + I*Sqrt)]*
> >  Log[16*(1 + x + x^2)^2] +
> >   (1/4)*I*Sqrt[(1/3)*(1 + I*Sqrt)]*
> >  Log[(1 + x + x^2)*(11*I + 4*Sqrt + (11*I + 4*Sqrt)*x^2 +
> >            10*I*Sqrt[1 - I*Sqrt]*Sqrt[1 - x + x^2] -
> >      x*(17*I + 4*Sqrt +
> >         8*I*Sqrt[1 - I*Sqrt]*Sqrt[1 - x + x^2]))] -
> >   (1/4)*I*Sqrt[(1/3)*(1 - I*Sqrt)]*
> >  Log[(1 + x + x^2)*(-11*I + 4*Sqrt + (-11*I + 4*Sqrt)*x^2 -
> >            10*I*Sqrt[1 + I*Sqrt]*Sqrt[1 - x + x^2] +
> >      x*(17*I - 4*Sqrt +
> >         8*I*Sqrt[1 + I*Sqrt]*Sqrt[1 - x + x^2]))]
>
> > In:= Length[ff]
>
> > Out= 6
>
> > In:= ({#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:= ff[]
> > ff[]
>
> > Out= (-(1/2))*Sqrt[(1/3)*(1 + I*Sqrt)]*
> > ArcTan[(3*(1 - x + x^2)*(I*Sqrt - x + x^2))/
> >       (3*I + Sqrt - 2*Sqrt*x^4 -
> >     Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2] +
> >          2*x^3*(3*I + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]=
) +
> >     x*(3*I - 3*Sqrt + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]) +
>
> >     x^2*(-3*I + Sqrt + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]))=
]
>
> > Out= (-(1/2))*Sqrt[(1/3)*(1 - I*Sqrt)]*
> > ArcTan[(3*(I*Sqrt + x - x^2)*(1 - x + x^2))/
> >       (3*I - Sqrt + 2*Sqrt*x^4 +
> >     Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2] +
> >          x^3*(6*I - 2*Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2]=
) +
> >     x*(3*I + 3*Sqrt - Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2]) -
>
> >     x^2*(3*I + Sqrt + Sqrt[3 + 3*I*Sqrt]*Sqrt[1 - x + x^2]))]
>
> > In particular the jump is at the point:
>
> > In:= (N[#1, 20] & )[Re[Select[
> >    x /. Solve[3*I + Sqrt - 2*Sqrt*x^4 -
> >        Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2] +
> >        2*x^3*(3*I + Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]) +
> >        x*(3*I - 3*Sqrt +
> >           Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]) +
> >        x^2*(-3*I + Sqrt +
> >           Sqrt[3 - 3*I*Sqrt]*Sqrt[1 - x + x^2]) == 0],=
-1 <
> >       Re[#1] < 0 && Abs[Abs[#1] - 0.7] < 0.1 & ][[
> >        1]]]]
>
> > Out= -0.72241927849008646432307690472106948257`20.
>
> > or in symbolic form
>
> > jump = 1/6 + Re[(I*(1 + I*Sqrt)*(1 -
> >         5*I*Sqrt))/(12*(-71*I + 9*Sqrt +
> >          Sqrt[-5022 - 918*I*Sqrt])^(1/3)) -
> >        (1/12)*I*(1 - I*Sqrt)*(-71*I + 9*Sqrt +
> >        Sqrt[-5022 - 918*I*Sqrt])^(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*(1 + x))/Sqrt[1 - x + x^2]]/Sqrt +
> > ArcTanh[(Sqrt[2/3]*(-1 + x))/Sqrt[1 - x + x^2]]/Sqrt;
>
> > 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:= FullSimplify[int /. x -> 0 - int /. x -> -1]
>
> > Out= (3/13)*Sqrt[3*(4 + Sqrt)]
>
> > In:= N[%] (*check*)
>
> > Out= 0.956959

I finally succeeded in obtaining a continuous antiderivative in the
whole real axis using Euler's substitution:

The integrand

In:= int = 1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]);

The substitution

In:= sub = x /. Solve[Sqrt[1 - x + x^2] == x + t, x][]

Out= (1 - t^2)/(1 + 2*t)

In:= 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= -((2*(1 + 2*t))/(3 + 6*t + t^2 - 2*t^3 + t^4))

Out= -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= -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= 0

Hence an antiderivative of 1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]) is

In:= antiDer

Out= -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:=Plot[antiDer, {x, -3, 3}]

and can be used for definite integration. E.g.

In:= NIntegrate[1/((x^2 + x + 1)*Sqrt[x^2 - x + 1]), {x, -
Infinity, Infinity}]

Out= 2.2869

In:= Chop[N[Limit[antiDer, x -> Infinity] - Limit[antiDer, x -> -
Infinity]]]

Out= 2.2869

Dimitris

```

• Prev by Date: Re: a bug in Integrate (2nd message)
• Next by Date: Parameter Fitting in coupled differential equation
• Previous by thread: Re: a bug in Integrate (2nd message)
• Next by thread: Integral points on elliptic curves