# Eigenvalues function not computing correct imaginary part

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

In:=
m4b2

Out= {{-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:= Eigenvalues[m4b2]

Out= {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:=
f[kx_] :=
0.5 (1.7537802112230758*^-15 - Sqrt[
3.0757450062177424*^-30 - 8.235683418798943*^-15 kx +
4. kx^2]);
f[1 10^-15]

Out= 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:= Eigenvalues[N[m4b2 /. kx -> 1 10^-15]]
Eigenvalues[m4b2 /. kx -> 1 10^-15]

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

Out= {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. – Samuel Flynn Feb 6 at 15:13
• You are right, I retracted the close vote. The problem is actually a precision issue. – Henrik Schumacher Feb 6 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! – Samuel Flynn Feb 6 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.