Quit[];

Mecánica Clásica 2019 2c

Algunos sistemas de la Guía 2 resueltos numéricamente. (Versión híper-comentada).

Si ya usaron el Mathematica, si son buenos entendedores o gente de poca paciencia, la explicación detallada que se da en este notebook les va a resultar cansadora: vean los notebooks mínimos subidos aparte para cada problema.

■ Recomendación general: la ayuda de los comandos del Mathematica es completísima. Incluye muchos ejemplos, y siempre se descubre una manera más eficaz o  prolija de hacer las cosas. Además es navegable. Nunca saben dónde van a terminar. Pueden apretar F1 para el menú principal de la ayuda o F1 luego de seleccionar una palabra o comando. Por ejemplo, ubiquen el cursor sobre la palabra NDSolve y tecleen F1.

■ Aquellos que no hayan usado nunca el Mathematica deben saber sólo esto: para evaluar una celda se usa el Enter del teclado numérico, o bien la combinación de teclas Shift+Enter.  Para que eso tenga efecto, el cursor tiene que estar dentro de la celda que quieren evaluar, o la celda tiene que estar seleccionada con el mouse usando las barritas que aparece en el extremo derecho de la pantalla. Pero no se apresuren aún.

■ No todas las celdas pueden evaluarse. Ésta, por ejemplo, es una celda de texto. 

■ Cada celda está abarcada por un corchete en el costado derecho de la pantalla. Así saben qué líneas incluye. Las celdas pueden agruparse jerárquicamente.

■ Un "punto y coma" luego de una serie de comandos hace que el resultado no se muestre en pantalla. Eso es útil para pasos intermedios que no estamos interesados en ver, por ejemplo asignarle cierto valor a una variable. No es útil si se trata de mostrar en pantalla el gráfico de una función.

■ En una primera lectura, evalúen las celdas con los parámetros y condiciones iniciales tal como figuran. Está chequeado que así las soluciones numéricas son bien comportadas y los resultados no producen mensajes de error. Luego experimenten.

■  Para los que no usaron antes el Mathematica, antes de pasar al par de problemas de la guía 1,  abajo tienen un ejemplo de celdas evaluables con unos cuantos comandos típicos. Prueben ejecutarlas una por una. Fíjense en la figura dónde tienen que hacer doble click para desplegar estas celdas.

Noten el formato que se usa para incluir comentarios:  (* blah blah blah *). Las celdas están agrupadas y ocultas, para no ocupar la pantalla.

Desplegar este grupo de celdas haciendo doble click en la flechita de la derecha de la pantalla.

Vayan parándose en cada celda y apretando el Enter del teclado numérico o Shift + Enter. También pueden seleccionar con el mouse el corchete que engloba a todo el grupo de celdas y ejecutarlas de una sola vez a todas.

(* Un clásico siempre vigente *)
1 + 1

(* Esta instrucción va a asignarle el valor 1 a la variable "numero" pero no va generar una salida en pantalla, debido al ";" final *)
numero = 1;

(* Aquí sí muestra el resultado. *)
(Sqrt[2])+numero

(* π entendido como número de precisión infinita. Para escribir el símbolo "π" usen [Tecla Esc]p[Tecla Esc]. Ídem para μ, ϵ, θ, y así. *)
π

(* Las primeras 50 cifras de π. *)
N[π, 50]

(* También pueden escribir la palabra Pi directamente. *)
Pi
N[Pi, 50]

(* Al multiplicar π por un número real de precisión finita, indicado por el "." luego del "1", el resultado también es de precisión finita *)
1  π
1.  π

(* Eso es válido en general. Comparen. *)
(100371/1000)

(100371./1000)

(* Esta celda contiene varias operaciones y muestra un resultado tras otro *) 

(* La constante de Euler es [Esc]ee[Esc] o E mayúscula. *) 
(* La unidad imaginaria es [Esc]ii[Esc] o I mayúscula. *)

10!
(2^(2^(2^2)))
(E^(-I π))+1 
Log[E]
(I^2)
ComplexExpand[(Sqrt[I])]

(* Una primitiva *)
Integrate[(1/u), u]

(* Una integral definida. El símbolo ∞ es [Esc]inf[Esc] *)
Integrate[((Sin[x])/x), {x, -Infinity, Infinity}]

(* Si tienen fijación con los símbolos, la integral definida se escribe con [Esc]dintt[Esc]. Raramente se usa así. Es más práctico el comando Integrate. *)
(∫_0^y)Exp[-x]\[DifferentialD]x
Integrate[Exp[-x], {x, 0, y}]

(* Los cálculos se simplifican si aclaran las cosas que el Mathematica puede asumir; si no, procede con morosa cautela. *)
(* Comparen *)

Integrate[((Exp[I x λ])/(1 + (x^2))), {x, -∞, ∞}]
Integrate[((Exp[I x λ])/(1 + (x^2))), {x, -∞, ∞}, Assumptions -> λ > 0]

Siguen varios ejemplos de cómo simplificar expresiones [ésta es una celda de texto].

(* Una expresión en bruto *)
((Cos[x])^2)+((Sin[x])^2)

(* La misma expresión pero aplicando Simplify. Ver también FullSimplify. *)
Simplify[((Cos[x])^2)+((Sin[x])^2)]

(* Para simplificar funciones trigonométricas *)
TrigReduce[Sin[x] Cos[x]]
TrigExpand[Cos[2x]]

(* El Mathematica es por demás melindroso. No va a cancelar la raíz cuadrada de una expresión al cuadrado así porque sí. *)
(* Pero pueden indicarle qué cosas asumir. Comparen los tres resultados siguientes *)
Simplify[(Sqrt[(x^2)])]  
Simplify[(Sqrt[(x^2)]), x ∈ Reals] (* El símbolo "pertenece" es [Esc]elem[Esc] *)
Simplify[(Sqrt[(x^2)]), x > 0]

