Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

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

Search the Archive

RootSum

  • To: mathgroup at smc.vnet.net
  • Subject: [mg75231] RootSum
  • From: dimitris <dimmechan at yahoo.com>
  • Date: Sat, 21 Apr 2007 23:09:45 -0400 (EDT)

Although 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
properly despite the generation of similar to RootSum objects I think
(and I hope!)
our lovely Mathematica will be able to do better!).

Returning back to the antiderivative

In[14]:=
Timing[(Limit[G[x], x -> #1] & ) /@ {0, Infinity}]
Out[14]=
{191.422*Second, {2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] +
2*Log[-#1]*#1 + Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ], 0}}

In[15]:=
(Plus[#2 - #1] & ) @@ %[[2]]
N[%]
Out[15]=
-2*RootSum[1 + 2*#1 + 3*#1^2 + #1^3 & , (Log[-#1] + 2*Log[-#1]*#1 +
Log[-#1]*#1^2)/(2 + 6*#1 + 3*#1^2) & ]
Out[114]=
0.37121697526024766 + 0.*I

Of course the result is correct but 190 seconds of timing!?!? Now you
understand
probably why I have sticked somehow with RootSum! Unecessarly very
long time
for evaluating something that at infinity is zero. (The value at zero
is returned almost
immediately).

Dimitris



  • Prev by Date: Re: how to get the table
  • Next by Date: Re: LegendreP Evaluation Mystery
  • Previous by thread: Re: RootSum
  • Next by thread: RowReduce and SparseArray