Quit[];

Mecánica Clásica 2019 2c - Guía 2, problema 1.

Notebook 'mínimo local' para el problema del péndulo con soporte móvil.

El sistema tiene sólo dos posibles movimientos cualitativamente diferentes: 

• oscilaciones periódicas, en donde ϕ avanza y retrocede,
• rotaciones periódicas, en donde ϕ es una función siempre creciente o siempre decreciente.

Todos los comentarios entre signos "(* ........*)" pueden borrarse. Están ahí como guía, pero entorpecen la visión.

Abajo: 

• la solución numérica 
• la extracción del período [detección de los sucesivos cruces de ϕ por 0]
• el gráfico en bruto de las funciones ϕ(t) e y(t)

(* Ejecutar esta celda con Shift+Enter o el Enter del teclado numérico [el cursor debe estar dentro de la celda] *)

Clear[ϕ, y]

(* Condiciones iniciales para y, ϕ y ϕ punto. Experimenten cambiando estos valores *)

y0 = 0;
ϕ0 =π/2;
ϕp0 = 0.2;

(*  Si cuentan como yo (con los dedos), verán que estas son sólo 3 condiciones. Falta una. *)
(*  La condición inicial suplementaria está implícita en las ecuaciones de movimiento dentro del NDSolve; a saber, (Subscript[P, y]) = 0. *)
(*  Es decir, el CM del sistema no se mueve en la dirección 𝓎. Recordar, además,que (Subscript[P, y]) se conserva. *)
(*  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. *)

m1 = 0.2; (* Es el único parámetro que importa. Estoy tomando m2 = 1, o, si se quiere, definiendo m1 en unidades de m2. Experimenten cambiando m1. *)

tmax = 30; (* Un tiempo elegido a ojo, lo suficientemente grande como para incluir varios períodos *)

(* Ahora sí: la solución numérica mediante NDSolve. *)
(* La conservación de (Subscript[P, y]) permitiría reducir todo a una ecuación diferencial para ϕ. *) 
(* De todas maneras, incluyo la ecuación de primer orden para y[t] así tienen un ejemplo con dos funciones incógnitas. *)
(* En realidad podríamos usar la conservación de (Subscript[P, y ])para escribir directamente y[t] = y0 -(1/(1+m1))(Sin[ϕ[t]] - Sin[ϕ0]). *)

{sol, CrucesPorCero} =

Reap[
NDSolve[
 {
(1-(((Cos[ϕ[t]])^2)/(1+m1))) ϕ''[t] +(( Sin[2ϕ[t]] ((ϕ'[t])^2))/(2(1+m1)))+ Sin[ϕ[t]]== 0, (* Ecuación para ϕ. ¡Deducir y revisar!*)
y'[t] ==-(1/(1+m1))Cos[ϕ[t]] ϕ'[t], (* Ecuación para y, habiendo tomado que el impulso total (Subscript[P, y]) = 0 *)

ϕ[0] == ϕ0,ϕ'[0] == ϕp0, y[0] == y0, (* Condiciones iniciales*)

WhenEvent[Sin[ϕ[t]/2] == 0 ,Sow[{t, Sign[ϕ'[t]]}]] (* Guarda los tiempos y la dirección de movimiento cuando la barra cruza el eje x *)

}, {ϕ, y}, {t, 0, tmax}  (* Funciones incógnitas e intervalo de integración *)
]
];


cruces = CrucesPorCero[[1]][[All, 1]]; (* Esta lista guarda los tiempos en los que la barra cruza por ϕ = 0 *)
sentido = CrucesPorCero[[1]][[All, 2]]; (* Esta lista guarda el sentido en que se mueve la barra en cada paso por ϕ = 0 *)

Print["Phi = ",Phi = sol[[1]][[1]][[2]], "\n"]
Print["Y = ", Y = sol[[1]][[2]][[2]], "\n"]
Print["Cruces por cero  -> ",  cruces, "\n"]
Print["Δt entre cruces consecutivos -> ", Differences[cruces], "\n"]
Print["El tipo de movimiento es de ", If[Length[Tally[sentido]] == 2, "oscilación.\n", "rotación.\n"]]

(*

Respecto al cálculo del período hecho aquí abajo: 

• Si el sistema oscila como un péndulo, hay cruces en uno y otro sentido. Por eso el período es 2 veces el intervalo entre cruces por 0. ¿Por qué?

• Si el sistema rota siempre en un mismo sentido, el período es directamente el intervalo entre sucesivos cruces por cero.

• La función Tally (¡F1 Tally!) lo que me dice es cuántos tipos de cruces hay. Si todos son en el mismo sentido, entonces Tally va a tener un sólo elemento. Si hay cruces en uno y otro sentido Tally va a tener 2 elementos. Entonces si multiplico por la longitud de Tally obtengo el período correcto, tanto para las oscilaciones como para las rotaciones.

*)

Print["Periodo T ≃ ", T = Differences[cruces][[1]]× Length[Tally[sentido]]]  

Plot[{Phi[t], Y[t]}, {t, 0,2T},PlotLegends->"Expressions", PlotStyle->Thick, ImageMargins->50, PlotLabel-> Style["ϕ(t), y(t) \n", 20]]

Abajo, la animación:

Ncuadros = 100;
Δ = T/Ncuadros;  (* Paso temporal entre configuraciones sucesivas. Habrá, en principio, Ncuadros gráficos distintos, desde t = 0 hasta t = T - Δ, a intervalos de Δ.*)

Manipulate[ (* La variable "t" se controla mediante este Manipulate. *)

f = Y[t]; 
g = Phi[t]; (* Evalúo las soluciones numéricas en t. *)

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

Show[
ejes,
orb,

Graphics[{AbsoluteThickness[2], Blue, Line[{{f, 0}, {f + Sin[g], -Cos[g]}}]}],
Graphics[{Green, Disk[ {f + Sin[g], -Cos[g]}, 0.1]}],
Graphics[{Red, Disk[{f, 0}, 0.1]}],

AspectRatio-> Automatic,
PlotRange-> 2{{-1, 1}, {-1, 1}},
Axes-> False,
Ticks-> False,
Background-> Black,
ImageSize-> 600
],
{t, 0, T-Δ, Δ, Appearance-> "Open"},(* El tiempo máximo llega hasta T - Δ, porque si incluyera T la animación repetiría el primer cuadro, ya que en t = 0 y t = T la configuración es la misma. El paso temporal por defecto es Δ. *)

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