Quit[];

(* para la salida de gráficos *)
antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]];

Péndulo doble: λ = l1/l2, μ = (m2/(m1 + m2)).

(*parámetros del péndulo. μ tiene que ser menor que 1*)
λ = 1/2;
μ = 2/3; 

(* condiciones iniciales *)
ϕ10 = 3;
ϕ20 = 2;
ϕ1p0 = 0;
ϕ2p0 = 0;

(*Este valor es para integrar las ecuaciones un poco en el pasado*)
tpast = -10;

(*Integrador numérico de las ecuaciones*)
sol = NDSolve[
{

(*ecuaciones diferenciales*)
λ ϕ1''[t] + μ (ϕ2''[t] Cos[ϕ2[t] - ϕ1[t]] - ((ϕ2'[t])^2)Sin[ϕ2[t] - ϕ1[t]]) +Sin[ϕ1[t]] == 0, 
ϕ2''[t] + λ (ϕ1''[t] Cos[ϕ2[t] - ϕ1[t]] + ((ϕ1'[t])^2)Sin[ϕ2[t] - ϕ1[t]]) +Sin[ϕ2[t]] == 0,

(*condiciones iniciales*)
ϕ1[0] ==  ϕ10,
ϕ2[0] == ϕ20,
ϕ1'[0] ==ϕ1p0,
ϕ2'[0] == ϕ2p0
},

{ϕ1, ϕ2}, {t, tpast, tmax = 50}

(* Se puede aumentar la precisión del cálculo descomentando la siguiente línea. Ajusten los valores como mejor les convenga. Esto se nota si grafican la conservación de la energía*)
(*WorkingPrecision-> 20, MaxSteps-> (10^6), *)

][[1]];

(*Extrae cada ángulo de la solución*)
{ϕ1sol, ϕ2sol} = sol[[All,2]];

(*contruye las funciones posición de cada masa*)
r1sol[t_]:= λ {Sin[ϕ1sol[t]], -Cos[ϕ1sol[t]]} ;
r2sol[t_]:= r1sol[t] + {Sin[ϕ2sol[t]], -Cos[ϕ2sol[t]]};

(*energía mecánica*)
e[t_]:= (1/2)(λ^2)((ϕ1sol'[t])^2) + (1/2)μ(((ϕ2sol'[t])^2)+ 2 λ ϕ1sol'[t] ϕ2sol'[t] Cos[ϕ2sol[t] - ϕ1sol[t]]) -  λ Cos[ϕ1sol[t]] -  μ Cos[ϕ2sol[t]]

(*un gráfico de la posición de la segunda masa*)
ParametricPlot[r2sol[tt], {tt, 0, tmax}, ColorFunction->Function[{x, y, u},Directive[Opacity[(Sqrt[u])]]]]

(*verificación de que la energía se conserva - bueno, más o menos. Descomenten las líneas del WorkingPresicion para mejorar la precisión.*)
e0 = e[0];
Plot[((e[t])/e0)-1, {t, 0, tmax}, PlotRange-> All, ImageSize->600]

(*Animación*)
DynamicModule[
(*Para que no se mezclen con variables dinámicas de otros Manipulate, las siguientes variables se definen localmente dentro de este "DynamicModule"*)
{ϕ1pr, ϕ2pr, r0, r1, r2, traza, tpr},

r0 = {0,0};

Manipulate[
(*La variable de control de este bloque "Manipulate" es el tiempo tpr*)

(*Evalúa los ángulos y las posiciones a tiempo tpr*)
ϕ1pr = ϕ1sol[tpr];
ϕ2pr = ϕ2sol[tpr];

r1 = λ{Sin[ϕ1pr],-Cos[ϕ1pr] };
r2 = r1 + { Sin[ϕ2pr],-Cos[ϕ2pr]};

(*Grafica parte de la trayectoria de la segunda masa, para mostrar la traza que va dejando*)
traza = ParametricPlot[r2sol[tt], {tt, tpr-4, tpr}, ColorFunction->Function[{x, y, u},Directive[Opacity[(Sqrt[u])]]], PlotPoints-> 50];

(*"Show" reúne en un mismo gráfico varios elementos*)
Show[
(* primero muestra la traza generada más arriba*)
traza,

(* después grafica las barras y las masas*)
Graphics[
{
AbsoluteThickness[3], 
Orange, Line[{r0, r1}],
Green, Line[{r1, r2}],

White, Disk[r1, radio = 0.05],
Red, Disk[r0, (2/3)radio],
Cyan, Disk[r2, radio]
}
],

(*Se da espacio suficiente para que las masas no se salgan de la zona graficada. En estas variables adimensionales, la longitud de la primera barra es λ y la de la segunda es 1.*)
PlotRange-> 1.2{{-1-λ, 1+λ}, {-1-λ, 1+λ}},
Axes-> False, 
Background-> Black, 
ImageSize-> 500
], 

(*Esto controla el Manipulate*)
{tpr, 0, tmax, 0.01, Appearance-> "Open", AnimationRate-> 2}, 
TrackedSymbols:>{tpr}
]
]

Salida de gráficos: lo mismo de antes pero para construir una tabla con los gráficos rasterizados  vía la función “antialias” definida al comienzo del notebook.

Monitor[
tgraf = 
Table[
ϕ1pr = ϕ1sol[tpr];
ϕ2pr = ϕ2sol[tpr];
r0 = {0,0};
r1 = λ{Sin[ϕ1pr],-Cos[ϕ1pr] };
r2 = r1 + { Sin[ϕ2pr],-Cos[ϕ2pr]};
traza = ParametricPlot[r2sol[tt], {tt, tpr-5, tpr},PlotStyle-> AbsoluteThickness[2],  ColorFunction->Function[{x, y, u},Directive[Opacity[(Sqrt[u])]]]];

gtemp = antialias[
Show[
traza,
Graphics[
{
AbsoluteThickness[3], 
Orange, Line[{r0, r1}],
Green, Line[{r1, r2}],

White, Disk[r1, radio = 0.05],
Red, Disk[r0, (2/3)radio],
Cyan, Disk[r2, radio]
}
],
PlotRange-> 1.2{{-1-λ, 1+λ}, {-1-λ, 1+λ}},
Axes-> False, 
Background-> Black, 
ImageSize-> (imsize = 500)
],
2 (*este número gobierna el rasterizado (resolución); 1 deja el gráfico tal como sale, 2 lo suaviza un poco, y de ahí para arriba. Un valor muy alto hace lenta la generación de los gráficos*)
], 

(*Gradúen el paso temporal. Empiecen con un Δt grande y si ven que todo va bien, disminúyanlo*)
{tpr, 0, tmax, 0.1}] ;,

(*Esta línea pertenece al bloque  "Monitor". Va mostrando las salidas a medida que se calculan. Sirve para ver que todo marche bien.*)
{Show[gtemp, ImageSize-> imsize], tpr}
];

(*Aquí pueden animar la tabla de gráficos*)
Manipulate[tgraf[[i]], {i, 1, Length@tgraf, 1, Appearance-> "Open", AnimationRate-> 25}]

Salida de archivos

(*Crea y selecciona el directorio de salida dentro del directorio en donde se encuentra guardado el notebook*)
SetDirectory[NotebookDirectory[]];
CreateDirectory[dirsal = "salida_pendulo_doble_animacion_1"];
SetDirectory[dirsal]

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

Dependencia en las condiciones iniciales. Construye una tabla de soluciones con diferentes condiciones iniciales. Se omiten los comentarios.

Monitor[
tablasol = 
Table[
λ = 1;
μ = 1/2;

ϕ10 = 2;
ϕ20 =3+n (10^(-6));
ϕ1p0 = 0;
ϕ2p0 = 0;

sol = NDSolve[
{
λ ϕ1''[t] + μ (ϕ2''[t] Cos[ϕ2[t] - ϕ1[t]] - ((ϕ2'[t])^2)Sin[ϕ2[t] - ϕ1[t]]) + Sin[ϕ1[t]] == 0, 
ϕ2''[t] + λ (ϕ1''[t] Cos[ϕ2[t] - ϕ1[t]] + ((ϕ1'[t])^2)Sin[ϕ2[t] - ϕ1[t]]) + Sin[ϕ2[t]] == 0,

ϕ1[0] ==  ϕ10,
ϕ2[0] == ϕ20,
ϕ1'[0] ==ϕ1p0,
ϕ2'[0] == ϕ2p0
},

{ϕ1, ϕ2}, {t, tpast = -10, tmax = 100}, WorkingPrecision->20, MaxSteps->(10^6)][[1]];
sol[[All,2]],
{n, 1, 5}
];, 
n]
r2soltabla[t_]:= Table[λ {Sin[tablasol[[i]][[1]][t]], -Cos[tablasol[[i]][[1]][t]]} + {Sin[tablasol[[i]][[2]][t]], -Cos[tablasol[[i]][[2]][t]]}, {i, 1, Length@tablasol}];

etabla[t_]:= Table[ (1/2)(λ^2)((tablasol[[i,1]]'[t])^2) + (1/2)μ(((tablasol[[i,2]]'[t])^2)+ 2 λ tablasol[[i,1]]'[t] tablasol[[i,2]]'[t] Cos[tablasol[[i,2]][t] - tablasol[[i,1]][t]]) - λ Cos[tablasol[[i,1]][t]] - μ Cos[tablasol[[i,2]][t]], {i, 1, Length@tablasol}];
e0tabla = etabla[0];

Table[Plot[((etabla[t][[ipr]])/(e0tabla[[ipr]]))-1, {t, 0, tmax}], {ipr, 1, Length@tablasol}]

len = Length@tablasol;

DynamicModule[
{ϕ1tabla, ϕ2tabla, r1tabla, r2tabla, trazatabla, tpr},
Manipulate[
ϕ1tabla = Table[tablasol[[i,1]][tpr], {i, 1, len}];
ϕ2tabla = Table[tablasol[[i,2]][tpr], {i, 1, len}];
r0 = {0,0};
r1tabla = λ Table[{Sin[ϕ1tabla[[i]]],-Cos[ϕ1tabla[[i]]]}, {i, 1, len}];
r2tabla = r1tabla + Table[{Sin[ϕ2tabla[[i]]],-Cos[ϕ2tabla[[i]]]}, {i, 1, len}];

(*trazatabla = ParametricPlot[r2soltabla[tt], {tt, tpr-4, tpr}, ColorFunction->Function[{x, y, u},Directive[Opacity[(Sqrt[u])]]]];*)

gtemp = Show[

(*trazatabla,*)
Graphics[
{
AbsoluteThickness[3],
Orange, Table[Line[{r0, r1tabla[[i]]}], {i, 1, len}],
Green, Table[Line[{r1tabla[[i]], r2tabla[[i]]}], {i, 1, len}],
Red, Disk[r0, 0.03],
Blue, Table[Disk[r2tabla[[i]], 0.05 (μ/(1 - μ))],  {i, 1, len}],
White, Table[Disk[r1tabla[[i]], 0.05], {i, 1, len}]
}
],
PlotRange-> 1.2{{-1-λ, 1+λ}, {-1-λ, 1+λ}},
Axes-> False, 
Background-> Black, 
ImageSize-> 500
], 
{tpr, 0, tmax, 0.01, Appearance-> "Open", AnimationRate-> 2}
]
]