4
$\begingroup$

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? The result of my code

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

$\endgroup$
8
$\begingroup$

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 orthoginal 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

I 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;

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ Apart from getting the eigensystem of the pencil constructed from the two fundamental forms, the other possibility for getting the eigensystem (equivalently, the principal curvatures and principal directions) is to get the eigensystem of the Weingarten map; there, you only need to solve a symmetric eigensystem. $\endgroup$ – J. M.'s technical difficulties Mar 26 at 22:39
  • $\begingroup$ Thanks for the immediate reply! To be honest, it is the first time that I ever heard "generalized eigensystem" in this sense, so it took me some minutes to figure out what it really is. I think that both your and my method for solving the eigenvalues are theorectically correct. I just compute the matrix of Weingarten map and find the eigenvalues. And I have used my code to run your example (torus), everything is right, the eigenvectors obtained are indeed orthogonal with respect to the metric matrix in your example. But the outputed eigenvectors of my surface are not orthogonal w the 1st f.f. $\endgroup$ – R.Z Chen Mar 26 at 23:07
  • $\begingroup$ @R.ZChen the eigenvectors returned are not necessarily normalized (for a very good reason), even if the input matrix is symmetric. $\endgroup$ – J. M.'s technical difficulties Mar 26 at 23:11
  • $\begingroup$ Well, eigenvectors are not normalized is ok. But in this case eigenvectors should have an "inner product" of zero, i.e. say $p_1$ and $p_2$ are eigenvectors belonging to different eigenvalues and $M$ is the first fundamental form, I think that we should have $p_2^T\cdot M \cdot p_1=0$, which is not the case here. I am wondering if I have misunderstood something? $\endgroup$ – R.Z Chen Mar 26 at 23:24
  • 1
    $\begingroup$ Ok! Problem is solved! Thank you! I find that in my textbook the Weingarten equation is defined as a row vector multiplied by the Weingarten matrix on the right without any warning, instead of the convention that a column vector multiplied by a matrix on the left. So I have to transpose the Weingarten matrix before solving the eigenvectors, which is the same as $FF1^{-1}\cdot FF2$. $\endgroup$ – R.Z Chen Mar 27 at 11:50
7
$\begingroup$

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 *)
|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thank you for your immediate reply, now I see that the outputed eigenvector is correct. However it still remains unknown for me why the eigenvectors are not orthogonal. $\endgroup$ – R.Z Chen Mar 26 at 22:54
  • $\begingroup$ Your WW matrix is not symmetric, which leads to your eigenvectors not being orthogonal. Throw in a Print[SymmetricMatrixQ@WW] to your function to see. $\endgroup$ – MikeY Mar 27 at 15:20
4
$\begingroup$

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[].

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ (N.B. Parts of this were based on this answer, where I had used the same technique for the hyperbolic quadratic eigenproblem. That answer also mentions an alternative based on Cholesky decomposition, which is slightly faster, but less reliable.) $\endgroup$ – J. M.'s technical difficulties Mar 27 at 15:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.