Quit[];

Se trata del sistema de 4 masas iguales unidas por resortes en un círculo. El lagrangiano (adimensionalizado) es:

L = (1/2)Σ ((OverDot[(Subscript[x, i])])^2)-(1/2)[((((Subscript[x, 1])-(Subscript[x, 2])))^2)+((((Subscript[x, 2])-(Subscript[x, 3])))^2)+((((Subscript[x, 3])-(Subscript[x, 4])))^2)+((((Subscript[x, 4])-(Subscript[x, 1])))^2)]

La matriz energía cinética:

n = 4;
T = IdentityMatrix[n];
T // MatrixForm

La matriz energía potencial:

Construida a mano:

V = {{2, -1, 0, -1}, {-1, 2, -1, 0}, {0, -1, 2, -1}, {-1, 0, -1, 2}};
V // MatrixForm

Construida usando comandos fáciles de generalizar a problemas con más masas:

n = 4;
V = 2 IdentityMatrix[n] + SparseArray[{Band[{1,2}]->-1,Band[{2,1}]->-1,{1,n}->-1, {n, 1}-> -1},{n,n}];
V // MatrixForm

El sistema de autovalores y autovectores: λ = (((ω/(Subscript[ω, o])))^2)

Primer método: buscando las soluciones de Det[V - λ T] = 0 y resolviendo para cada autovalor el sistema lineal.

n = 4;
T = IdentityMatrix[n];
V = 2 IdentityMatrix[n] + SparseArray[{Band[{1,2}]->-1,Band[{2,1}]->-1,{1,n}->-1, {n, 1}-> -1},{n,n}];

solpre = Solve[Det[V - λ T] == 0, λ]
solλ = solpre[[All, 1,2]]

Para los autovalores con multiplicidad 1, fijamos la primera componente igual a 1

(* λ = 0 *)
x = {1, b, c, d};
solx = Solve[(V - solλ[[1]] T) . x == 0, {b, c, d}]
x1 = x /.  solx[[1]]

(* λ = 4 *)
x = {1, b, c, d};
solx = Solve[(V - solλ[[4]] T) . x == 0, {b, c, d}]
x2 = x /.  solx[[1]]

Para el autovalor con multiplicidad 2 debemos fijar dos componentes para conseguir dos soluciones L.I. 

Lo más cómodo es elegir una componente igual a 1 y otra a 0.

(* λ = 2 - Primer autovector *)
x = {1, 0, c, d};
solx = Solve[(V - solλ[[2]] T) . x == 0, {c, d}]
x3 = x /.  solx[[1]]

(* λ = 2 - Segundo autovector*)
x = {0, 1, c, d};
solx = Solve[(V - solλ[[2]] T) . x == 0, {c, d}]
x4= x /.  solx[[1]]

Finalmente, estos son los autovectores:

x1
x2
x3
x4

El sistema de autovalores y autovectores, segundo método: reescribiendo la ecuación característica como 

V.(T^(-1))= λ x

y usando el comando Eigensystem.  Se lee todo directamente y es fácil de generalizar para cualquier número de masas.

T = IdentityMatrix[n];
V = 2 IdentityMatrix[n] + SparseArray[{Band[{1,2}]->-1,Band[{2,1}]->-1,{1,n}->-1, {n, 1}-> -1},{n,n}];
solESpre = Eigensystem[V . Inverse[T]]

Como el orden en que aparecen los autovalores no suele ser siempre el mismo (el autovalor 0 no siempre es el último valor que da Eigensystem), aquí abajo reordeno autovalores y autovectores de modo que el primero de todos corresponda siempre al modo de rotación uniforme, es decir, al que tiene λ = 0,

solLambda =Sort[solESpre[[1]]]
solEquis = solESpre[[2,Ordering[solESpre[[1]]]]]

La solución general: 

Notar que depende de 2n constantes. Eligiendo ci[[ j, . ]] ∝ (Subscript[δ, jl]) se aisla el modo l.

Notar que sumamos  el vector con las posiciones de equilibrio.

Notar que la solución del modo con frecuencia cero no es una trigonométrica sino una lineal.

ωs = Sqrt[solLambda]; (*las frecuencias*)
Xequilibrio = Table[(i-1) ((2π)/n), {i,1, n}]; (*las posiciones de equilibrio*)

X[t_, ci_]:=
 Xequilibrio + 
 ci[[1]]solEquis[[1]] t +
  Sum[solEquis[[i]] (ci[[i,1]] Cos[ωs[[i]]t]+ ci[[i,2]] Sin[ωs[[i]] t]), {i, 2, n}];

Gráficos:

Primero preparamos el escenario:

Show[
Graphics[
Table[{Hue[(i/n)], Disk[{Cos[(i-1)((2.π)/n)], Sin[(i-1)((2.π)/n)]},0.1]}, {i, 1, n}]],

Graphics[{AbsoluteThickness[2],Circle[]}],
Axes-> True,
Ticks-> False,
PlotRange-> 1.2 {{-1,1}, {-1,1}}
]

Ahora reemplazamos las posiciones de los discos por las funciones de t y hacemos la animación:

Fijo todas las constantes que acompañan a los Sin[ω t] iguales a cero, de modo que los controles se reducen a 4 constantes, cada una corresponde a la amplitud de un modo. 

