Re: a bug in Integrate (2nd message)
- To: mathgroup at smc.vnet.net
- Subject: [mg122301] Re: a bug in Integrate (2nd message)
- From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
- Date: Mon, 24 Oct 2011 05:13:53 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201110221007.GAA29655@smc.vnet.net> <j80q6r$ad0$1@smc.vnet.net> <7e925e27-3825-431d-a054-b259f0c6bdf5@n18g2000vbv.googlegroups.com>
Well done. It seems that the method is sufficiently algorithmic that it
ought to be used by Mathematica for this types of integrals (instead of
the Risch algorithm). Now that you pointed it out, I remembered that I
taught (or rather mentioned) this method last semester in our first year
analysis course for computer science students - but after teaching it I
immediately completely forgot all about it (integration is not
something I think about unless I have to teach it).
Andrej
On 23 Oct 2011, at 16:27, dimitris wrote:
> 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] + Sqrt[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] + Sqrt[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)