(* Un gráfico simple *)
Plot[Sin[(x^2)], {x, 0, 2π}]

RandomReal[{-1, 1}] (*Un número real pseudo random entre -1 y 1. A ver si adivinan cuál será. *)

(*Una lista de 1000 números random graficados con el comando ListPlot *)

tabla = Table[RandomReal[{-1, 1}],{i, 1, 1000}]; (*Alternativamente tabla = RandomReal[{-1, 1}, 1000]; *)
ListPlot[tabla]

(* La definición de una función cuya variable no tiene atributos determinados. Es lo más usual. Notar el guion bajo y el símbolo := *)

funcion[x_]:= (x^2); 
funcion[y]
funcion[4]

(* Una función que sólo se evalúa cuando x tiene el atributo de ser un entero *)

uncion[x_Integer]:= x+(Sqrt[2]) ;
uncion[π]
uncion[1]

■ Algunos shortcuts útiles: raíz, potencia y cociente.

Ctrl + 2  -> (Sqrt[□]) 

Ctrl + 6  -> (□^□)

Ctrl + 7 -> (□/□)

La flechita es  ->[espacio].

Primer ejemplo: problema 3.



Desplegar este grupo de celdas desde la flecha de la derecha, como muestra la figura
siguiente:

Se trata del sistema de dos masas que se mueven cada una sobre un eje, unidas por una barra de longitud ℓ. Elegimos como coordenada el ángulo ϕ que forma el vector (Subscript[(Overscript[r,  ]), 21])= (Subscript[r, 2])-(Subscript[r, 1]) con el eje x.

El lagrangiano es (¡revisen!):

ℒ(ϕ, ((dϕ)/dt)) = (1/2)(((ℓ^2)(((dϕ)/dt)))^2) ((Subscript[m, 1 ])(cos^2) ϕ + (Subscript[m, 2])(sin^2)ϕ) + (Subscript[m, 2 ])g ℓ cos ϕ,

y la ecuación de movimiento:

[(Subscript[m, 1])+((Subscript[m, 2])-(Subscript[m, 1])) (sin^2)ϕ](((d^2)ϕ)/(dt^2))+((Subscript[m, 2])-(Subscript[m, 1])) sin (2ϕ)  (((((dϕ)/dt)))^2)+ (g/ℓ)   (Subscript[m, 2])sin ϕ = 0.

Introduciendo la relación μ = (Subscript[m, 1])/(Subscript[m, 2]) y el tiempo adimensional τ = (Sqrt[g/ℓ])t, queda una ecuación que depende de un único parámetro:

