Goldstein, problem #6.8

For a 3-atomic molecule, let's assume that the three atoms have coordinates {(Subscript[x, 1]), (Subscript[y, 1])}, {a + (Subscript[x, 2]), (Subscript[y, 2])} and {(Subscript[x, 3]), a + (Subscript[y, 3])}. Atom 1 will represent the right-corner atom. Then, the expression for kinetic and potential energies will become:
T = (m/2)(((\(\(x\_1\)\& . \))^2)+((\(\(y\_1\)\& . \))^2)+((\(\(x\_2\)\& . \))^2)+((\(\(y\_2\)\& . \))^2)+((\(\(x\_3\)\& . \))^2)+((\(\(y\_3\)\& . \))^2))
\(TraditionalForm\`V\  = \ \(k\/2\) \((\((L\_12 - a)\)\^2 + \((L\_13 - a)\)\^2 + \((L\_23 - a)\)\^2)\)\) 
Here (Subscript[d, ij])is a distance between i-th and j-th particles.

L12 = (Sqrt[(((a + x2 - x1))^2) + (((y2 - y1))^2)]); 
L13 = (Sqrt[(((x3 - x1))^2) + (((a + y3 - y1))^2)]);
L23 = (Sqrt[(((x3 - a - x2))^2) + (((a + y3 - y2))^2)]);

Assuming (Subscript[x, j])are small, we can expand these expressions to the first order and have:

vars = {x1, x2, x3, y1, y2, y3};
\({"δL12 = ", 
      dL12 = Plus @@ \((Simplify[\(∂\_# L12 &\) /@ 
                    vars /. \((\(\((# -> 0)\) &\) /@ vars)\), a > 0] 
              vars)\)}\)
\({"δL13 = ", \ 
      dL13 = Plus @@ \((Simplify[\(∂\_# L13 &\) /@ 
                    vars /. \((\(\((# -> 0)\) &\) /@ vars)\), a > 0] 
              vars)\)}\)
\({"δL23 = ", 
      dL23 = Plus @@ \((Simplify[\(∂\_# L23 &\) /@ 
                    vars /. \((\(\((# -> 0)\) &\) /@ vars)\), a > 0] 
              vars)\)}\)

Now we can easily find kinetic and potential energy matrices near the equilibrium:

v[x1_, x2_, x3_, y1_, y2_, 
          y3_] := (k/2) ((dL12^2) + (dL13^2) + 
              (dL23^2));
vars = {x1, y1, x2, y2, x3, y3};
vars
der = 
        Table[{vars[[i]], 
            vars[[j]]}, {i, 6}, {j, 
            6}];
Vij = 
        Map[(∂_(#[[1]], #[[2]]))v[x1, x2, x3, y1, y2, y3] &, 
          der, {2}];
Tij = 
        DiagonalMatrix[{m, m, m, m, m, m}]; 
MatrixForm /@ {Tij, Vij}

As soon as (Subscript[T, ij])is a diagonal matrix, to find normal modes of the system it's suffice to find the eigenvalues of (T^(-1)) V:

A = Inverse[Tij] . Vij;
eig = Eigensystem[A];
Print["Eigenvalues:", 
        eig[[1]]];
Print["Eigenvectors:", 
        Transpose[eig[[2]]] // 
          MatrixForm];

The first three frequences represent the zero degenerate frequency mode. We can also see that displacements, corresponding to rotations and translations, belong to zero eigenspace of the operator A:

A . 
      Transpose[{
{1, 0, 1, 0, 1, 0},    
{0, 1, 0, 1, 
            0, 1},    
{0, 0, 0, 
            1, -1, 0},    
{0, 1, 0, 0, 1, 
            1},    
{1, 
            0, 1, 1, 0, 0}   
}]

Three other modes of oscillation are:

\(TraditionalForm\`\(ω\_1 = \(\(\@\(k\/m\)\) : \ 
              e\&_\_1 = \(-x\&^\_1\) - y\&^\_1 + y\&^\_2 + x\&^\_3\);\)\)
\(TraditionalForm\`\(ω\_2 = \(\(\@\(\(2  k\)\/m\)\) : \ 
              e\&_\_2 = x\&^\_1 - y\&^\_1 - x\&^\_2 + y\&^\_3\);\)\)
\(TraditionalForm\`\(ω\_3 = \(\(\@\(\(3  k\)\/m\)\) : \ 
              e\&_\_1 = \(-x\&^\_1\) - y\&^\_1 + 2  x\&^\_2 - y\&^\_2 - 
              x\&^\_3 + 2  y\&^\_3\);\)\)

It's not difficult to visulalize the modes found. To see the animations, double-click the graphic.

visSystem[lstPoints_] := 
        Module[{ret, x1, x2, x3, y1, y2, 
            y3}, 
{x1, y1, x2, y2, x3, y3} = 
            lstPoints; 
ret = {
Line[{{x1, y1}, {x2, y2}, {x3, y3}, {x1, 
                    y1}}], 
Circle[{x1, y1}, 
                0.05], 
Circle[{x2, y2}, 
                0.05], 
Circle[{x3, y3}, 
                0.05]
}; 
Return[
            ret];
];

    
Show[Graphics[visSystem[{0, 0, 1, 0, 0, 1}]], 
      DisplayFunction -> $DisplayFunction]

pos0 = {0, 0, 1, 0, 0, 1};
ϵ = 0.1;
tMax = 24;
For[t = 0, 
        t < tMax, t++, 
Show[
GraphicsArray[
Table[Graphics[
                visSystem[
                  pos0 + ϵ eig[[2, ω + 
                          3]] Sin[((2 π ω t)/tMax)]], AspectRatio -> 1, 
                PlotRange -> {{-0.5, 1.5}, {-0.5, 1.5}}, 
                Frame -> True, FrameTicks -> False, 
                DefaultColor -> RGBColor[0.5, 0, 0.5]], {ω, 
                3}]]
]
];

The last animation represents a random oscillation of the triangle, obtained by superposition.

pos0 = {0, 0, 1, 0, 0, 1};
ϵ = 0.03;
coef = Table[ϵ Random[], {3}];
tMax = 24;
For[t = 0, 
        t < tMax, t++, 
For[i = 1; ev2 = Table[0, {6}], 
          i ≤ 3, i++, 
          ev2 += coef[[i]] 
              Sin[((2 π i t)/tMax)] 
              eig[[2, 
                i + 3]]]; 
Show[
Table[
            Graphics[visSystem[pos0 + ev2], AspectRatio -> 1, 
              PlotRange -> {{-0.5, 1.5}, {-0.5, 1.5}}, 
              Frame -> True, FrameTicks -> False, 
              DefaultColor -> RGBColor[0.5, 0, 0.5]], {ω, 
              3}]]
];