Hay un control extra, Amp, que multiplica la amplitud de todos los modos por el mismo factor. 

La animación empieza con el modo '4' activo y los otros apagados.

Manipulate[
Show[
Graphics[
Table[{Hue[(i/n)], Disk[{Cos[punto = X[t, Amp{c1, {c2,0}, {c3,0},{c4,0}}][[i]]], Sin[punto]},0.1]}, {i, 1, n}]],

Graphics[{AbsoluteThickness[2],Circle[]}],
Axes-> True,
Ticks-> False,
PlotRange-> 1.2 {{-1,1}, {-1,1}}
],
{t, 0., ∞, (π/200.), Appearance-> "Open", AnimationRate-> 2},
{{Amp, 0.25}, 0., 4, 0.01, Appearance-> "Open"},
{{c1, 0}, 0, 1, 0.01, Appearance-> "Open"},
{{c2, 0}, 0, 1, 0.01, Appearance-> "Open"},
{{c3, 0}, 0, 1, 0.01, Appearance-> "Open"},
{{c4, 1}, 0, 1, 0.01, Appearance-> "Open"},
TrackedSymbols:> {Amp, t, c1, c2, c3, c3},
ControlPlacement->Left]

Para terminar, construyo una tabla con valores predeterminados de las condiciones iniciales, para hacer una animación, y en todo caso exportar el archivo gif.

Monitor[tpuntos = Table[
temp = Show[
Graphics[
Table[{Hue[(i/n)], Disk[{Cos[punto = X[t, 0.25{1, {.4,0}, {.9,0},{1,0}}][[i]]], Sin[punto]},0.1]}, {i, 1, n}]],

Graphics[{AbsoluteThickness[2],Circle[]}],
Axes-> True,
Ticks-> False,
PlotRange-> 1.2 {{-1,1}, {-1,1}}
],
{t, 0.,20 π, (π/50.)}];, {t, temp}]

Manipulate[tpuntos[[i]], {i, 1, Length[tpuntos], 1, Appearance-> "Open", AnimationRate-> 25}]

Aquí exporto los frames y hago el gif aparte. También se puede exportar directamente como archivo gif toda la tabla de gráficos.

dir = CreateDirectory["..\\Desktop\\masas\\"]

Monitor[For[i = 1, i ≤ Length[tpuntos], ++i,
itemp = tpuntos[[i]];
Export[StringJoin[dir, "frame_", If[i < 10, "0000", If[i < 100, "000", If[i < 1000, "00", "0"]]], ToString[i], ".png"], itemp, "PNG"];
];,i]

En 3D

antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]];

aro = Graphics3D[{RGBColor[0.16,0.5,0.83],Specularity[White,100], Tube[Table[{Cos[x], Sin[x], 0}, {x, 0, 2.π, 2(π/200.)}], 0.02]}];
color = {Orange, Red, Yellow, Green};

marco = Graphics3D[{{RGBColor[0.16,0.5,0.83], Line[0.75{{-1,0,0}, {1,0,0}}],
Line[0.75{{0,-1,0}, {0,1,0}}]},
{RGBColor[0.16,0.5,0.83], Line[{{-1.5,0,0}, {-1.25,0,0}}],
Line[{{0,-1.5,0}, {0,-1.25,0}}]},
{RGBColor[0.16,0.5,0.83], Line[{{1.5,0,0}, {1.25,0,0}}],
Line[{{0,1.5,0}, {0,1.25,0}}]}}];

tpr = 0.2;

graf = Show[
marco,
aro,
Graphics3D[
Table[{ {color[[i]],Specularity[White,100], Sphere[{Cos[punto = X[tpr, 0.25{1, {.4,0}, {.9,0},{1,0}}][[i]]], Sin[punto],0},0.1]}}, {i, 1, n}]],
Axes-> False,
Ticks-> False,
Boxed-> False,
Background-> Black,
PlotRange-> 2.5 {{-1,1}, {-1,1}, {-1,1}},
ViewVector-> {5,2, 3},
ViewVertical-> {0,0,1},
ViewAngle->30°,
ImageSize-> 800
]

$HistoryLength = 0; (* para que no guarde resultados y sature la memoria al hacer los gráficos 3D*)

dir = SetDirectory["..\\Desktop\\masas\\"]

nframes = 500;
Monitor[For[i = 80, i ≤ nframes, ++i,
tpr = ((10.π)/nframes)i;
graf = Show[
marco,
aro,
Graphics3D[
Table[{ {color[[i]],Specularity[White,100], Sphere[{Cos[punto = X[tpr, 0.25{1, {.4,0}, {.9,0},{1,0}}][[i]]], Sin[punto],0},0.1]}}, {i, 1, n}]],
Axes-> False,
Ticks-> False,
Boxed-> False,
Background-> Black,
PlotRange-> 2.5 {{-1,1}, {-1,1}, {-1,1}},
ViewVector-> {5,2, 3},
ViewVertical-> {0,0,1},
ViewAngle->30°,
ImageSize-> 800
];

itemp = antialias[graf, 2];

Export[StringJoin[dir, "\\frame_3D_", If[i < 10, "0000", If[i < 100, "000", If[i < 1000, "00", "0"]]], ToString[i], ".png"], itemp, "PNG"];
];,{i, itemp}]