our good (?!) friend RootSum
- To: mathgroup at smc.vnet.net
- Subject: [mg75022] our good (?!) friend RootSum
- From: "dimitris" <dimmechan at yahoo.com>
- Date: Sat, 14 Apr 2007 01:11:39 -0400 (EDT)
I think with the last threads it is apparent that when RootSum is
generated by Integrate there are some
possibilities that Mathematica may fail in the quest of a definite
integration. Especially (and almost mostly) if there is involved
Infinity in the integration limit(s). However this failure hasn't to
do with Integrate but rather than with Limit. That is, I believe, it
has to do with the difficulty of Mathematica to evaluate sometimes
limits of
the structure Limit[RootSum[...],x->Infinity].
For example consider the definite integral of the function
[(2*Sqrt[x])/(x^3 + 3*x^2 + 2*x + 1)
in the interval [0,Infinity). The integral is convergent as someone
can easily proved
In[791]:=
((2*Sqrt[x])/(x^3 + 3*x^2 + 2*x + 1) + O[x, #1]^3 & ) /@ {0,
Infinity}
Reduce[x^3 + 3*x^2 + 2*x + 1 == 0 && x >= 0, x, Reals]
Out[791]=
{SeriesData[x, 0, {2, 0, -4, 0, 2}, 1, 6, 2], SeriesData[x, Infinity,
{2}, 5, 6, 2]}
Out[792]=
False
Let now
In[789]:=
f = HoldForm[Integrate[(2*Sqrt[x])/(x^3 + 3*x^2 + 2*x + 1), x]]
Mathematica can easily get an antiderivative
In[796]:=
g=ReleaseHold[f]
D[g, x]
Out[796]=
2*RootSum[1 + 2*#1^2 + 3*#1^4 + #1^6 & , (Log[Sqrt[x] - #1]*#1)/(2 +
6*#1^2 + 3*#1^4) & ]
Out[797]=
(2*Sqrt[x])/(1 + 2*x + 3*x^2 + x^3)
However (but not surprising considering some recent posts) the
definite integral is not an easy
task. But it was explained this is due to the incapability of Limit
and not because of Integrate
algorithm.
In[804]:=
ReleaseHold[f /. Integrate[g_, x_] :> Integrate[g, {x, 0, Infinity}]]
Out[804]=
Integrate[(2*Sqrt[x])/(1 + 2*x + 3*x^2 + x^3), {x, 0, Infinity}]
For this integral Mathematica finds an antiderivative and then (after
verifying convergence) applies
the Newton-Leibniz formula.
It is obvious that the form of Newton Leibniz formula that Mathematica
uses for this integral fails
at Infinity as the following code shows
In[805]:=
Unprotect[Limit];
Limit[a___] := Null /; (Print[InputForm[limit[a]]]; False)
In[807]:=
ReleaseHold[f /. Integrate[g_, x_] :> Integrate[g, {x, 0, Infinity}]]
(*output is ommited due to its length*)
Of course in zero there is no such a problem.
In[800]:=
Unprotect[Limit];
Clear[Limit];
In[822]:=
(Limit[g, x -> #1] & ) /@ {0, Infinity}
Out[822]=
{2*RootSum[1 + 2*#1^2 + 3*#1^4 + #1^6 & , (Log[-#1]*#1)/(2 + 6*#1^2 +
3*#1^4) & ],
Limit[2*RootSum[1 + 2*#1^2 + 3*#1^4 + #1^6 & , (Log[Sqrt[x] -
#1]*#1)/(2 + 6*#1^2 + 3*#1^4) & ], x -> Infinity]}
Let see the antiderivative
In[823]:=
Plot[g, {x, 0, 20}]
It is obvious that g->0 as x->Infinity.
The question is of course WHY Limit can't recognise this limit
behavior!
Note that as it was pointed out in another post
In[825]:=
g + O[x, Infinity]^4
Series::nmer :...
Series::nmer :...
Out[825]=
2*RootSum[1 + 2*#1^2 + 3*#1^4 + #1^6 & , (Log[Sqrt[x] - #1]*#1)/(2 +
6*#1^2 + 3*#1^4) & ] + SeriesData[x, Infinity, {}, 4, 4, 1]
Taking Limit[g->Infinity]=0 by ourselves and applying the Newton
Leibniz formula
results in the correct result.
In[838]:=
0 - g /. x -> 0
{N[%], ReleaseHold[f /. Integrate[h_, x_] :> NIntegrate[h, {x, 0,
Infinity}]]}
Out[838]=
-2*RootSum[1 + 2*#1^2 + 3*#1^4 + #1^6 & , (Log[-#1]*#1)/(2 + 6*#1^2 +
3*#1^4) & ]
Out[839]=
{0.8689140125914716+0.*I, 0.8689140125951508}
As a comparison, in another CAS (due to the policy of the forum I
cannot
mention its name) I simply took
int((2*sqrt(x))/(x^3 + 3*x^2 + 2*x + 1),x=0..infinity);
/ ----- \
| \ 50784 5 2760 3 22 |
-4 | ) _R ln(------ _R + ---- _R - -- _R)|
| / 19 19 19 |
| ----- |
\_R = %1 /
6 4 2
%1 := RootOf(33856 _Z + 1104 _Z + 32 _Z + 1)
evalf(%);
0.8689140124 - 0. I
Very good performance indeed.
In the same system I got
int((2*sqrt(x))/(x^3 + 3*x^2 + 2*x + 1),x);
/ ----- 5 3 \
| \ 1/2 50784 _R 2760 _R 22 _R |
4 | ) _R ln(x - --------- + -------- - -----)|
| / 19 19 19 |
| ----- |
\_R = %1 /
6 4 2
%1 := RootOf(33856 _Z + 1104 _Z + 32 _Z + 1)
simplify(diff(%,x));
1/2
2 x
-------------------
3 2
x + 3 x + 2 x + 1
In the field of integration I do consider Mathematica by far more
superior than other CASs
but I would like if someone could explain me why, even though both
Mathematica and the other CAS
I use, generating a quite similar root sum object Mathematica has such
a low performance.
Note that in all of my recent posts dealing with integrals of rational
functions in which Integrate generates a RootSum and Limit fails to
get the limit behavior of the antiderivative at Infinity (ussually
zero) the other
CAS simply do not face problems.