Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • 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