You can't rely on floating point calculations to give exact results. See this link:
You need to use tolerances for these types of floating point comparisons. E.g.,
>> A = [3 -3 4;2 3 4;0 -1 1];
>> A2 = A*A;
>> eig(A).^2
ans =
-0.999999999999988 +18.973665961010287i
-0.999999999999988 -18.973665961010287i
1.000000000000000 + 0.000000000000000i
>> eig(A2)
ans =
-1.000000000000000 +18.973665961010266i
-1.000000000000000 -18.973665961010266i
1.000000000000001 + 0.000000000000000i
>> eig(A).^2 == eig(A2)
ans =
3×1 logical array
0
0
0
>> abs(eig(A).^2 - eig(A2)) < 1e-10
ans =
3×1 logical array
1
1
1
But, this isn't really a "proof" of course ... just a demonstration within numerical tolerances.
A proof for this case would be to use the Symbolic Toolbox. E.g., start with this:
As = sym(A);
And then work everything with As instead of A. Then see how the results compare ... (hint: you may need to use the simplify( ) function and you may need to reorder the eigenvalues before the comparison)
A general proof for an arbitraty matrix A would be starting with A*x = lambda*x on paper, multiplying both sides by A, and then seeing what you get as a result using hand calculations and simplifying.
