Re: Strange answer from Eigensystem[]!
- To: mathgroup at christensen.cybernetics.net
- Subject: [mg928] Re: [mg839] Strange answer from Eigensystem[]!
- From: Xavier Joubert <xavier at seas.ucla.edu>
- Date: Wed, 3 May 1995 00:07:13 -0400
On Tue, 25 Apr 1995, GLandry wrote:
> I have a problem w/ the eigenvectors not being correct for
> the simple nonsingular matrix A (shown below). All
> eigenvalues seem correct. Two eigenvectors check, 2 don't..
>
> WHY!!??
>
> Also, using Mathcad, I get the same eigenvalues, but two
> repeated eigenvectors (rounding off 10^-4 parts):
> {.463I,0.117-0.773I,0.354+0.54I,0.212}
> {-.463I,0.117+0.773I,0.354-0.54I,0.212}
>
> The bizarre thing to me is these eigenvectors seem to check
> too!?
>
> Help! I must be missing some basic degeneracy. Unfortunately,
> so in Mma.
>
> ****************************************
> Gary Landry
> glandry at aol.com
> eelandry at utacnvx.uta.edu
> ________________
> The University of Texas at Arlington
> Department of Electrical Engineering
> Box 19016
> Arlington, TX 76019
> ****************************************
>
> A={{0,0,0,1},{0,0,-1,0},{0,0.21,0,0},{-.21,0,0,0}};
>
> MatrixForm[A]
>
> 0 0 0 1
> 0 0 -1 0
> 0 0.21 0 0
> -0.21 0 0 0
>
> {val,vec}=Eigensystem[A];
>
> (Table[{A.vec[[i]],val[[i]] vec[[i]],
> A.vec[[i]]==val[[i]] vec[[i]]},{i,4}])
>
> -17
> {{{4.17175 10 + 0. I, 0. + 0.416598 I, 0.190909 + 0. I,
>
> 0. + 0. I}, {0. + 0. I, 0. + 0.416598 I, 0.190909 + 0. I,
>
> -17
> 0. + 1.91173 10 I}, False},
>
> -17
> {{4.17175 10 + 0. I, 0. - 0.416598 I, 0.190909 + 0. I,
>
> 0. + 0. I}, {0. + 0. I, 0. - 0.416598 I, 0.190909 + 0. I,
>
> -17
> 0. - 1.91173 10 I}, False},
>
> -312
> {{0. + 1.95628 10 I, -1. + 0. I, 0. + 0. I, 0. + 0. I},
>
> {0. + 0. I, 0. + 0. I, 0. - 0.458258 I, 0. + 0. I}, False},
>
> -312
> {{0. - 1.95628 10 I, -1. + 0. I, 0. + 0. I, 0. + 0. I},
>
> {0. + 0. I, 0. + 0. I, 0. + 0.458258 I, 0. + 0. I}, False}}
>
> Chop off "little" values, then the FIRST two work as solutions....
>
> Cvec=Chop[vec];
> (Table[{A.Cvec[[i]],val[[i]] Cvec[[i]],
> A.Cvec[[i]]==val[[i]] Cvec[[i]]},{i,4}])
>
> {{{0. + 0. I, 0. + 0.416598 I, 0.190909 + 0. I, 0. + 0. I},
>
> {0, 0. + 0.416598 I, 0.190909 + 0. I, 0}, True},
>
> {{0. + 0. I, 0. - 0.416598 I, 0.190909 + 0. I, 0. + 0. I},
>
> {0, 0. - 0.416598 I, 0.190909 + 0. I, 0}, True},
>
> {{0., -1., 0., 0.}, {0, 0, 0. - 0.458258 I, 0}, False},
>
> {{0., -1., 0., 0.}, {0, 0, 0. + 0.458258 I, 0}, False}}
>
>
Here is an answer. Expecting it will help.
-------Cut here for MMa-----------
aStrange={{0,0,0,1},{0,0,-1,0},{0,0.21,0,0},{-0.21,0,0,0}};
MatrixForm[aStrange]
(*
The "strange" results arises when MMa attempts to solve
2 4
0.0441 + 0.42 lambda + lambda == 0
*)
Solve[Det[aStrange - lambda IdentityMatrix[4]] == 0]
(*
We will only get an approximation. (note : for real & complex part)
-12
{{lambda -> -1.36825 10 - 0.458258 I},
-12
{lambda -> -1.36825 10 + 0.458258 I},
-12
{lambda -> 1.36826 10 - 0.458258 I},
-12
{lambda -> 1.36826 10 + 0.458258 I}}
So that when using
*)
{aStrangeEig,aStrangeVec}=Eigensystem[aStrange]
(* one may only reach :
{{0. + 0.458258 I, 0. - 0.458258 I, 0. + 0.458258 I, 0. - 0.458258 I},
-17
{{2.18218 + 0. I, 0. + 4.5916 10 I, 0. + 0. I, 0. + 1. I},
-17
{2.18218 + 0. I, 0. - 4.5916 10 I, 0. + 0. I, 0. - 1. I},
-17
{0. + 0. I, 0. + 1. I, 0.458258 + 0. I, 0. + 4.5916 10 I},
-17
{0. + 0. I, 0. - 1. I, 0.458258 + 0. I, 0. - 4.5916 10 I}}}
(and note again that MMa truncates the values for display)
But this result is not totally bad.
The real question : How much is it bad ?
Well... to answer it we need to accept that we should not expect more than
we give. (sometimes we can but we need to know before)
We accept to approximate. Ok so let's approximate but let's do it always.
And in fact it's enough (and better) to do it only during the last step :
*)
Table[ Chop [ {aStrange.aStrangeVec[[i]] - aStrangeEig[[i]] aStrangeVec[[i]]} ]
== {{0,0,0,0}}, {i,4}]
(* returns
{True, True, True, True}
*)
(*
In fact MMa can easily reach the good answer.
Just supply the good information !
*)
a={{0,0,0,1},{0,0,-1,0},{0,21/100,0,0},{-21/100,0,0,0}};
MatrixForm[a]
Solve[Det[a - lambda IdentityMatrix[4]] == 0]
(*
I -I
{{lambda -> -- Sqrt[21]}, {lambda -> -- Sqrt[21]},
10 10
I -I
{lambda -> -- Sqrt[21]}, {lambda -> -- Sqrt[21]}}
10 10
Aaaaaaaaaah that's better !
*)
{aEig,aVec}=Eigensystem[a]
(* I bet that the following will suit more
-I -I I I
{{-- Sqrt[21], -- Sqrt[21], -- Sqrt[21], -- Sqrt[21]},
10 10 10 10
{{10 I Sqrt[21], 0, 0, 21},
{0, -10 I Sqrt[21], 21, 0},
{-10 I Sqrt[21], 0, 0, 21},
{0, 10 I Sqrt[21], 21, 0}}}
If you need to find eigenvalues, eigenvectors in many
more tricky cases, I recommend you to create a package.
You will avoid many approximations as long as you keep
a symbolic expression.
But the more complicated, the longer...
For example if you deal a lot with this kind matrices,
you'd rather create a package using some ideas from sthg like :
*)
specialDiag={{0,0,0,u},{0,0,-u,0},{0,v,0,0},{-v,0,0,0}};
secondDiag={{0,0,0,e},{0,0,f,0},{0,g,0,0},{h,0,0,0}};
Solve[Det[specialDiag - lambda IdentityMatrix[4]] == 0, lambda]
(*
{{lambda -> I Sqrt[u] Sqrt[v]}, {lambda -> -I Sqrt[u] Sqrt[v]},
{lambda -> I Sqrt[u] Sqrt[v]}, {lambda -> -I Sqrt[u] Sqrt[v]}}
*)
Solve[Det[secondDiag - lambda IdentityMatrix[4]] == 0, lambda]
(*
{{lambda -> Sqrt[e] Sqrt[h]}, {lambda -> -(Sqrt[e] Sqrt[h])},
{lambda -> Sqrt[f] Sqrt[g]}, {lambda -> -(Sqrt[f] Sqrt[g])}}
*)
{specialDiagEig,specialDiagVec}=Eigensystem[specialDiag]
(*
{{-I Sqrt[u] Sqrt[v], -I Sqrt[u] Sqrt[v], I Sqrt[u] Sqrt[v],
I Sqrt[u] -I Sqrt[u]
I Sqrt[u] Sqrt[v]}, {{---------, 0, 0, 1}, {0, ----------, 1, 0},
Sqrt[v] Sqrt[v]
*)
{secondDiagEig,secondDiagVec}=Eigensystem[secondDiag]
(*
-I Sqrt[u] I Sqrt[u]
{----------, 0, 0, 1}, {0, ---------, 1, 0}}}
Sqrt[v] Sqrt[v]
{{-(Sqrt[f] Sqrt[g]), Sqrt[f] Sqrt[g], -(Sqrt[e] Sqrt[h]),
Sqrt[f] Sqrt[f]
Sqrt[e] Sqrt[h]}, {{0, -(-------), 1, 0}, {0, -------, 1, 0},
Sqrt[g] Sqrt[g]
Sqrt[e] Sqrt[e]
{-(-------), 0, 0, 1}, {-------, 0, 0, 1}}}
Sqrt[h] Sqrt[h]
And then apply to your real numbers.
This idea can be extended as far as you want ; complete matrices,
n x n size... and not only matrices. (what is the best rule in order
to solve a x^2 + b x + c = 0 and by extension a x^4 + b x + c = 0...
or why not any polynomial... just say what are the rules...
as long as the computer will proceed the computation in an
acceptable delay...
*)
Xavier.