It seems Eigensystem[m] returns vectors that are not eigenvectors

Question

I am new to here so please forgive me if I do something wrong carelessly. I have faced a serious problem in eigensystem method, or more particular, eigenvalue. It seems that the following codes that called Eigensystem[] do not get the correct eigenvector:

 FundForm[r_, u_, v_] := 
  Module[{ru, rv, E1, F1, G1, ruu, ruv, rvv, n0, n, L2, M2, N2, FF1, 
    FF2, WW, K, H},
   ru = D[r, u];
   rv = D[r, v];
   E1 = Simplify[Dot[ru, ru]];
   F1 = Simplify[Dot[ru, rv]];
   G1 = Simplify[Dot[rv, rv]];
   ruu = D[ru, u];
   ruv = D[ru, v];
   rvv = D[rv, v];
   n0 = Cross[ru, rv];
   n = n0/Norm[n0];
   L2 = Simplify[Dot[ruu, n]];
   M2 = Simplify[Dot[ruv, n]];
   N2 = Simplify[Dot[rvv, n]];
   Print[E1, ";", F1, ";", G1];
   Print[L2, ";", M2, ";", N2];
   FF1 = ( {
      {E1, F1},
      {F1, G1}
     } );
   FF2 = ( {
      {L2, M2},
      {M2, N2}
     } );
   WW = FF2.Inverse[FF1];
   K = Simplify[(L2*N2 - M2^2)/(E1*G1 - F1^2)  ];
   H = Simplify[(E1 N2 - 2 F1 M2 + G1 L2)/(2 (E1*G1 - F1^2)) ];
   Print[K, ";", H];
   Print[Simplify[Eigensystem[WW]]];
   Print[WW.Transpose[{Eigenvectors[WW][[1]]}] === 
     Eigenvalues[WW][[1]] Transpose[{Eigenvectors[WW][[1]]}]];
   ];
$Assumptions = 
  Element[a, Reals] && Element[b, Reals] && Element[u, Reals] && 
   Element[v, Reals];
FundForm[{a (u + v), b (u - v), 4  u v}, u, v]

Please let me explain this code first, I am learning differential geometry now and I wanted to solve the first and second fundamental form of a 3-D surface, moreover I wanted to compute the Gauss/Mean/Principal curvatures and the principal directions. So I used the Eigensystem[m] to solve the principal directions. However I found that the outputed principal directions are not orthogonal under the metric matrix $\left(\begin{array}{cc}E& F\\F&G\end{array}\right)$ (they should be orthogonal!).

So I was wondering what is wrong with my codes. After some debugging, it turned out that Eigensystem[m] may returned the wrong eigenvectors (I tested if $W\cdot v=\lambda v$ and Mathematica returned a False)! I tried shutting down the Mathematica and open it again but thing remains the same. Now I am really frustrated, I do not know whether I am crazy or there is a bug lying in Eigensystem[m] indeed.

The version of mathematica installed on my PC is 12.0 and my PC runs Windows 10. "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)"

Any help would be appreciated! Also can anybody choose the right tag for this post?

EDIT: Thanks to @HenrikSchumacher, I find the problem lying in the code, the Weingarten matrix $WW$ should be $FF1^{-1}\cdot FF2$ instead of the reverse, also, I am not testing the eigenvector correctly, thanks to @MikeY

Answer 1

Classical problem. You want to compute the eigensystem of the second fundamental form with respect to the first fundamental form. Thus you have to solve a generalized eigensystem. This can be done with Eigensystem[{FF2, FF1}], but it does not work very well with symbolic functions. Also, one has to normalize the eigenvectors for some reason I don't get. But they are orthogonal with respect to FF1.

Here a concrete example for a torus with radii $2$ and $3$:

