Máquina de Atwood con dos grados de libertad.

Para los que nunca hayan usado el Mathematica, lo indispensable que necesitan saber es que para ejecutar cada grupo de comandos deben usar las teclas shift+enter, o el enter del teclado numérico.  Prueben para empezar por la celda de abajo.

unafuncion[x_]:= Sin[x];
Plot[unafuncion[x], {x, 0, 4π}]

Las ecuaciones:

El sistema:

El lagrangiano:

L= (1/2)(m + M) ((OverDot[r])^2)+ (1/2)m (r^2)((OverDot[θ])^2)+ g r (m Cos[θ] - M)

En lugar de la variable t es útil cambiar a una variable adimensional τ= (Sqrt[g/(Subscript[r, 0])])t, donde (Subscript[r, 0]) es el radio inicial. Si además se escribe m = ϵ M, queda

L =M g (Subscript[r, 0]) [(1/2)(1 + ϵ) ((OverDot[x])^2)+ (1/2)ϵ (x^2)((OverDot[θ])^2)+  x (ϵ Cos[θ] - 1)],

donde x = r/(Subscript[r, 0]), y las derivadas están tomadas respecto de τ, no de t. La ventaja de escribir el lagrangiano de ese modo es que 'universaliza' el problema. Las condiciones iniciales serán siempre x(0) = 1,  θ(0) = (π/2),  (OverDot[x])(0) = (OverDot[θ])(0) = 0.  Cada problema en particular se recupera rescaleando la variable temporal. 

Multiplicar el lagrangiano por una constante tampoco afecta las ecuaciones de movimiento, de manera que podemos concentrarnos en estudiar el siguiente lagrangiano 

L =(1/2)(1 + ϵ) ((OverDot[x])^2)+ (1/2)ϵ (x^2)((OverDot[θ])^2)+  x (ϵ Cos[θ] - 1).

La energía se conserva, de modo que podemos escribir  W = (1/2)(1 + ϵ) ((OverDot[x])^2)+ (1/2)ϵ (x^2)((OverDot[θ])^2)-  x (ϵ Cos[θ] - 1) = constante.  Notar que para las condiciones iniciales W = 1.

Las ecuaciones de movimiento son 

(1 + ϵ)  (Overscript[x, ..]) = ϵ x ((OverDot[θ])^2) + ϵ Cos[θ] - 1,
x(Overscript[θ, ..]) + 2 (OverDot[x]) (OverDot[θ]) = - Sin[θ]. 

Este sistema de ecuaciones no puede resolverse de manera sencilla. Más abajo figura la solución numérica.

A la derecha de la cada recuadro, abajo, aparece una pequeña barra con una flecha. Usar doble click para expandir o colapsar. 

El último recuadro es autosuficiente y puede ejecutarse sin haber leído ninguno de los anteriores.

La solución numérica en el caso ϵ = 1/3.

Se fija el valor de ϵ y el tiempo máximo para la solución, llamado tmax. La solución se guarda en las funciones rsol y θsol. Además se define la función energía con el símbolo W. Como W es una cantidad conservada, graficar W en función del tiempo es un modo de comprobar si la integración numérica deja de ser precisa en algún momento. 

En el ejemplo de abajo ϵ = 1/3 y tmax =300.  Variando esos dos parámetros más o menos controlan todo.  Si tmax es muy alto, se corre el riesgo de que el Mathematica avise que la integración falla, pero por lo menos avisa.

ϵ = 1/3;
tmax =300;

