Developer Forum

Working with integration of rational functions I came to the
conclusion
that when RootSum is generated in the course of the computation
strange behavior of Mathematica (at least for me) takes place.

In a recent thread I show one such example.

As another example consider the integral of (x^3 + 6*x + 7)/(x^6 +
4*x^4 + 3*x^3 + 6)
in [0,Infinity).

The indefinite integral obtained by Mathematica is

F=Integrate[(x^3 + 6*x + 7)/(x^6 + 4*x^4 + 3*x^3 + 6), x]

RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x - #1] + 6*Log[x -
#1]*#1 + Log[x - #1]*#1^3)/(9*#1^2 + 16*#1^3 + 6*#1^5) & ]

which is of course correct.

Together[D[F, x]]

(7 + 6*x + x^3)/(6 + 3*x^3 + 4*x^4 + x^6)

Application of the Newton-Leibniz formula to get from the
andiderivative the requested
definite integral fails.

Limit[RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x - #1] + 6*Log[x
- #1]*#1 + Log[x - #1]*#1^3)/
(9*#1^2 + 16*#1^3 + 6*#1^5) & ], x -> Infinity] -
Limit[RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x - #1] +
6*Log[x - #1]*#1 + Log[x - #1]*#1^3)/
(9*#1^2 + 16*#1^3 + 6*#1^5) & ], x -> 0]

Limit[RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x - #1] + 6*Log[x
- #1]*#1 + Log[x - #1]*#1^3)/
(9*#1^2 + 16*#1^3 + 6*#1^5) & ], x -> Infinity] - RootSum[6 +
3*#1^3 + 4*#1^4 + #1^6 & ,
(7*Log[-#1] + 6*Log[-#1]*#1 + Log[-#1]*#1^3)/(9*#1^2 + 16*#1^3 +
6*#1^5) & ]

My question comes now...

Since for integrals like this Mathematica uses (as it appears in a
recent forum) the Newton-Leibniz formula
how it can evaluate the definite integral?

Integrate[(x^3 + 6*x + 7)/(x^6 + 4*x^4 + 3*x^3 + 6), {x, 0, Infinity}]

-((Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 1]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 1] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 1]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 1]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 1] + 6*Root[6 +
3*#1^3 + 4*#1^4 + #1^6 & , 1]^3))) -
(Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 2]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 2] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 2]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 2]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 2] + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 2]^3)) -
(Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 3]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 3] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 3]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 3]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 3] + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 3]^3)) -
(Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 4]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 4] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 4]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 4]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 4] + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 4]^3)) -
(Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 5]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 5] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 5]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 5]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 5] + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 5]^3)) -
(Log[-Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 6]]*(7 + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 6] +
Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 6]^3))/(Root[6 + 3*#1^3 +
4*#1^4 + #1^6 & , 6]^2*
(9 + 16*Root[6 + 3*#1^3 + 4*#1^4 + #1^6 & , 6] + 6*Root[6 + 3*#1^3
+ 4*#1^4 + #1^6 & , 6]^3))


What am I missing here?

Dimitris

ello,

I think, here is nothing very mysterious...

The antiderivative is
F = Integrate[(x^3 + 6*x + 7)/(x^6 + 4*x^4 + 3*x^3 + 6), x]

RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x - #1] + 6*Log[x - #1]*#1 + Log[x - #1]*#1^3)/(9*#1^2 + 16*#1^3 + 6*#1^5) & ]

The limit for x->Infinity is 0 as is easily obtained

RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[x] + 6*Log[x]*#1 + Log[x]*#1^3)/(9*#1^2 + 16*#1^3 + 6*#1^5) & ]

0

Therefore, the integral F is given by

-RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (7*Log[-#1] + 6*Log[-#1]*#1 + Log[-#1]*#1^3)/(9*#1^2 + 16*#1^3 + 6*#1^5) & ]

according to Newton-Leibnitz.

That this is the correct result is easily checked numerically.

Directly calculating the definite integral yields

F = Integrate[(x^3 + 6*x + 7)/(x^6 + 4*x^4 + 3*x^3 + 6), {x, 0, Infinity}]

-RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , (Log[-#1]*#1)/(9 + 16*#1 + 6*#1^3) & ] -

6*RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , Log[-#1]/(9*#1 + 16*#1^2 + 6*#1^4) & ] -

7*RootSum[6 + 3*#1^3 + 4*#1^4 + #1^6 & , Log[-#1]/(9*#1^2 + 16*#1^3 + 6*#1^5) & ]

which is equivalent to the pedestrian Newton-Leibnitz result given above.

The only surprising thing may be the fact that the integrator in this example gets the limit x-> Infinity right, while in some other
examples it does not do that. And also your way to calculate this limit using the build in command Limit[] does not seem to work
properly... But this in not a problem of RootSum, but of Limit, I think...

Regards Michael





"dimitris" < XXXX@XXXXX.COM > schrieb im Newsbeitrag news:evfkhu$73g$ XXXX@XXXXX.COM ...

lthough I understand the presence of RootSum in cases of rational
functions P(x)/Q(x) where Q(x) is a polynomial of degree 5 or higher
(recall the Abel's impisibility theorem) I believe that for degree
less than
5 (in particular for 3 and 4) Integrate could not generate RootSum.

In[307];
Quit[]

In[1]:=
$Version
Out[1]=
"5.2 for Microsoft Windows (June 20, 2005)"

Consider the function

In[2]:=
g[x_] := (2*x)/((x + 1)*(x^3 + 3*x^2 + 2*x + 1))

The following analysis show that the function is integrable
in the range [0,Infinity).

In[3]:=
Plot[g[x], {x, 0, 10}]
Out[3]=
Graphics[]

In[4]:=
(g[x] + O[x, #1]^2 & ) /@ {0, Infinity}
Out[4]=
{SeriesData[x, 0, {2}, 1, 2, 1], SeriesData[x, Infinity, {}, 2, 2, 1]}

In[5]:=
Reduce[Denominator[g[x]] == 0 && x > 0, x, Reals]
Out[5]=
False

Mathematica tries really hard but it returns the definite integral
back.

In[6]:=
Timing[Integrate[g[x], {x, 0, Infinity}]]
Out[6]=
{142.078*Second, Integrate[(2*x)/((1 + x)*(1 + 2*x + 3*x^2 + x^3)),
{x, 0, Infinity}]}

Let's obtain the (Mathematica's) antiderivative of g(x)

In[6]:=
G[x_] = Integrate[g[x], x]
Out[6]=
2*(-Log[1 + x] + RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[x - #1] +
2*Log[x - #1]*#1 + Log[x - #1]*#1^2)/
(2 + 6*#1 + 3*#1^2) & ])

Note at this point how Integrate works. First Apart is called in order
to write the integrand
as follows

In[8]:=
Apart[g[x]]
Out[8]=
-(2/(1 + x)) + (2*(1 + 2*x + x^2))/(1 + 2*x + 3*x^2 + x^3)

and then Integrate...integrates each term of the sum.

In[9]:=
(Integrate[#1, x] & ) /@ %
Out[9]=
-2*Log[1 + x] + 2*RootSum[1 - #1 + #1^3 & , (Log[1 + x - #1]*#1^2)/(-1
+ 3*#1^2) & ]

Note also that version 5.2 is more careful and doesn't return Infinity
for the definite
integral as version 4.0 does (I think because it considers the
integrals of the parts in the last output, which
both are divergent at Infinity).

So, the generation of the RootSum object in the antiderivative is due
to the failure if the Apart function to do the partial fraction
decomposition of (2*(1 + 2*x + x^2))/(1 + 2*x + 3*x^2 + x^3).

Note that I lack serious knowledge on symbolics algebra and so, I may
miss something fundamental but (although a few months ago I have a
query about this behavior of Apart) I still don't understand why there
are cases that Apart returns its argument back.

Here are two examples

In[12]:=
Apart /@ {x^2/(x^2 + 3*x + 2), x/(x^3 + 2*x^2 + x + 1)}
Out[12]=
{1 + 1/(1 + x) - 4/(2 + x), x/(1 + x + 2*x^2 + x^3)}

It seems that Apart does the PFD provided it can evaluate the roots of
the denominator of the
rational function. In the second case it can't evaluate and that's why
it fails to do the PFD.
Could/Should Apart have an option to make it do something like?

In[13]:=
Simplify[x /. Solve[1 + x + 2*x^2 + x^3 == 0, x]]
Simplify /@ Factor[1 + x + 2*x^2 + x^3, Extension -> %]
Simplify /@ Apart[x/%]

(*output is ommited*)

I think if Apart could do this, then the generation of RootSum objects
could be avoided
(especially for definite integrals where one or both limits of
integration is infinity,
and undoubtfully Integrate faces problems. Of course may be what I ask
is too difficult
and in view of the output of In[13] it seems quite complicated, and in
the very next version
RootSum at infinity does not face problems; If I knew that I would
strop right now the message! If the other CAS I use can already do it
pr