R1 = 2;
R2 = 3;
r = {Cos[u] (R2 + R1 Cos[v]), (R2 + R1 Cos[v]) Sin[u], R1 Sin[v]};
Dr = D[r, {{u, v}, 1}];
FF1 = Transpose[Dr].Dr // Simplify;
n = Simplify[#/Sqrt[#.#] &[Cross @@ Transpose[Dr]]];
FF2 = n.D[r, {{u, v}, 2}] // Simplify;
{κ, e} = Simplify[Eigensystem[{FF2, FF1}]];
e = #/Sqrt[#.FF1.#] & /@ e;

Now:

e.FF1.Transpose[e]
κ

{{1, 0}, {0, 1}}

{-(1/2), -(Cos[v]/(3 + 2 Cos[v]))}

Sanity check:

For v -> 0, the principal curvatures should be -1/R1 and -1/(R2+R1) (because the normal points outward). Let's see:

(κ /. v -> 0) == {-1/R1, 1/(-R1 - R2)}

True

For v -> Pi, the principal curvatures should be -1/R1 and 1/(R2-R1) (because the normal points outward):

(κ /. v -> Pi) == {-1/R1, 1/(R2-R1)}

True

Checking Gauß and mean curvature:

Det[FF2]/Det[FF1] == Times @@ κ // Simplify
Tr[FF2.Inverse[FF1]] == Total[κ] // Simplify

True

True

Also the principal curvature direction are tangent to the coordinate lines. And e does show that:

e

{{0, 1/2}, {1/Sqrt[(3 + 2 Cos[v])^2], 0}}

How to avoid generalized eigensystems

It matters from which side you multiply the inverse of the metric:

Here a random positive definite matrix `FF1` and symmetric matrix `FF2`:


d = 2;
SeedRandom[1];
FF1 = #.#\[Transpose] &@RandomReal[{-1, 1}, {d, d}];
FF2 = # + #\[Transpose] &@RandomReal[{-1, 1}, {d, d}];

Now compare

{λwrong, ewrong} = Eigensystem[FF2.Inverse[FF1]];
ewrong.FF1.Transpose[ewrong] // Chop

{{0.372482, 0.80143}, {0.80143, 1.73043}}

to

{λ, e} = Eigensystem[Inverse[FF1].FF2];
e.FF1.Transpose[e] // Chop

{{0.00166469, 0}, {0, 1.35961}}

Only the latter is diagonal. So your code should work as expect if you replace WW = FF2.Inverse[FF1]; by WW = Inverse[FF1].FF2;

Answer 2

I just made a small edit to your code. You weren't properly testing your eigensystem. See the MMA help on it.

FundForm[r_, u_, v_] := 
  Module[{ru, rv, E1, F1, G1, ruu, ruv, rvv, n0, n, L2, M2, N2, FF1, 
  FF2, WW, K, H, evals, evecs}, ru = D[r, u];
  rv = D[r, v];
   E1 = Simplify[Dot[ru, ru]];
  F1 = Simplify[Dot[ru, rv]];
  G1 = Simplify[Dot[rv, rv]];
  ruu = D[ru, u];
  ruv = D[ru, v];
  rvv = D[rv, v];
  n0 = Cross[ru, rv];
  n = n0/Norm[n0];
  L2 = Simplify[Dot[ruu, n]];
  M2 = Simplify[Dot[ruv, n]];
  N2 = Simplify[Dot[rvv, n]];
  FF1 = ({{E1, F1}, {F1, G1}});
  FF2 = ({{L2, M2}, {M2, N2}});
  WW = FF2.Inverse[FF1];
 {evals, evecs} = Eigensystem[WW];
 WW.Transpose@evecs == (Transpose@evecs).DiagonalMatrix[evals] // Simplify
 ];

and

$Assumptions = Element[a, Reals] && Element[b, Reals] && Element[u, Reals] && 
 Element[v, Reals];
FundForm[{a (u + v), b (u - v), 4 u v}, u, v]

(* True *)

Answer 3

This was supposed to be a comment to Henrik's addendum in his answer on how to get the eigensystem of a symmetric-definite pencil, but it got too long.

To keep things concrete, here is the pencil I will use in the following demo:

{m1, m2} = {HilbertMatrix[2], Array[Min, {2, 2}]};

Of course, Eigensystem[] can handle this pencil directly:

Eigensystem[{m1, m2}] // RootReduce
   {{1/6 (4 + Sqrt[13]), 1/6 (4 - Sqrt[13])},
    {{1/3 (-5 - Sqrt[13]), 1}, {1/3 (-5 + Sqrt[13]), 1}}}

and of course, one could consider the equivalent solution (using LinearSolve[m2, m1] instead of Inverse[m2].m1], as is good linear algebra practice):

Eigensystem[LinearSolve[m2, m1]] // RootReduce

which has the same result as above.

The usual concern about this procedure is that even if m1 is symmetric, and m2 is symmetric positive-definite, LinearSolve[m2, m1] is not symmetric at all!

LinearSolve[m2, m1]
   {{3/2, 2/3}, {-1/2, -1/6}}

Since methods for symmetric eigenproblems are (generally) more efficient and reliable than methods for unsymmetric eigenproblems, there is interest in constructing a symmetric eigenproblem that is equivalent to a symmetric definite pencil. I will present one such method.

First, compute the eigensystem of m2:

{vat, vet} = Eigensystem[m2];

Then, construct the following matrix:

mt = Transpose[Orthogonalize[vet]].DiagonalMatrix[Sqrt[vat]];

Construct a LinearSolveFunction[] out of this matrix:

lf = LinearSolve[mt];

and then perform the following similarity transformation on m1:

mt = lf[Transpose[lf[m1]]] // RootReduce
   {{1/15 (10 - 2 Sqrt[5]), -(7/(6 Sqrt[5]))},
    {-(7/(6 Sqrt[5])), 1/15 (10 + 2 Sqrt[5])}}

which is manifestly symmetric. We can now use Eigensystem[] on it:

{vals, vecs} = Eigensystem[mt] // RootReduce
   {{1/6 (4 + Sqrt[13]), 1/6 (4 - Sqrt[13])},
    {{1/7 (4 - Sqrt[65]), 1}, {1/7 (4 + Sqrt[65]), 1}}}

Notice that we have already obtained the eigenvalues. (If that's all you want, use Eigenvalues[] instead, of course.) To get the eigenvectors of the pencil, we need to do some more work:

vecs = #/Last[#] & /@ Transpose[lf[Transpose[vecs], "T"]] // RootReduce
   {{1/3 (-5 - Sqrt[13]), 1}, {1/3 (-5 + Sqrt[13]), 1}}

and we get the same results as the one where we directly got the pencil's eigensystem. Here, I use the undocumented method to solve a transposed linear system using LinearSolve[].