solucion = Flatten[NDSolve[{r''[t] == (ϵ/(1+ϵ)) r[t] ((θ'[t])^2)+  ((ϵ/(1+ϵ)) Cos[θ[t]]-(1/(1+ϵ))),r[t] θ''[t] + 2 r'[t]θ'[t] == -   Sin[θ[t]], r[0] == 1, r'[0] == 0, θ[0] ==  (π/2), θ'[0] == 0}, {r, θ}, {t, 0, tmax}, MaxSteps-> 100000]];

rsol[t_]:=Evaluate[r[t]/.solucion];
θsol[t_]:=Evaluate[θ[t]/.solucion];

W[t_]:= (1/2)(1+ϵ) ((rsol'[t])^2)+ (1/2)ϵ ((rsol[t])^2) ((θsol'[t])^2)-  rsol[t] (ϵ Cos[θsol[t]] - 1);

Verificaciones

La diferencia entre la  energía en función del tiempo y su valor inicial. Idealmente debería ser nula. (Como los valores que se quieren graficar son muy pequeños y cubren varios órdenes de magnitud, el gráfico está en escala logarítmica.)

Plot[((Log[Abs[W[t]-1]])/(Log[10])), {t,0,tmax}, PlotRange-> All]

Notar que a medida que pasa el tiempo, se va deteriorando la precisión de la integración numérica, aunque dentro del intervalo de tiempo estudiado siempre es muy pequeña la variación de W.

Detalle del gráfico anterior para tiempos entre 140 y 150.

Plot[((Log[Abs[W[t]-1]])/(Log[10])), {t,140,150}, PlotRange-> All]

También pueden graficarse las propias ecuaciones de movimiento. Por ejemplo

Plot[rsol''[t] -( (ϵ/(1+ϵ)) rsol[t] ((θsol'[t])^2)+  ((ϵ/(1+ϵ)) Cos[θsol[t]]-(1/(1+ϵ)))), {t, 0, tmax}, PlotRange-> {-(10^(-5)), (10^(-5))}]

Pueden incluirse algunos comandos para aumentar la precisión del cálculo, al costo de hacer todo un poco más lento.

wp = 25;

solucionWP = Flatten[NDSolve[{r''[t] == (ϵ/(1+ϵ)) r[t] ((θ'[t])^2)+  ((ϵ/(1+ϵ)) Cos[θ[t]]-(1/(1+ϵ))),r[t] θ''[t] + 2 r'[t]θ'[t] == -   Sin[θ[t]], r[0] == 1, r'[0] == 0, θ[0] ==  (π/2), θ'[0] == 0}, {r, θ}, {t, 0, tmax}, MaxSteps-> 100000, AccuracyGoal->Automatic,PrecisionGoal->Automatic, WorkingPrecision-> wp]];

rsolWP[t_]:=Evaluate[r[t]/.solucionWP];
θsolWP[t_]:=Evaluate[θ[t]/.solucionWP];

WWP[t_]:= (1/2)(1+ϵ) ((rsolWP'[t])^2)+ (1/2)ϵ ((rsolWP[t])^2) ((θsolWP'[t])^2)-  rsolWP[t] (ϵ Cos[θsolWP[t]] - 1);

Abajo se grafican las energías calculadas por los dos métodos.  (Toma su tiempo.) El error disminuye unos 4 órdenes de magnitud.

Plot[{((Log[Abs[WWP[t]-1]])/(Log[10])),((Log[Abs[W[t]-1]])/(Log[10]))}, {t,0,tmax}, PlotRange-> All]

Gráficos de las coordenadas en función del tiempo.

Pueden graficar θ(t) y r(t) y sus derivadas. Abajo, el gráfico de θ(t) y el de dθ/dt.

Plot[θsol[t], {t, 0, tmax/3}, Axes-> False, Frame-> True]

Plot[θsol'[t], {t,0,tmax/2}, PlotRange-> All, Axes-> False, Frame-> True]

Abajo, los gráficos de r, dr/dt .

Plot[rsol[t], {t,0,tmax/2}, PlotRange-> All,Axes-> False, Frame-> True]

Plot[rsol'[t], {t,0,tmax/2}, PlotRange-> All, Axes-> False, Frame-> True]

Es interesante graficar la trayectoria en el diagrama de fases asociado a cada coordenada. Esto es, gráficos de r y (Subscript[p, r])= (OverDot[r]) , por un lado, y de θ y (Subscript[p, θ])= (r^2)(OverDot[θ]), por otro.

Abajo, la trayectoria en el subespacio definido por   r y (OverDot[r])

fasegrafr = ParametricPlot[{rsol[t], rsol'[t]}, {t, 0, tmax}, AspectRatio-> Automatic, Axes-> False, Frame-> True, PlotPoints-> 100, PlotStyle-> {Opacity[0.5],RGBColor[0,0,.5]}, ImageSize-> 300]

Abajo, la trayectoria en el subespacio definido por   θ y (r^2)(OverDot[θ]).

fasegrafθ = ParametricPlot[{θsol[t],((rsol[t])^2) θsol'[t]}, {t, 0, tmax}, PlotRange-> All, AspectRatio->1, Axes-> False, Frame-> True, PlotPoints-> 200, PlotStyle-> {Opacity[0.5],RGBColor[.5,0,0]}, ImageSize-> 300]

Abajo, superpuesta a la trayectoria, la animación muestra el punto que ocupa el sistema instante a instante en el espacio definido por  r y (OverDot[r]) y por θ y (OverDot[θ]).

Animate[{Show[fasegrafr, Graphics[{Orange, Disk[{rsol[t], rsol'[t]}, 1/40]}], ImageSize-> 300],Show[fasegrafθ, Graphics[{Orange, Disk[{θsol[t], ((rsol[t])^2)θsol'[t]}, {1/5, 1/20}]}], PlotRange-> {{-8.2,2}, {-1.5, 1.5}}, ImageSize-> 420]}, {t, 0, tmax, 1/100}, AnimationRunning-> False, AnimationRate-> 500]

Una variación del gráfico anterior, que muestra la parte de la trayectoria en el espacio de fases ya recorrida. Ejecutar antes la primera celda a continuación.

recorridor = Table[{rsol[t], rsol'[t]}, {t, 0, tmax, 1/100}];
recorridoθ = Table[{θsol[t], ((rsol[t])^2)θsol'[t]}, {t, 0, tmax, 1/100}];

Animate[{
Show[ListPlot[Take[recorridor, i], Joined-> True, PlotStyle-> {RGBColor[1,0,0], Thickness[0.006]}, Axes-> False, Frame-> True],fasegrafr, Graphics[{Orange, Disk[ {rsol[(i/100)], rsol'[(i/100)]}, 1/40]}], AspectRatio-> Automatic,PlotRange-> {{0.1,1.5}, {-1.0, 1.0}}, ImageSize-> 300],Show[ListPlot[Take[recorridoθ, i], Joined-> True, PlotStyle-> {RGBColor[1,0,0], Thickness[0.006]}, Axes-> False, Frame-> True],fasegrafθ, Graphics[{Orange, Disk[ {θsol[(i/100)],((rsol[(i/100)])^2)θsol'[(i/100)]}, {1/5, 1/20}]}], AspectRatio-> 1,PlotRange-> {{-8.2,2}, {-1.5, 1.5}},ImageSize-> 420]}, {i, 1,100 tmax, 1},AnimationRate->500, AnimationRunning-> False]

La trayectoria de la masa oscilante puede graficarse en el plano xy.  La envolvente corresponde a la curva en donde dr/dt = dθ/dt = 0, y cuya ecuación (pueden demostrar) es r = 1/(1-ϵ Cos[θ]).

x[t_]:= rsol[t] Sin[θsol[t]];
y[t_]:= -rsol[t] Cos[θsol[t]];

trayectoria = Show[ParametricPlot[{{x[t], y[t]}},{t, 0,tmax}, PlotPoints-> 80], ParametricPlot[(1/(1-(1/3)Cos[θ])){Sin[θ], -Cos[θ]}, {θ, 0, 2π}, PlotStyle-> RGBColor[1,0,0]], PlotRange-> {{-1.2, 1.2}, {-1.6, .8}}]

El gráfico anterior muestra la trayectoria, pero no permite distinguir la evolución temporal.

Animaciones:

El siguiente gráfico  muestra la evolcuión del sistema. Para graficar la posición de la masa más grande se supuso que la longitud de la cuerda es L =  2. [La estela que deja la masa pequeña puede borrarse, es el primer ParametricPlor dentro del primer Show.]

L = 2;

Animate[{Show[ParametricPlot[{x[τ], y[τ]},{τ, Max[(10^(-4)), t - (tmax/100)], t},PlotRange-> {{-1.2, 1.2}, {-1.6, .8}}, Axes-> False], Graphics[{Line[{{0,0}, {0,-L+rsol[t]}}], {Gray,Disk[{0, -L+rsol[t]}, 1/12]}, Line[{{0,0}, {x[t],y[t]}}], {Orange,Disk[{x[t], y[t]}, 1/20]}}], PlotRange-> {{-1.2, 1.2}, {-2.2, .62}}, ImageSize-> 400],Show[trayectoria, Graphics[{Orange, Disk[{x[t], y[t]}, 1/25]}], ImageSize-> {400,400}]}, {{t, 0.,"tiempo"},0, tmax, (1./500)},   AnimationRunning -> False,  AnimationRate->1]

Caso M/m = 100.

L = 2;
ϵ = 1/100;
tmax =200;

solucion = Flatten[NDSolve[{r''[t] == (ϵ/(1+ϵ)) r[t] ((θ'[t])^2)+  ((ϵ/(1+ϵ)) Cos[θ[t]]-(1/(1+ϵ))),r[t] θ''[t] + 2 r'[t]θ'[t] == -   Sin[θ[t]], r[0] == 1, r'[0] == 0, θ[0] ==  (π/2), θ'[0] == 0}, {r, θ}, {t, 0, tmax}, MaxSteps->100000]];

rsol100[t_]:=Evaluate[r[t]/.solucion];
θsol100[t_]:=Evaluate[θ[t]/.solucion];

x100[t_]:= rsol100[t] Sin[θsol100[t]];
y100[t_]:= -rsol100[t] Cos[θsol100[t]];

W100[t_]:= (1/2)(1+ϵ) ((rsol100'[t])^2)+ (1/2)ϵ ((rsol100[t])^2) ((θsol100'[t])^2)-  rsol100[t] (ϵ Cos[θsol100[t]] - 1);

Gráfico de la energía. (Pierde precisión mucho más rápidamente que en el caso ϵ = 1/3.)

Plot[W100[t]-1, {t,0,tmax}, PlotRange-> All]

Los gráficos de θ y (OverDot[θ]) muestran variaciones mucho más bruscas; intervalos de comportamiento suave separados por saltos repentinos.

Plot[θsol100[t], {t, 0, tmax/4},Axes-> False, Frame-> True]

Plot[θsol100'[t], {t, 0, tmax/3}, PlotRange-> All,Axes-> False, Frame-> True]

Los gráficos en el espacio de fases:

ParametricPlot[{rsol100[t], rsol100'[t]}, {t, 0, tmax}, AspectRatio-> 1, Axes-> False, Frame-> True, PlotPoints-> 100, PlotStyle-> {Opacity[0.5],RGBColor[0,0,.5]}, ImageSize-> 400]

ParametricPlot[{θsol100[t],((rsol100[t])^2) θsol100'[t]}, {t, 0, tmax}, AspectRatio-> 1, Axes-> False, Frame-> True, PlotPoints-> 100, PlotStyle-> {Opacity[0.5],RGBColor[0,0,.5]}, ImageSize-> 400]

Gráfico de la trayectoria

trayectoria100 = Show[ParametricPlot[{x100[t], y100[t]},{t, 0,tmax}, PlotPoints-> 500], ParametricPlot[(1/(1-(1/100)Cos[θ])){Sin[θ], -Cos[θ]}, {θ, 0, 2π}, PlotStyle-> RGBColor[1,0,0]], PlotRange-> {{-1.2, 1.2}, {-1.1, 1}}]

Animación del movimiento de las dos masas.

Animate[{Show[ParametricPlot[{x100[τ], y100[τ]},{τ, Max[(10^(-4)), t - (tmax/400)], t},PlotRange-> {{-1.2, 1.2}, {-2.2, 1}}, Axes-> False], Graphics[{Line[{{0,0}, {0,-L+rsol100[t]}}], {Gray,Disk[{0, -L+rsol100[t]}, 1/12]}, Line[{{0,0}, {x100[t],y100[t]}}], {Orange,Disk[{x100[t], y100[t]}, 1/20]}}], PlotRange-> {{-1.2, 1.2}, {-2.2, 1}}, ImageSize-> 400],Show[trayectoria100, Graphics[{Orange, Disk[{x100[t], y100[t]}, 1/25]}], ImageSize-> {400,400}]}, {t, 0., tmax, (1./500)},AnimationRunning-> False,  AnimationRate->1]

La comparación entre los casos ϵ = 1/100 y ϵ = 1/3.

Animate[{Show[ParametricPlot[{x100[τ], y100[τ]},{τ, Max[(10^(-4)), t - (tmax/400)], t},PlotRange-> {{-1.2, 1.2}, {-2.2, .6}}, Axes-> False], Graphics[{Line[{{0,0}, {0,-L+rsol100[t]}}], {Gray,Disk[{0, -L+rsol100[t]}, 1/12]}, Line[{{0,0}, {x100[t],y100[t]}}], {Orange,Disk[{x100[t], y100[t]}, 1/50]}}], PlotRange-> {{-1.2, 1.2}, {-2.2, .62}}, ImageSize-> 250],Show[ParametricPlot[{x[τ], y[τ]},{τ, Max[(10^(-4)), t - (tmax/400)], t},PlotRange-> {{-1.2, 1.2}, {-2.2, .6}}, Axes-> False],
 Graphics[{Line[{{0,0}, {0,-L+rsol[t]}}], {Gray,Disk[{0, -L+rsol[t]}, 1/12]}, Line[{{0,0}, {x[t],y[t]}}], {Orange,Disk[{x[t], y[t]}, 1/20]}}], PlotRange-> {{-1.2, 1.2}, {-2.2, .62}}, ImageSize-> 250]}, {t, 0., tmax, (1./500)},AnimationRunning-> False, AnimationRate->3]

Un ejemplo de cómo guardar los resultados en archivos gif o avi.

Manipulate[{Show[ParametricPlot[{x100[τ], y100[τ]},{τ, Max[(10^(-4)), t - (tmax/400)], t},PlotRange-> {{-1.2, 1.2}, {-2.2, .6}}, Axes-> False], Graphics[{Line[{{0,0}, {0,-L+rsol100[t]}}], {Gray,Disk[{0, -L+rsol100[t]}, 1/12]}, Line[{{0,0}, {x100[t],y100[t]}}], {Orange,Disk[{x100[t], y100[t]}, 1/20]}}], PlotRange-> {{-1.2, 1.2}, {-2.2, .62}}, ImageSize-> 250],Show[trayectoria100, Graphics[{Orange, Disk[{x100[t], y100[t]}, 1/25]}], ImageSize-> {250,250}]}, {t, 0., tmax, (1./500)}]

gtable100 = Table[Show[ParametricPlot[{x100[τ]-3, y100[τ]},{τ, Max[(10^(-4)), t - (tmax/400)], t}, Axes-> False], Graphics[{Line[{{-3,0}, {-3,-L+rsol100[t]}}], {Gray,Disk[{-3, -L+rsol100[t]}, 1/12]}, Line[{{-3,0}, {x100[t]-3,y100[t]}}], {Orange,Disk[{x100[t]-3, y100[t]}, 1/20]}}], trayectoria100, Graphics[{Orange, Disk[{x100[t], y100[t]}, 1/25]}], PlotRange-> {{-4.2,1.2},{-2.2,.62}}, ImageSize-> 600], {t, 0,11.5, 1/25}];

Esto exporta la tabla de gráficos como una película avi.

Export["masitas.avi", gtable100, "AVI"]

Esto exporta la tabla de gráficos como una animación gif.

Export["masitas.gif", gtable100, "GIF"]

El directorio dónde se guardan los archivos puede consultarse con Directory[].

Directory[]

El directorio puede modificarse con SetDirectory: SetDirectory["ruta"].

Un recorrido visual por las trayectorias a medida que se varía la relación entre las dos masas. Aquí ϵ = 1/η.

Primero se construye una tabla con todas las soluciones, con η variando entre 1.3 y 10. [El comando Monitor permite ver el grado de avance del cálculo; cuando llega a cero termina.]

tmax =150;
Monitor[tsolucion = Table[Flatten[NDSolve[{r''[t] == ((1/η)/(1+1/η)) r[t] ((θ'[t])^2)+  (((1/η)/(1+1/η)) Cos[θ[t]]-(1/(1+1/η))),r[t] θ''[t] + 2 r'[t]θ'[t] == -  Sin[θ[t]], r[0] == 1, r'[0] == 0, θ[0] ==  (π/2), θ'[0] == 0} , {r, θ}, {t, 0, tmax}, MaxSteps->100000]], {η, (13/10), 10,(1/100)}], 10-N[η]];

Module[{i}, For[i = 1, i ≤ Length[tsolucion],trsol[t_, i]:= Evaluate[r[t]/.tsolucion[[i]]];++i];
For[i = 1, i ≤ Length[tsolucion],tθsol[t_, i]:= Evaluate[θ[t]/.tsolucion[[i]]];++i]];

Ahora se construye una tabla con los gráficos de todas las trayectorias. [Lleva memoria y tiempo, pero sólo hay que hacerlo una vez.]

tgraftray = Monitor[Table[Show[ParametricPlot[{trsol[t,i] Sin[tθsol[t,i]],-trsol[t,i] Cos[tθsol[t,i]]} , {t, 0, 100}],Graphics[Text[Style[N[(13/10)+ ((i-1)/100), 3], Large],{1.1,0.7}]], PlotRange-> {{-1.6, 1.6}, {-2, 1}}], {i, 1,Length[tsolucion], 1}], Length[tsolucion]-i];

Para no tener que calcular todo de nuevo, se guardan todos los gráficos en un archivo.

tgraftray >> familia_de_trayectorias

Ejemplo para recuperar el archivo:

tgraftray  = << familia_de_trayectorias;

El directorio dónde se guardan los archivos puede consultarse con Directory[].

Directory[]

El gráfico de todas las trayectorias. Cada trayectoria está encerrada por la curva en donde dr/dt = dθ/dt = 0.

Manipulate[Show[ParametricPlot[(1/(1-(1/((13/10)+ ((i-1)/100)))Cos[θ])){Sin[θ], -Cos[θ]}, {θ, 0, 2π}, PlotStyle-> RGBColor[1,0,0]],tgraftray[[i]], PlotRange-> {{-1.6, 1.6}, {-1.5, 1}}], {i, 1, Length[tgraftray],1}]

Se pueden estudiar los casos interesantes refinando el valor de η; por ejemplo, para ver lo que pasa en la vecindad de η = 4.75. 

La variable "paso" abajo gobierna la escala en que se varía η, partiendo desde el valor 4.75 y con incrementos de ±(10^(-paso)).

L = 2;
tmax = 100;
Manipulate[
Module[
{soln,theta, radio, x, y,t},soln = Flatten[NDSolve[{r''[t] == ((1/η)/(1+1/η)) r[t] ((θ'[t])^2)+  (((1/η)/(1+1/η)) Cos[θ[t]]-(1/(1+1/η))),r[t] θ''[t] + 2 r'[t]θ'[t] == -   Sin[θ[t]], r[0] == 1, r'[0] == 0, θ[0] ==  (π/2), θ'[0] == 0} , {r, θ}, {t, 0, tmax}, MaxSteps->100000]];
theta[t_]:=Evaluate[θ[t]/.soln];
radio[t_]:= Evaluate[r[t] /. soln];

ParametricPlot[{radio[t] Sin[theta[t]],-radio[t] Cos[theta[t]]} , {t, 0, tmax}, PlotRange-> {{-1.2, 1.2}, {-1.6, .62}}, ImageSize-> 550, PerformanceGoal-> "Quality"]
],
{{paso, 2}, 1, 10,1,Appearance->"Open"},
{{η,4.75,"relación entre las masas"},2,10,(10^(-paso)),Appearance->{"Open"}},SynchronousUpdating->False]

Animación que puede ejecutarse sin haber leído nada de lo anterior.

La celda de abajo permite variar la relación M/m = η entre las masas sin interrumpir la simulación. La condición inicial se actualiza cada vez que se cambia el valor de η, de modo que la trayectoria es continua, pero en general  no corresponderá a la misma condición en t= 0 para distintos valores de η. 

Para ejecutar, situarse sobre cualquier punto de la celda y apretar shift+enter o el enter del teclado numérico.

DynamicModule[{L = 2, tmax = 200, radio = 1,  theta = π/2, radioprima = 0, thetaprima = 0, ηinicio = 3, cero = 0,   run, x, y,  lastη, soln, r, θ, trayectoria, tru = False, reset=True},

Manipulate[

If[η ≠ lastη ||  reset , soln = Flatten[NDSolve[{r''[t] == ((1/η)/(1+1/η)) r[t] ((θ'[t])^2)+  (((1/η)/(1+1/η)) Cos[θ[t]]-(1/(1+1/η))),r[t] θ''[t] + 2 r'[t]θ'[t] == -   Sin[θ[t]], r[run -cero] == radio, r'[run -cero] == radioprima, θ[run -cero] ==  theta, θ'[run -cero] == thetaprima} , {r, θ}, {t, 0, tmax}, MaxSteps->100000]];

trayectoria = ParametricPlot[{ (radio =r[t]/.soln) Sin[(theta = θ[t]/.soln)], -radio  Cos[theta]}, {t, 0, tmax},PlotPoints-> 200, PerformanceGoal-> "Quality", PlotStyle-> Opacity[0.5]];
lastη = η;
reset = False; ];

If[tr ==  True  && tru ≠ tr, trayectoria = ParametricPlot[{ (radio =r[t]/.soln) Sin[(theta = θ[t]/.soln)], -radio  Cos[theta]}, {t, 0, tmax},PlotPoints-> 200, PerformanceGoal-> "Quality", PlotStyle-> Opacity[0.5]]; tru = True];

theta=θ[run -cero]/.soln;
radio = r[run -cero] /. soln;
thetaprima=θ'[run -cero]/.soln;
radioprima = r'[run -cero] /. soln;
x = radio Sin[theta];
y = -radio Cos[theta];
Switch[setη, 1, η = 1.66500000;  radio = 1;  theta = π/2; radioprima = 0; thetaprima = 0;  cero = 0; reset = True, 
2, η = 1.55500000; radio = 1;  theta = π/2; radioprima = 0; thetaprima = 0;  cero = 0; reset = True,
3,  η = 4.17700000; radio = 1;  theta = π/2; radioprima = 0; thetaprima = 0;  cero = 0; reset = True,
 4, η = 7.71000000 ; radio = 1;  theta = π/2; radioprima = 0; thetaprima = 0;  cero = 0; reset = True];
setη = Null;

If[run == tmax, cero = 0; Abort[]];

Show[ Graphics[{If[tr, trayectoria[[1]], {}], {Line[{{0,0}, {0,-L+radio}}], {Gray,Disk[{0, -L+radio}, 1/12]}, Line[{{0,0}, {x,y}}], {Orange,Disk[{x, y}, 1/(12(Sqrt[η]))+(1/50)]}}}], PlotRange-> {{-1.2, 1.2}, {-2.2, 1.5}}, ImageSize-> 300],

Style["% Animación", Bold],
Dynamic[N[100 run/tmax]],
{{run,0, ""},0,tmax, ControlType->Trigger,   AnimationRate-> 1, AppearanceElements->{"PlayPauseButton","FasterSlowerButtons"}},
Style["Al llegar a 100% se detiene"],

Delimiter,
{{η,ηinicio,"η = M/m"},1.1,100,(10^(-paso)), Appearance->"Open", AppearanceElements->{"InputField", "StepLeftButton", "StepRightButton"}},
Style["Algunos valores interesantes de η", FontSize-> 10],
{{setη,Null, ""},{ 1-> "1.665", 2->"1.555", 3-> "4.177", 4-> "7.710"},ControlType->Setter, Enabled-> True},


Delimiter,
{{paso, 1, "Δη =  ± (10^(-x))" }, 0, 10, 1, Appearance -> "Open",  AppearanceElements->{"InputField", "StepLeftButton", "StepRightButton"}},
Style["Δη controla la presición en el ajuste de η"],

Delimiter,
{{tr, False, Style["Mostrar trayectoria", Bold]}, {True, False}},

Delimiter,
Style["Resetear a la condición inicial por defecto",Bold],
Button[Graphics[{Line[{{0,0}, {0,-L+1}}], {Gray,Disk[{0, -L+1}, 1/12]}, Line[{{0,0}, {1,0}}], {Orange,Disk[{1, 0},  1/(12(Sqrt[η]))+(1/50)]}}], 
cero = run;
radio = 1;
theta = π/2;
radioprima = 0;
thetaprima = 0;
reset = True,
ImageSize->75],

 SynchronousUpdating->False, ControlPlacement->Left, SaveDefinitions-> False]]