It seems Eigensystem[m] returns vectors that are not eigenvectors
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;
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 *)
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[]
.