Re: Root function
- To: mathgroup at smc.vnet.net
- Subject: [mg51033] Re: [mg51009] Root function
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 2 Oct 2004 03:18:07 -0400 (EDT)
- References: <200410010849.EAA11569@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Ben wrote:
> When I try to diagonalize this matrix in mathematica 5:
> Eigensystem[{{-1, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 0}, {0, 0, 1, 0}}]
> I get this ugly mess
>
> \!\({{Root[1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 1], Root[1 + 6\ #1 +
> 5\ #1\^2 + #1\^3 &, 2], Root[
> 1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 3], 0}, {{\(-1\) - 3\
> Root[1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 1] - Root[1 + 6\ #1 +
> etc etc etc
>
> Mathematica 4 gives me the right answer.
> How do I get v5 to not spit this out and what does it mean?
The Root form is a concise way of expressing algebraic numbers via the
minimal polynomial they satisfy, along with a canonical ordering in the
complex plane (hence those numbers 1-3 in the second field).
Let's look a bit more closely at your example. Start with
mat = {{-1, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 0}, {0, 0, 1, 0}};
In[2]:= InputForm[Timing[es = Eigensystem[mat]]]
Out[2]//InputForm=
{0.05*Second, {{Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0],
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0], Root[1 + 6*#1 + 5*#1^2 +
#1^3 & ,
3, 0], 0}, {{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0], 1},
{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0], 1},
{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0], 1}, {0, 0, 0, 1}}}}
Fast, and concise. To get a radical form, you can use ToRadicals.
In[3]:= InputForm[esrad = ToRadicals[es]]
Out[3]//InputForm=
{{-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6,
-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6,
-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3, 0},
{{-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/
(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) -
(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
2*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 +
(3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) +
(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6, 1},
{-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/
(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) -
(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
2*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 +
(3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) +
(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6, 1},
{-1 - 3*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3) -
(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3)^2,
2*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3) +
(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3)^2,
-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3, 1}, {0, 0, 0, 1}}}
In[4]:= LeafCount[es]
Out[4]= 414
In[5]:= LeafCount[esrad]
Out[5]= 1284
This, by the way, is actually slightly smaller than the radical version
4 result.
In addition to size, there are other reasons to prefer the Root[] form.
(i) It is typically faster to obtain. In this example, which is fairly
easy, I see a factor of three or so gain in version 5 vs 4.2.
(ii) It is numerically more stable to evaluate. In general radical
formulations are prone to numeric problems. Root[] objects do not have
this liability.
(iii) When the roots of an irreducible cubic are all real, but not
rational, the so-called "casus irreducibilus" shows that they still must
be expressed in terms of Sqrt[-1]. This means that numeric evaluation
will give small imaginary parts unless, by happenstance, they exactly
cancel. Small numeric error from round off makes this unlikely. For example:
In[6]:= InputForm[N[es]]
Out[6]//InputForm=
{{-3.246979603717467, -1.554958132087371, -0.19806226419516174, 0.},
{{-1.801937735804838, 4.048917339522306, -3.246979603717467, 1.},
{1.2469796037174667, -0.6920214716300959, -1.554958132087371, 1.},
{-0.4450418679126288, -0.3568958678922094, -0.19806226419516174, 1.},
{0., 0., 0., 1.}}}
In[7]:= InputForm[N[esrad]]
Out[7]//InputForm=
{{-3.246979603717467 + 0.*I, -1.5549581320873713 -
2.220446049250313*^-16*I,
-0.19806226419516182 + 1.1102230246251565*^-16*I, 0.},
{{-1.801937735804838 + 0.*I, 4.048917339522306 + 0.*I,
-3.246979603717467 + 0.*I, 1.},
{1.2469796037174672 - 2.4406313453516104*^-17*I,
-0.6920214716300959 + 2.464509183785474*^-16*I,
-1.5549581320873713 - 2.220446049250313*^-16*I, 1.},
{-0.4450418679126286 - 2.8908825018377506*^-16*I,
-0.35689586789220956 + 1.7806594772125943*^-16*I,
-0.19806226419516182 + 1.1102230246251565*^-16*I, 1.}, {0., 0., 0.,
1.}}}
(iv) For sufficiently complicated algabraics, it is often faster to
evaluate the Root form numerically, at least at high precision. For example:
In[8]:= Timing[N[esrad,100];]
Out[8]= {0.02 Second, Null}
In[9]:= Timing[N[es,100];]
Out[9]= {0. Second, Null}
To answer your question, one can avoid the Root form by using ToRadicals
but for virtually all purposes I can think of you'd be best off leaving
it as is.
Daniel Lichtblau
Wolfram Research
- References:
- Root function
- From: Ben <bmauer@wisc.edu>
- Root function