# Eigenvalues function not computing correct imaginary part

I have the following matrix, with a free variable kx.

In[234]:=
m4b2

Out[234]= {{-1.5257*10^-16 + 1. kx, 0., 0,
1.5257*10^-16}, {-1.90635*10^-15, 1.90635*10^-15 - kx, 0., 0}, {0,
0., -1.5257*10^-16 + 1. kx, 1.5257*10^-16}, {-1.90635*10^-15, 0, 0.,
1.90635*10^-15 - kx}}


I then compute the eigenvalues as a function of kx. Two have the possibility of being imaginary. I made a function of the 3rd one. It has an imaginary part for the given value of kx

In[196]:= Eigenvalues[m4b2]

Out[196]= {1.90635*10^-15 - kx, -1.5257*10^-16 + 1. kx,
0.5 (1.75378*10^-15 - Sqrt[
3.07575*10^-30 - 8.23568*10^-15 kx + 4. kx^2]),
0.5 (1.75378*10^-15 + Sqrt[
3.07575*10^-30 - 8.23568*10^-15 kx + 4. kx^2])}

In[215]:=
f[kx_] :=
0.5 (1.7537802112230758*^-15 - Sqrt[
3.0757450062177424*^-30 - 8.235683418798943*^-15 kx +
4. kx^2]);
f[1 10^-15]

Out[216]= 8.7689*10^-16 - 5.38502*10^-16 I


But if I substitute this value of kx into the matrix and then Solve for the eigenvalues, I get the correct real part, but no imaginary part.

In[236]:= Eigenvalues[N[m4b2 /. kx -> 1 10^-15]]
Eigenvalues[m4b2 /. kx -> 1 10^-15]

Out[236]= {9.06351*10^-16, 8.7689*10^-16, 8.7689*10^-16, 8.4743*10^-16}

Out[237]= {9.06351*10^-16, 8.7689*10^-16, 8.7689*10^-16,
8.4743*10^-16}


What's going on here? The Eigenvalue method "direct" is the only one that's applicable, so it can't be that. Why does it think the imaginary part is negligible when it's the same order of magnitude as the real part?

• I don't think so? From what I understand that problem has to do with a natural ambiguity in how eigenvectors are defined and normalized. I'm willing to accept that the algebraic solution for the eigenvalues here isn't correct, but I don't see why based on that threads discussion. Commented Feb 6, 2020 at 15:13
• You are right, I retracted the close vote. The problem is actually a precision issue. Commented Feb 6, 2020 at 15:15

You work very close to $MachineEpsilon and that leads to problems. You can circumvent the issue by first scaling to "natural" sizes: A = {{-1.5257*10^-16 + 1. kx, 0., 0, 1.5257*10^-16}, {-1.90635*10^-15, 1.90635*10^-15 - kx, 0., 0}, {0, 0., -1.5257*10^-16 + 1. kx, 1.5257*10^-16}, {-1.90635*10^-15, 0, 0., 1.90635*10^-15 - kx}}; kx0 = 1. 10^-15; a = Sort@Eigenvalues[A /. kx -> kx0]; b = Sort@Eigenvalues[A] /. kx -> kx0; Min /@ DistanceMatrix[a, b]/Max[Abs[a]]  {0., 0.032504, 0.032504, 0.} Catastrophic loss of precision. Here the rescaled version: eps =$MachineEpsilon;
As = Expand[(A/eps) /. kx -> eps kxs];
kx0s = kx0/eps;
as = eps Eigenvalues[As /. kxs -> kx0s];
bs = eps Eigenvalues[As] /. kxs -> kx0s;
Min /@ DistanceMatrix[as, bs]/Max[Abs[as]]


{9.5825*10^-17, 9.5825*10^-17, 0., 0.}

• Yep, this does indeed seem to be the issue. I had an inkling precision may have been the problem, and this scaling is an elegant fix. Thank you! Commented Feb 6, 2020 at 15:26

Alternatively, with precision issues, use exact values until the last evaluation.

m4b2 = {{-1.5257*10^-16 + 1. kx, 0., 0, 1.5257*10^-16}, {-1.90635*10^-15,
1.90635*10^-15 - kx, 0., 0}, {0, 0., -1.5257*10^-16 + 1. kx,
1.5257*10^-16}, {-1.90635*10^-15, 0, 0., 1.90635*10^-15 - kx}} //
Rationalize[#, 0] &;

ev1 = Eigenvalues[m4b2] /. kx -> 10^-15 // N // Sort

(* {8.4743*10^-16, 8.7689*10^-16 - 5.38502*10^-16 I,
8.7689*10^-16 + 5.38502*10^-16 I, 9.0635*10^-16} *)

ev2 = Eigenvalues[m4b2 /. kx -> 10^-15] // N // Sort

(* {8.4743*10^-16, 8.7689*10^-16 - 5.38502*10^-16 I,
8.7689*10^-16 + 5.38502*10^-16 I, 9.0635*10^-16} *)


The results are identical

ev1 === ev2

(* True *)
`

However, exact calculations are generally much slower.