[μ+(1-μ) (sin^2)ϕ]ϕ'' +(1-μ) sin (2ϕ)  ((ϕ')^2)+ sin ϕ = 0.

Aquí las derivadas de ϕ son respecto a la variable τ. Si las masas son iguales, μ = 1,  resulta la ecuación del péndulo. 

Al haber introducido el tiempo τ, la escala de tiempo característica del problema es 2π, salvo casos extremos.

Para explorar las soluciones posibles ya no es necesario variar por separado (Subscript[m, 1]), (Subscript[m, 2]), g y ℓ. Todo lo que importa es μ, la relación entre las masas.

La solución numérica de la ecuación de movimiento: 

El comando esencial es NDSolve. Si buscan en la ayuda, van a encontrar numerosos ejemplos de aplicación a problemas físicos, incluidos problemas de mecánica.

Lo más concreto es proponer condiciones iniciales y resolver directamente. Va un ejemplo de prueba (pr).

La celda de abajo es la primera celda que tienen que ejecutar sí  o sí para poder seguir adelante. 

Recuerden: con el cursor dentro de la celda, tecleen Enter del teclado numérico o Shift + Enter para ejecutar.

(* Ángulo inicial de prueba *)
ϕ0pr = ((5π)/6); 

(* Velocidad angular inicial de prueba: lo leo así "phi punto 0 de prueba" *)
ϕp0pr = 0; 

(* Un valor de μ de prueba *)
μpr = 2; 

 (* Un tiempo máximo de prueba hasta el cual integrar las ecuaciones; elegido en realidad después de algunos tanteos preliminares. *)
tmaxpr = 20;

(* La cosa en sí. *)
(* Notar el uso del doble "=" para significar "igualdad" entre expresiones. Un solo "=" asigna. *) 
(* Por ejemplo, "A = 1" le asigna a "A" el valor "1". *)
(* Un doble "=" significa "igualdad" sin asignación. "A == 1" puede ser cierta o falsa, pero no le asigna a "A" el valor "1". *)

sol = 
NDSolve[
{
(μpr + (1-μpr) ((Sin[ϕ[t]])^2)) ϕ''[t] + (1-μpr) Sin[2ϕ[t]] ((ϕ'[t])^2)+ Sin[ϕ[t]] == 0, (* Ec. de movimiento, copiada de más arriba [revisen]*)
ϕ[0] == ϕ0pr, ϕ'[0] == ϕp0pr (* Condiciones iniciales *)
},
ϕ, {t, 0, tmaxpr}  (* La función para la que va a resolver y en qué variable e intervalo. Aquí es la función ϕ, de la variable t en el intervalo entre 0 y tmaxpr *)
]

Si todo salió como es debido, debería aparecer justo arriba el siguiente cuadrito

{{ϕ->(CopyTag[InterpolatingFunction[]])}}

Es la solución de NDSolve, y la hemos asignado a la variable "sol".

La solución está dada en términos de una función numérica, de la que no se dan otros detalles que su intervalo de definición y el tipo de salida que arroja, en este caso, un escalar.

Para extraer de aquí la función ϕ(t) hay que seleccionar la parte adecuada de "sol". Paso por paso nos vamos abriendo camino en la lista {{...}} hasta llegar a la función propiamente dicha:

sol[[1]]
sol[[1]][[1]]
sol[[1]][[1]][[2]] (* Con este nos quedamos. *)

La operación de seleccionar partes de una lista se hace mediante los dobles corchetes. Lo único que nos interesa es el último resultado. Es decir, en general nos bastaría con escribir:

Phi = sol[[1]][[1]][[2]];

El gráfico de la solución ϕ(t):

Phi = sol[[1]][[1]][[2]]; (* Repito la definición por las dudas de que se olviden de ejecutar la celda anterior *)

Plot[
Phi[t], {t, 0, tmaxpr},

PlotStyle->AbsoluteThickness[3],
Background->Black,
AxesStyle->White,
TicksStyle->Directive[White, 16],
ImageSize->600, 
AxesLabel->{"t", "ϕ"},
LabelStyle->20,
ImageMargins->40
]

Lo único importante en la sintaxis del Plot son las comas. Los cortes de línea no tienen ningún efecto, salvo facilitar la lectura. Todas las opciones después de "{t, 0, tmaxpr}" pueden omitirse y permutarse.

Plot[Phi[t], {t, 0, tmaxpr}]

La animación. 

Primero probamos graficar las cosas en un instante t fijo. Elegido un t de prueba, la celda de abajo combina varios gráficos: los ejes, las masas y la barra entre las masas.

(* Un valor de tiempo de prueba. Tiene que estar entre 0 y tmaxpr *)
(* Cambien luego este valor para obtener gráficos a distintos tiempos *)
tpr = 3.75; 

(* Para ese valor del tiempo, la configuración del sistema queda dada por el valor de ϕ que arroja la función Phi[tpr] *)
(* Conocido ϕpr = Phi[tpr] debería poder graficar la configuración del sistema *)
ϕpr = Phi[tpr]; 


(* El comando Show combina varios gráficos que comparten la misma región del plano *)
(* Empiezo por graficar un par de ejes, y luego agrego las masas y la barra *)
Show[ 

(* Ejes *)
ejes = Graphics[{ 
{AbsoluteThickness[2], LightBlue, Arrow[{{-1.5, 0}, {1.5, 0}}], Arrow[{{0, 1.5}, {0,-1.5}}]},
{White,Text[Style["𝓍", 30],{0.2,-0.95×1.5}]},
{White,Text[Style["𝓎", 30],{0.99 1.5,0.22}]}}
], 

(* Masas y barra *)
Graphics[{
{AbsoluteThickness[2], Blue, Line[{{-Sin[ϕpr], 0}, {0, -Cos[ϕpr]}}]},
{Green, Disk[ {-Sin[ϕpr], 0}, 0.05]},
{Red, Disk[{0, -Cos[ϕpr]}, 0.05]}}
],

(* Retoques opcionales. *)
PlotRange->1.6{{-1, 1}, {-1, 1}},
Axes-> False,
Background-> Black,
ImageSize-> 400
]

OK. Si resultó para un t dado, tiene que funcionar para un t cualquiera. Más que nada, lo de arriba sirvió para ajustar los tamaños de los discos, la extensión de los ejes, etc.

En lugar de andar cambiando el valor de tpr una y otra vez, lo más práctico es usar el comando Manipulate, que nos permite variar el tiempo y graficar la configuración de las masas de manera interactiva.

Manipulate[
ϕpr = Phi[t];

Show[
 ejes, (* El gráfico de los ejes se definió más arriba. Ahora lo uso directamente. *)

Graphics[{
{AbsoluteThickness[2], Blue, Line[{{-Sin[ϕpr], 0}, {0, -Cos[ϕpr]}}]},
{Green, Disk[ {-Sin[ϕpr], 0}, 0.05]},
{Red, Disk[{0, -Cos[ϕpr]}, 0.05]}}
],

PlotRange-> 1.6{{-1, 1}, {-1, 1}},
Axes-> False,
Background-> Black,
ImageSize-> 400
],

(* Abajo está contenida la información sobre la variable que se contrala con Manipulate. El tiempo t va desde 0 hasta tmaxpr y el paso por defecto es "tmaxpr/1000.". El "." luego del "1000." es para que no use números de precisión infinita. *)

{t, 0, tmaxpr, (tmaxpr/1000.), Appearance-> "Open", AnimationRate-> 2}  
]

¡Cautivador! Pueden mover el tiempo usando el mouse, apretar el botón Play, o ingresar un tiempo cualquiera en el cuadrito debajo de la barra. En tal caso usen el Enter normal para forzar el cambio.

Si le dan Play, al llegar al final del recorrido temporal la animación vuelve a t = 0. Sería raro que justo el tiempo tmaxpr  fuera igual a un múltiplo entero del período del movimiento. De manera que, en general, la animación no se cerrará sobre sí misma, sino que tendrá un salto desde la configuración en tmaxpr hasta la configuación en t = 0. 

Si quieren hacer un gif prolijo, esos saltos son indeseables. Es necesario que la animación termine justo un paso antes de completarse un período del movimiento.

Podríamos calcular el período del movimiento de manera analítica (lo veremos más adelante en la materia). En lugar de eso, vamos a incluir dentro de NDSolve un comando que vigile los cruces por cero del ángulo ϕ. Ese comando, que sirve también para muchas otras cosas, es WhenEvent. Ya que ϕ sí o sí cruza por cero (¿por qué?), esto servirá para medir el intervalo entre pasos sucesivos.

Cada vez que sea cierto que Sin[ϕ(t)/2] = 0, lo que implica que ϕ(t) es un multiplo entero y par de 2π, WhenEvent va a "sembrar", vía Sow,  el valor de t correspondiente. Cuando el cálculo esté completo, el comando Reap va a "cosechar" la lista de esos valores. La diferencia temporal entre estos sucesivos cruces por 0  debería ser constante e igual al período del movimiento.

En realidad hay que hacer una salvedad. La barra puede rotar (siempre en el mismo sentido) u oscilar como un péndulo. Si oscila, el período es igual a dos veces el intervalo entre pasos sucesivos por ϕ = 0 (¿por qué?). Si rota, el período es igual al intervalo entre dos pasos sucesivos. Para cubrirnos en ambos casos, definimos T = 2 veces el intervalo entre pasos sucesivos, lo que es redundante para el movimiento de rotación. Hay maneras de corregir esto.

He aquí entonces el NDSolve que nos devuelve la solución de la ecuación de movimiento y la lista de tiempos en los que ϕ cruza por cero.

(Notar que el tiempo tmax tiene que ser lo suficientemente largo para que haya varios pasos por la vertical. Si en la primera elección no funciona, tendrán que cambiarlo.)

ϕ0pr = ((5π)/6);
ϕp0pr = 0; 
μpr = 2;
tmaxpr = 50; 

sol2 =
Reap[ (* El cosechador de lo que siembra Sow más abajo*)
NDSolve[ (* El mismo sistema de ecuaciones y condiciones iniciales que antes *)
{
(μpr + (1-μpr) ((Sin[ϕ[t]])^2)) ϕ''[t] + (1-μpr) Sin[2ϕ[t]] ((ϕ'[t])^2)+ Sin[ϕ[t]] == 0, 
ϕ[0] == ϕ0pr, ϕ'[0] == ϕp0pr,
WhenEvent[Sin[ϕ[t]/2] == 0, Sow[t]] (* WhenEvent para detectar los pasos por la vertical en combinación con el sembrador Sow *)
},
ϕ, {t, 0, tmaxpr}  
]
]

Si todo salió bien, el resultado de la celda de arriba debería tener esta forma:

{{{ϕ->(CopyTag[InterpolatingFunction[]])}},{{3.3428455269521837`,10.028534627349119`,16.714221051476557`,23.399904800155852`,30.085585873357363`,36.77126427101624`,43.45693999306614`}}

Ven que el resultado contiene dos partes: la solución propiamente dicha, ϕ(t), y una serie de tiempos en los que ϕ cruza por cero. La celda de abajo extrae cada cosa por separado:

Phi = sol2[[1]][[1]][[1]][[2]]
Tiempos = sol2[[2]][[1]]

Verifiquemos que los sucesivos pasos por cero estén equiespaciados

lista = Differences[Tiempos]

ListPlot[lista, PlotRange-> {Min[lista]-((Max[lista]-Min[lista])/4),Max[lista]+((Max[lista]-Min[lista])/4)}, Frame-> True, ImageSize->500, FrameTicksStyle->16]

¡OK dentro del error! Miren la escala vertical.

Entonces para hacer la animación alcanza con variar t entre 0 y el período,  obtenido como la diferencia de dos valores consecutivos dentro de la lista de tiempos Tiempos (y multiplicado por 2 para abarcar los casos de oscilación y rotación).

Periodo =2 × Differences[Tiempos][[1]]; (* El símbolo × es redundante, pero a veces ayuda a la lectura. Se escribe como [Esc]*[Esc] *)

Ncuadros = 100.;
Δ = (Periodo/Ncuadros); 

ejes = Graphics[{ 
{AbsoluteThickness[2], LightBlue, Arrow[{{-1.5, 0}, {1.5, 0}}], Arrow[{{0, 1.5}, {0,-1.5}}]},
{White,Text[Style["𝓍", 30],{0.2,-0.9×1.5}]},
{White,Text[Style["𝓎", 30],{0.9 1.5,0.2}]}}
];

Manipulate[ (* Ahora el tiempo máximo está elegido igual al Período. El paso temporal por defecto entre cuadros sucesivos es Δ. *)

ϕpr = Phi[t];

Show[

ejes, 

Graphics[{
{AbsoluteThickness[2], Blue, Line[{{-Sin[ϕpr], 0}, {0, -Cos[ϕpr]}}]},
{Green, Disk[ {-Sin[ϕpr], 0}, 0.05]},
{Red, Disk[{0, -Cos[ϕpr]}, 0.05]}}
],

AspectRatio-> Automatic,
PlotRange-> 1.1× 1.5{{-1, 1}, {-1, 1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 400
],

{t, 0, Periodo-Δ,Δ, Appearance-> "Open", AnimationRate-> 5} 
]

¡Ahora la animación es un loop cerrado! ¡Vaya que hay diferencia!

Para exportar la animación como un gif usando el propio Mathematica. 

Primero creamos una tabla con todos los gráficos. En esencia, es sustituir en la celda de arriba Manipulate por Table.

Ncuadros = 200.;
T = Periodo;  (* Si el movimiento fuera de rotación, habría que dividir el Periodo por 2, para no hacer un gif el doble de largo de lo necesario *)
Δ = (T/Ncuadros);  

tgraficos = Table[
ϕpr = Phi[t];
Show[
ejes, 
Graphics[{
{AbsoluteThickness[2], Blue, Line[{{-Sin[ϕpr], 0}, {0, -Cos[ϕpr]}}]},
{Green, Disk[ {-Sin[ϕpr], 0}, 0.05]},
{Red, Disk[{0, -Cos[ϕpr]}, 0.05]}}
],
AspectRatio-> Automatic,
PlotRange-> 1.1× 1.5{{-1, 1}, {-1, 1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 400
],
{t, 0, T-Δ,Δ} 
];

En efecto, todos los cuadros están en la lista tgraficos, y se pueden animar directamente:

Manipulate[tgraficos[[i]], {i, 1, Length[tgraficos], 1, Appearance-> "Open", AnimationRate-> 50}]

Ahora exportamos el archivo gif, que aparecerá en el mismo directorio que el Notebook. Puede tomar hasta 1 minuto.

SetDirectory[NotebookDirectory[]];
Export["animacion_1.gif", tgraficos, "GIF", "DisplayDurations"-> 1/25];

Mientras se hace el gif,  aprovecho para decir que lo mejor suele ser exportar cada figura a un archivo, y luego hacer el gif en otro programa.

El tamaño en kB del gif generado:

Round[FileByteCount["animacion_1.gif"]/1000]

Segundo ejemplo:  problema 1.



Incluye comandos apenas más avanzados, opciones 3D y un exportador de gráficos. Se da por conocido lo explicado en el problema anterior.

Expandir las celdas haciendo doble click en la flechita (por demás pequeñita) de la derecha.

Lo que sigue: 

• La solución numérica. 
• La extracción del período.
• El gráfico en bruto de las funciones ϕ[t] e y[t].

Antes que nada un Quit. Reinicia el kernel del Mathematica, libera la memoria y borra las definiciones y resultados anteriores.

Quit[];

El comando Module aquí abajo encierra todo el cálculo de la solución numérica de la ecuación de movimiento. 

Module sirve para definir variables locales. Así se evita el peligro de definiciones múltiples para un mismo símbolo. Es práctica recomendada. Todo lo que está definido dentro de las primeras llaves tiene significado únicamente dentro del bloque encerrado por Module.

Module[ 
{ϕ, 
y,

(*Ccondiciones iniciales*)
ϕ0 =π/2,
ϕp0 = -1.42,
y0 = 0,

(*  Hay una condición inicial suplementaria que está implícita en las ecuaciones de movimiento dentro del NDSolve; a saber, (Subscript[P, y]) = 0. *)
(*  La elección (Subscript[P, y]) = 0 no significa ninguna pérdida de generalidad. Corresponde a tomar un SR que se mueve con la velocidad horizontal del CM. *)

(* Estoy tomando m2 = 1  y g/l = 1 *)
m1 = 5, 

tmax = 30
},

(* Solución numérica y cosecha de los tiempos de paso por la vertical, para medir el período *)
{sol, cruces} =Reap[NDSolve[
 {Sin[ϕ[t]]+(( Sin[2ϕ[t]] ((ϕ'[t])^2))/(2(1+m1)))+(1-(((Cos[ϕ[t]])^2)/(1+m1))) ϕ''[t] == 0,y'[t] ==-(1/(1+m1))Cos[ϕ[t]] ϕ'[t], ϕ[0] == ϕ0, ϕ'[0] == ϕp0, y[0] == y0,
WhenEvent[Sin[ϕ[t]/2] == 0 , Sow[t]]}, {ϕ, y}, {t, 0, tmax}]];
];

Print["Phi = ",Phi = sol[[1]][[1]][[2]], "\n"];
Print["Y = ", Y = sol[[1]][[2]][[2]], "\n"];
Print["Cruces por cero (a menos de múltiplos de 2π) en t -> ",  cruces[[1]], "\n"];
Print["Δt entre cruces consecutivos -> ", Differences[cruces[[1]]], "\n"];
Print["Periodo (o 2 veces el período si el movimiento es de rotación): ", T =2 Differences[cruces[[1]]][[1]]];
{Plot[Phi[t], {t, 0,2T}, PlotStyle->Thick, ImageMargins->30, ImageSize->400, PlotLabel-> Style["ϕ(t) \n", 18]],
Plot[Y[t], {t, 0,2T}, PlotStyle->Thick, ImageMargins->30,ImageSize->400, PlotLabel-> Style["y(t) \n", 18]]}

Abajo, la animación:

(* El bloque DynamicModule  es igual que el Module pero para celdas que se evalúan dinámicamente. Puede omitirse sin peligro si no repiten definiciones de los símbolos. *) 

DynamicModule[ 

{y, ϕ, t, ran = 2, Δ = T/100, ejes, orb}, (* Todas estas variables son locales. Fijo en T/100 el intervalo entre dos cuadros sucesivos. *)

ejes = Graphics[{AbsoluteThickness[2], White,Line[{{-ran, 0}, {ran, 0}}],Line[{{0, -ran}, {0,ran}}]}];
orb = ParametricPlot[{Y[τ]+Sin[Phi[τ]], -Cos[Phi[τ]]}, {τ, 0, T}, PlotStyle->{Dotted,Orange}]; (* La órbita de la segunda masa *)

Manipulate[

y = Y[t]; 
ϕ = Phi[t];

Show[
ejes,
orb,
Graphics[{AbsoluteThickness[2], Blue, Line[{{y, 0}, {y + Sin[ϕ], -Cos[ϕ]}}]}],
Graphics[{Green, Disk[ {y + Sin[ϕ], -Cos[ϕ]}, 0.1]}],
Graphics[{Red, Disk[{y, 0}, 0.1]}],
AspectRatio-> Automatic,
PlotRange-> ran{{-1, 1}, {-1, 1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 600
],

(*Abajo: el tiempo máximo llega hasta T - Δ, porque si llegase hasta T la animación repetiría el primer cuadro, ya que en t = 0 y t = T la configuración es la misma. *)

{t, 0, T-Δ, Δ, Appearance-> "Open"},

TrackedSymbols:> {t}   (*Sólo recalcula las cosas cuando cambia t.*)
]
]

Ya está ajustado para que el loop sea cerrado. Lo podemos subir a Tumblr sin que nos avergüence (al menos no por esa razón).

En (SuperStar[glorioso]) 3D, pero nada óptimo.

DynamicModule[
{y, ϕ, t, ran = 2, Δ = T/200, ejes, orb}, 
ejes = Graphics3D[{ White,Specularity[White,100],Tube[{{0,-ran, 0}, {0,ran, 0}}, 0.012],Tube[{{0,0, -ran}, {0,0,ran}}, 0.012]}];
orb = ParametricPlot3D[{0, Y[τ]+Sin[Phi[τ]], -Cos[Phi[τ]]}, {τ, 0, T}, PlotStyle->{ AbsoluteDashing[{5, 8}],Orange}];

Manipulate[
y = Y[t]; 
ϕ = Phi[t];

Show[
ejes,
orb,

Graphics3D[{AbsoluteThickness[1], Blue, Specularity[White,20], Tube[{{0,y, 0}, {0,y + Sin[ϕ], -Cos[ϕ]}}, 0.03]}],
Graphics3D[{Green,Specularity[White,100],Sphere[ {0,y + Sin[ϕ], -Cos[ϕ]}, 0.08]}],
Graphics3D[{Red, Specularity[White,100],Sphere[{0,y, 0}, 0.08]}],

AspectRatio-> Automatic,
PlotRange-> ran{{-0.5, 0.5}, {-1, 1}, {-1,1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 600,
ViewVector-> {8{2, -0.9, 0.5}, {0,0,0}},
ViewAngle-> 20°
],

{t, 0, T-Δ, Δ, Appearance-> "Open", AnimationRate-> 2}, 
TrackedSymbols:> {t}
]
]

(SuperStar[])Al menos hasta la versión 10.0, que es la que estoy usando, los gráficos 3D del Mathematica no se caracterizan por un rendering muy delicado ni muy veloz. 

Si todo lo que interesa es la velocidad, ésta puede mejorarse si calculan antes todos los cuadros y los guardan como imágenes 2D. Eso es lo que hace la tabla de abajo. Las imágenes se guardan en la tabla llamada tgraf. Para un gráfico con tan pocos elementos, la ganancia no es muy notable.

(* Puede tardar varios segundos*) 

Module[
{y, ϕ, t,τ,  ran = 2, Δ = T/100, ejes, orb, imagentemp}, 

ejes = Graphics3D[{ White,Specularity[White,100],Tube[{{0,-ran, 0}, {0,ran, 0}}, 0.012],Tube[{{0,0, -ran}, {0,0,ran}}, 0.012]}];
orb = ParametricPlot3D[{0, Y[τ]+Sin[Phi[τ]], -Cos[Phi[τ]]}, {τ, 0, T}, PlotStyle->{ AbsoluteDashing[{5, 8}],Orange}];

tgraf = 
Table[
y = Y[t]; 
ϕ = Phi[t];

Image[Show[
ejes,
orb,

Graphics3D[{AbsoluteThickness[1], Blue,Specularity[White,20], Tube[{{0,y, 0}, {0,y + Sin[ϕ], -Cos[ϕ]}}, 0.03]}],
Graphics3D[{Green,Specularity[White,100],Sphere[ {0,y + Sin[ϕ], -Cos[ϕ]}, 0.08]}],
Graphics3D[{Red, Specularity[White,100],Sphere[{0,y, 0}, 0.08]}],

AspectRatio-> Automatic,
PlotRange-> ran{{-0.5, 0.5}, {-1, 1}, {-1,1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 600,
ViewVector-> {8{2, -0.9, 0.5}, {0,0,0}},
ViewAngle-> 20°
]],

{t, 0, T-Δ, Δ}
];
]

Ahora puede hacerse la animación con todo ya calculado y guardado previamente:

Manipulate[
tgraf[[i]], 
{i, 1, Length[tgraf], 1, Appearance-> "Open"}
]

El movimiento es más fluido pero los gráficos siguen siendo espantosos.

Para obtener buenos gráficos 3D hay que hacer primero los gráficos en un tamaño exagerado y luego reducirlos. Así obtienen una resolución aceptable.

Eso es lo que hace la función antialias definida aquí abajo. Hace un imagen 2D n veces más grande de lo que queremos, vía el comando Rasterize, y luego la reduce.

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

Ejemplo: primero va el gráfico en bruto y luego el antialiasiado. Comparen.

antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]]; (*repito la definición por las dudas*)

Module[
{y, ϕ, t =2, ran = 2,  orb, ejes},

{y = Y[t], ϕ = Phi[t]};
orb = ParametricPlot3D[{0, Y[τ]+Sin[Phi[τ]], -Cos[Phi[τ]]}, {τ, 0, T}, PlotStyle->{ AbsoluteDashing[{5,8}],Orange}];
ejes = Graphics3D[{ White,Specularity[White,100],Tube[{{0,-ran, 0}, {0,ran, 0}}, 0.012],Tube[{{0,0, -ran}, {0,0,ran}}, 0.012]}];

bruto = Show[
orb,
ejes,
Graphics3D[{AbsoluteThickness[1], Blue,Specularity[White,20], Tube[{{0,y, 0}, {0,y + Sin[ϕ], -Cos[ϕ]}}, 0.03]}],
Graphics3D[{Green,Specularity[White,100],Sphere[ {0,y + Sin[ϕ], -Cos[ϕ]}, 0.08]}],
Graphics3D[{Red, Specularity[White,100],Sphere[{0,y, 0}, 0.08]}],
AspectRatio-> Automatic,
PlotRange-> ran{{-0.5, 0.5}, {-1, 1}, {-1,1}},
Axes-> False,
Ticks-> False,
Background-> Black,
BoxStyle-> {White, AbsoluteThickness[3]},
ImageSize-> 600,
ViewVector-> {8{2, -0.9, 0.5}, {0,0,0}},
ViewAngle-> 17°, 
PlotRangePadding-> None
]
];

antibruto = antialias[bruto, 3];

ImageAssemble[{{Image[bruto], antibruto}}]

Noten sobre todo cómo cambian las líneas rectas entre una figura y la otra. Noten también que rescalear y antialiazar el gráfico no es un proceso rápido.

No exageren con la ampliación auxiliar. Aquí un valor de n = 3 significa que en los pasos intermedios el Mathematica maneja imágenes de 1800 px, ya que el original tiene 600 px

Línea de producción para hacer un gif como la gente.

Cada gráfico es sometido al proceso de antialias.

Presten atención al uso del comando Monitor, que les permite ir viendo los resultados, y, en caso de que esté dando cualquier cosa, pueden interrumpir la evaluación ("Alt + .") sin necesidad de esperar a que llegue a su término natural. 

Puede tardar varios minutos.

ATENCIÓN: A esta altura ya sabrán si el movimiento que van a animar es de rotación o de oscilación. Para no calcular cosas de más, en el primer caso habrá que redefinir Tverdadero  = T/2.

antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]];
$HistoryLength = 0; (* Prevención válida para que no queden resultados innecesarios almacenados en el Historial *)
Clear[AnimaTabla]; (* En prevención de que se haya calculado antes y ocupe mucha memoria. *)


Module[
{Tverdadero = T/2, (*dejar "T" si el movimiento es de oscilación*)
 Ncuadros = 100, (* Disminuyan el número de cuadros si tienen poca paciencia. *)
 y, ϕ, ran = 2, Δ, pasos = 0, orb,ejes,  bruto=Graphics[], rescaleado=Graphics[]},

Δ = Tverdadero/Ncuadros; (*el paso temporal entre cuadros sucesivos*)

orb = ParametricPlot3D[{0, Y[τ]+Sin[Phi[τ]], -Cos[Phi[τ]]}, {τ, 0, Tverdadero}, PlotStyle->{ AbsoluteDashing[{5, 8}],Orange}];
ejes = Graphics3D[{ White,Specularity[White,100],Tube[{{0,-ran, 0}, {0,ran, 0}}, 0.012],Tube[{{0,0, -ran}, {0,0,ran}}, 0.012]}];

Monitor [ (*este comando permite ir viendo las variables que uno elija [que se indican al final de todo] a medida que se ejecuta la celda*)

AnimaTabla =  (*en esta tabla voy almacenando los cuadros para hacer luego la animación*)

 Table[ 
{y = Y[t], ϕ = Phi[t]};

bruto = Show[
orb,
ejes,
Graphics3D[{AbsoluteThickness[1], Blue,Specularity[White,20], Tube[{{0,y, 0}, {0,y + Sin[ϕ], -Cos[ϕ]}}, 0.02]}],
Graphics3D[{Green,Specularity[White,100],Sphere[ {0,y + Sin[ϕ], -Cos[ϕ]}, 0.08]}],
Graphics3D[{Red, Specularity[White,100],Sphere[{0,y, 0}, 0.08]}],
AspectRatio-> Automatic,
PlotRange-> ran{{-0.5, 0.5}, {-1, 1}, {-1,1}},
Axes-> False,
Ticks-> False,
Background-> Black,
BoxStyle-> {White, AbsoluteThickness[3]},
ImageSize-> 600,
ViewVector-> {8{2, -0.9, 0.5}, {0,0,0}},
ViewAngle-> 17°, 
PlotRangePadding-> None
];

++pasos; (*voy llevando la cuenta de los cuadros para que los muestre Monitor*)

rescaleado = antialias[bruto, 3],

{t, 0, Tverdadero-Δ,Δ}
]; ,

(* Abajo: las variables que muestra el comando Monitor. Incluyo la Memoria en Uso para tener una idea de si va a alcanzar la memoria RAM. Caso contrario, habrá que ir guardando las figuras en archivos.*)

{pasos, Show[bruto, ImageSize-> 400], Show[rescaleado, ImageSize-> 400], MemoryInUse[]}  
]
]
$HistoryLength = ∞;

Abajo: la animación con los gráficos remasterizados.

Manipulate[
AnimaTabla[[i]], {i, 1, Length[AnimaTabla], 1, Appearance-> "Open"}]

¡Qué diferencia, señores!

Exporto un gif con los cuadros de la animación guardados en AnimaTabla. Puede tomar de 1 a 3 minutos.

SetDirectory[NotebookDirectory[]];
AbsoluteTiming[
Export["animacion_3D.gif", AnimaTabla, "GIF",  "DisplayDurations"-> 1/25];
]

Alternativamente, para poder editar los cuadros por fuera del Mathematica, los guardo en un directorio aparte. Luego puedo manipularlos en otro programa.

(* Presten atención al uso de las dobles barras para indicar la ruta. Si están usando Linux deberán cambiarlas por / o quizá  // *)

SetDirectory["C:\\users\\Juan Zanella\\Desktop\\cuadros\\"]; 

Monitor[
For[i = 1, i ≤ Length[AnimaTabla], ++i, (* El comando For tiene una sintáxis transparente*)

Export[StringJoin["frame_", If[i < 10, "000", If[i < 100, "00", "0"]], ToString[i], ".png"], AnimaTabla[[i]], "PNG"];
];, 
i (* "i" es la variable que muestra el comando Monitor; la cosa no se acaba hasta que "i" no llega al número de cuadros guardados en AnimaTabla*)
]

Dentro del Export, la función del StringJoin es construir el nombre de cada archivo. El par de If agrega ceros antes del número de cuadro para que luego se ordenen correctamente. El comando Monitor nos permite planificar en qué usaremos los próximos segundos, minutos u horas, en el peor de los casos.

Detrás de la escena: el Lagrangiano y las ecuaciones de movimiento a través del Mathematica. (expandir la celda)

Problema 5:

Quit[];

(* Posiciones en términos de la coordenada generalizada ϕ *)
r1[ϕ_]:=l {-Sin[ϕ], 0};
r2[ϕ_]:= l{0, -Cos[ϕ]};

(* Energía cinética y energía potencial *)
T[ϕ_, ϕp_]:=Module[{Phi}, (* variable auxiliar para hacer reemplazos *)
 (1/2)m1 (D[r1[Phi[t]], t] . D[r1[Phi[t]], t]) +  (1/2)m2 (D[r2[Phi[t]], t] . D[r2[Phi[t]], t])  /. {Phi[t]-> ϕ, Phi'[t]-> ϕp}];
V[ϕ_]:= m2 g l  Cos[ϕ];

(* Lagrangiano *)
L[ϕ_, ϕp_]:= T[ϕ, ϕp] - V[ϕ];

Collect[L[ϕ, ϕp], {ϕp, m1, m2}, Simplify]

Ecuación de movimiento en forma directa. Simplificada paso a paso como hicimos en clase.

ec = Simplify[D[D[L[ϕ[t], ϕp[t]], ϕp[t]], t] - D[L[ϕ[t], ϕp[t]], ϕ[t]] /. {ϕp[t]-> ϕ'[t], ϕp'[t]-> ϕ''[t]}]

(* Para hacerla más legible abrevio ϕ[t] como ϕ, ϕ'[t] como ϕ' y ϕ''[t] como ϕ'' *)

ecsimple = ec /. {ϕ[t]-> ϕ, ϕ'[t]-> ϕ', ϕ''[t]-> ϕ''}

(* Defino la frecuencia g = (ω^(2 ))l, y divido la ecuación por (l^2).  Se va l. *)
ecsimple2 = Simplify[(1/(l^2))(ecsimple/. g-> l (ω^2))]

(* Luego defino m1 = μ m2 y divido la ecuación por m2. Se va m2 *)
ecsimple3 = Simplify[(1/m2)(ecsimple2/. {m1-> μ m2})]

(* Por último, usar τ = ω t en lugar de t equivale a definir ω = 1. Se va ω. *)
ecfinal=ecsimple3 /. ω-> 1

Ecuación de movimiento obtenida usando el comando EulerEquations:

Needs["VariationalMethods`"]

ecEE = EulerEquations[L[ϕ[t], ϕ'[t]], ϕ[t], t][[1]]

(* La simplificamos como antes *)

ecEEfinal = Simplify[(1/((l^2)m2))(ecEE /.  {ϕ[t]-> ϕ, ϕ'[t]-> ϕ', ϕ''[t]-> ϕ'', g-> l (ω^2), m1-> μ m2})/. ω-> 1]

(*Salvo el signo, es el mismo resultado de antes *)

Simplify[ecEEfinal + ecfinal]

Problema 17:

Quit[];

(* Posiciones en términos de las coordenadas generalizadas {y, ϕ} *)
r1[y_, ϕ_]:= {y, 0};
r2[y_, ϕ_]:= r1[y, ϕ] + l{-Cos[ϕ], Sin[ϕ]};

(* Energía cinética y energía potencial *)
T[y_, ϕ_, yp_, ϕp_]:=Module[{Phi, Y}, (* variables auxiliares para hacer reemplazos *)
 (1/2)m1 (D[r1[Y[t], Phi[t]], t] . D[r1[Y[t], Phi[t]], t])  +  (1/2)m2 (D[r2[Y[t], Phi[t]], t] . D[r2[Y[t], Phi[t]], t])  /. {Y[t]-> y, Y'[t]-> yp, Phi[t]-> ϕ, Phi'[t]-> ϕp}];
V[y_, ϕ_]:= m2 g l  Cos[ϕ];

(* Lagrangiano *)
L[y_, ϕ_, yp_, ϕp_]:= T[y, ϕ, yp, ϕp] - V[y, ϕ];

Collect[L[y, ϕ, yp, ϕp], {m1, m2}, Simplify]

Ecuación de movimiento obtenida usando directamente EulerEquations:

Needs["VariationalMethods`"]

(* Ahora son dos ecuaciones *)
ecEE = EulerEquations[L[y[t], ϕ[t],y'[t],  ϕ'[t]], {y[t], ϕ[t]}, t][[All, 1]]

(* Las simplificamos como antes,  definiendo adicionalmente l = 1. *)
ecEEfinal = Simplify[(1/m2)(ecEE /.  {y[t]-> y, ϕ[t]-> ϕ, y'[t]-> y', ϕ'[t]-> ϕ', y''[t]-> y'', ϕ''[t]-> ϕ'', g-> l (ω^2), m1-> μ m2})/. {ω-> 1, l-> 1}]

(* Despejando y'' de la primera ecuación y reemplazando en la segunda, se obtiene una ecuación sólo para ϕ. *)

ypp = Solve[ecEEfinal[[1]] == 0, y''][[1]][[1]][[2]];
ecϕ =Collect[ecEEfinal[[2]] /. y'' -> ypp, ϕ'', FullSimplify]

Identidades trigonométricas mediante, con μ = m1, es la misma ecuación que usamos al resolver numéricamente el problema.