3
$\begingroup$

I'm solving the following eigensystem, and I get result which looks complex-valued. I expect the result to have 0 imaginary part, can anyone see a way to simplify it away?

B = {{17/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/
    3}, {1/3, 1/3, 1/3, 82/3}};
A = {{5/3, 0, 1/3, 0}, {0, 5/3, 0, 1/3}, {1/3, 0, 10/3, 0}, {0, 1/3, 
    0, 10/3}};
First /@ Eigensystem[{B, A}, 1]
$\endgroup$
5
$\begingroup$
B = {{17/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/3}, {1/3, 
    1/3, 1/3, 82/3}};
A = {{5/3, 0, 1/3, 0}, {0, 5/3, 0, 1/3}, {1/3, 0, 10/3, 0}, {0, 1/3, 0, 
    10/3}};

val = First /@ Eigensystem[{B, A}, 1];

All of the values are real

valR = val // RootReduce

enter image description here

Element[valR, Reals]

(* True *)

valR // N

(* {8.33542, {0.0311157, -0.179391, 0.00719285, 1.}} *)

However, if represented using radicals, they must use complex numbers (see Casus irreducibilis)

valR // ToRadicals

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ aha, RootReduce was the magic keyword! $\endgroup$ Sep 11 '20 at 3:36
5
$\begingroup$
Chop[N[First /@ Eigensystem[{B, A}, 1]]]

returns

{8.33542,{0.0311157,-0.179391,0.00719285,1.}}
$\endgroup$
2
  • $\begingroup$ I'm not 100% sure the imaginary part is exactly 0, or just close to 0, so was trying to use Mathematica to confirm my expectation $\endgroup$ Sep 11 '20 at 0:01
  • 1
    $\begingroup$ N[First /@ Eigensystem[{B, A}, 1],256] tries to give you the first 256 digits of each number and that shows if there is an imaginary part that it is preceded by about 256 zeros. $\endgroup$
    – Bill
    Sep 11 '20 at 1:55
4
$\begingroup$

This takes a long while, but it works if you just want to prove that the imaginary part is $0$.

B = {{17/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/
    3}, {1/3, 1/3, 1/3, 82/3}};
A = {{5/3, 0, 1/3, 0}, {0, 5/3, 0, 1/3}, {1/3, 0, 10/3, 0}, {0, 1/3, 
    0, 10/3}};
result = First /@ Eigensystem[{B, A}, 1]

Im[result] // FullSimplify
(*{0,{0,0,0,0}}*)

I let the computer run overnight to get the answer. Simplify did not get there.

$\endgroup$
2
  • 1
    $\begingroup$ I'm eagerly waiting for computers to get fast enough for FullSimplify to run faster than overnight $\endgroup$ Sep 11 '20 at 18:49
  • 1
    $\begingroup$ FullSimplify@Im@Numerator@Together@ComplexExpand@result will wrap up before you can fall asleep ;) $\endgroup$ Sep 11 '20 at 20:08
3
$\begingroup$

Try numerical evaluation with Mathematica's N[]

B = {{17/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 1/3}, {1/3, 1/3, 1/3, 
1/3}, {1/3, 1/3, 1/3, 82/3}};
A = {{5/3, 0, 1/3, 0}, {0, 5/3, 0, 1/3}, {1/3, 0, 10/3, 0}, {0, 1/3, 
0, 10/3}};
First /@ Eigensystem[{N@B, N@A}, 1]

The result should be

{8.33542, {0.0306117, -0.176485, 0.00707634, 0.983802}}
$\endgroup$
1
  • $\begingroup$ interesting, it seems applying N before Eigenvalues is better than applying N after $\endgroup$ Sep 11 '20 at 4:09
2
$\begingroup$

To piggyback off Bill's answer, one can just use CountRoots[] on the characteristic polynomial of the given matrix pencil, if one only wishes to show that the eigenvalues are all real:

CountRoots[CharacteristicPolynomial[{B, A}, x], x]
   4

One can then use RootIntervals[] to find brackets for the roots:

RootIntervals[CharacteristicPolynomial[{B, A}, x], Reals]
   {{{0, 0}, {0, 1}, {3, 4}, {4, 10}}, {{1}, {1}, {1}, {1}}}

Note that the root at $x=0$ was exactly isolated. The largest eigenvalue of the pencil would correspond to the last entry with the isolating interval $(4,10)$, which you can then give to Solve[]:

Solve[CharacteristicPolynomial[{B, A}, x] == 0 && 4 < x < 10, x, 
      Cubics -> False, Quartics -> False]
   {{x -> Root[-19440 + 76898 #1 - 28959 #1^2 + 2401 #1^3 &, 3]}}

Bob has already mentioned casus irreducibilis; to summarize, if you insist on a radical representation, then the use of a complex representation is (often) unavoidable, even if all the roots are real.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.