Quit[];

Integración numérica del problema 5 de la guía 4.

Esta celda hace todo el trabajo de integración y cálculo de la solución aproximada de pequeñas oscilaciones. Pueden experimentar con cilindros de otros radios.

(*Los radios de los cilindros y los parámetros (Subscript[Δ, i])*)

r0 =8;
r1 =5;
r2 =1;
Δ1 = ((r0-r1)/r2);
Δ2 = ((r1-r2)/r2);

(*Las frecuencias y autovectores del problema de pequeñas oscilaciones*)

M = {{4 ((Δ1)^2), 2 Δ1 Δ2}, {2Δ1 Δ2, 2((Δ2)^2)}};
ω1 =(Sqrt[ ((Δ1 + Δ2 + (Sqrt[((Δ1)^2)+((Δ2)^2)]))/(2Δ1 Δ2))]);
ω2 =(Sqrt[ ((Δ1 + Δ2 - (Sqrt[((Δ1)^2)+((Δ2)^2)]))/(2Δ1 Δ2))]);
A1 = {1, -((4((ω1)^2)Δ1-2)/(2((ω1)^2)Δ2))};
A2 = {1, -((4((ω2)^2)Δ1-2)/(2((ω2)^2)Δ2))};
X1 = (A1/(Sqrt[A1 . M . A1]));
X2 = (A2/(Sqrt[A2 . M . A2]));

(*Integración numérica de las ecuaciones de movimiento exactas*)

Module[{x1, x2, t}, 
sol = NDSolve[
{
(*Ecuaciones*)
4((Δ1)^2)ϕ1''[t] + Δ1 Δ2 (ϕ2''[t] (1 + Cos[ϕ1[t]- ϕ2[t]]) +  ((ϕ2'[t])^2) Sin[ϕ1[t] - ϕ2[t]]) == - 2 Δ1 Sin[ϕ1[t]], 
2((Δ2)^2)ϕ2''[t] + Δ1 Δ2(ϕ1''[t] (1 + Cos[ϕ1[t] - ϕ2[t]]) - ((ϕ1'[t])^2) Sin[ϕ1[t] - ϕ2[t]]) == - Δ2 Sin[ϕ2[t]],

(*Condiciones iniciales*)
ϕ1[0] == (ϕ10 = 0), 
ϕ2[0] == (ϕ20 = 0),
ϕ1'[0] ==(ϕ1p0 =5/10), 
ϕ2'[0] == (ϕ2p0 =  -1)
}, 

{ϕ1, ϕ2},
{t, 0, tmax = 400}, 

WorkingPrecision->20, MaxSteps->(10^6)
][[1]];
]

{ϕ1sol = sol[[1]][[2]], ϕ2sol = sol[[2]][[2]]};

(*La energía mecánica -- para verificar a grandes rasgos que la integración numérica respete su conservación*)

T[t_]:= 2((Δ1)^2)((ϕ1sol'[t])^2)+ ((Δ2)^2) ((ϕ2sol'[t])^2)+ Δ1 Δ2 ϕ1sol'[t] ϕ2sol'[t] ( 1 + Cos[ϕ1sol[t] - ϕ2sol[t]]) - 2Δ1 Cos[ϕ1sol[t]] - Δ2 Cos[ϕ2sol[t]];
meanT = (1/tmax)NIntegrate[T[t], {t, 0, tmax}];

(* La solución de pequeñas oscilaciones *)

ϕ0 = {ϕ10, ϕ20};
ϕp0 = {ϕ1p0, ϕ2p0};
ϕosc[t_]:= (X1 . M . ϕ0 Cos[ω1 t] X1+ (1/(ω1))X1. M . ϕp0 Sin[ω1 t] X1)+ (X2 . M . ϕ0 Cos[ω2 t] X2+ (1/(ω2))X2. M . ϕp0 Sin[ω2 t] X2);
ϕ1osc[t_]:= ϕosc[t][[1]];
ϕ2osc[t_]:= ϕosc[t][[2]];

(*Un gráfico de la solución "exacta" y verificación de que la energía se conserva en la integración numérica *)

{Plot[{ϕ1sol[t], ϕ2sol[t]}, {t, 0, 200}, ImageSize->600], Plot[((T[t]-meanT)/meanT), {t, 0, tmax}, ImageSize-> 600, PlotRange-> All]}

Esta celda muestra una animación del sistema para la solución “exacta”.

(*Se incluyen distintas marcas para seguir las posiciones y orientaciones de los cilindros*)

mr = 0.15;
DynamicModule[{t},
Manipulate[
ϕ1t = ϕ1sol[t];
ϕ2t = ϕ2sol[t];
α1t= -((Δ1)/r1)ϕ1t;
α2t = -(1/r2)(Δ1 ϕ1t + Δ2 ϕ2t);
r1t = (r0 - r1) {Sin[ϕ1t], -Cos[ϕ1t]};
r2t = r1t + (r1 - r2) {Sin[ϕ2t], -Cos[ϕ2t]};
marca11 = r1t + (punto1 = (r1-mr) {Sin[α1t], -Cos[α1t]});
marca12 = r1t - punto1;
marca13 = r1t - (punto2 =(r1-mr) {Cos[α1t], Sin[α1t]});
marca14 = r1t + punto2;
marca21 = r2t + (punto3 = (r2-mr) {Sin[α2t], -Cos[α2t]});
marca22 = r2t - punto3;
marca23 = r2t + (punto4 = (r2-mr) {Cos[α2t], Sin[α2t]});
marca24 = r2t - punto4;

base = Graphics[{Red, Thick, Circle[{0,0}, r0], Disk[{0,0}, 0.2]}];

Show[
base,
Graphics[{Magenta,Thick, Dashed, Line[{{0,0} , r1t}], Line[{r1t, r2t}]}], 

Graphics[
{Cyan, Thick, Circle[r1t, r1], White, Circle[r2t, r2], 
Cyan, Line[{marca11, marca12}], Line[{marca13, marca14}],
Pink, Disk[marca11, mr], Cyan, Disk[marca12, mr], Disk[marca13, mr], Disk[marca14, mr],
White, Line[{marca21, marca22}], Line[{marca23, marca24}], 
Green, Disk[marca21, mr], White, Disk[marca22, mr], Disk[marca23, mr], Disk[marca24, mr]
}
],

PlotRange-> {{-r0-2, r0+2}, {-r0-2, r0+2}}, 
ImageSize-> 600,
Background-> Black
], 
{t, 0, tmax, 0.01, Appearance-> "Open", AnimationRate-> 5}
]
]

De nuevo, el gráfico de la solución “exacta”.

Plot[{ϕ1sol[t], ϕ2sol[t]}, {t, 0, tmax}, ImageSize->600]

Comparación entre la solución “exacta” y la de pequeñas oscilaciones.

{Plot[{ϕ1sol[t], ϕ1osc[t]}, {t, 0, 50}, ImageSize->600], Plot[{ϕ2sol[t], ϕ2osc[t]}, {t, 0, 100}, ImageSize->600]}

Precisión de la integración

Integra numéricamente el sistema variando la precisión “WorkingPrecision”, para una condición inicial que da amplitudes de movimiento grandes.
Puede tomar algunos minutos

Module[{r0, r1, r2, Δ1, Δ2,ϕ10, ϕ20, ϕ1p0, ϕ2p0, tmax, ϕ1sol, sol, ϕ2sol},
r0 =8;
r1 =5;
r2 =1;
Δ1 = ((r0-r1)/r2);
Δ2 = ((r1-r2)/r2);

Monitor[
tplot = Table[
Module[{ϕ1, ϕ2, t}, 
sol = NDSolve[
{
4((Δ1)^2)ϕ1''[t] + Δ1 Δ2 (ϕ2''[t] (1 + Cos[ϕ1[t]- ϕ2[t]]) +  ((ϕ2'[t])^2) Sin[ϕ1[t] - ϕ2[t]]) == - 2 Δ1 Sin[ϕ1[t]], 
2((Δ2)^2)ϕ2''[t] + Δ1 Δ2(ϕ1''[t] (1 + Cos[ϕ1[t] - ϕ2[t]]) - ((ϕ1'[t])^2) Sin[ϕ1[t] - ϕ2[t]]) == - Δ2 Sin[ϕ2[t]],
ϕ1[0] == (ϕ10 = 0), 
ϕ2[0] == (ϕ20 = 0),
ϕ1'[0] ==(ϕ1p0 =12/10), 
ϕ2'[0] == (ϕ2p0 =  -1)
}, 
{ϕ1, ϕ2},
{t, 0, tmax = 400}
, WorkingPrecision->n, MaxSteps->(10^8)
][[1]];
];
{ϕ1sol = sol[[1]][[2]], ϕ2sol = sol[[2]][[2]]};
gtemp = Plot[{ϕ1sol[t], ϕ2sol[t]}, {t, 0, tmax}, ImageSize->600],
{n, 20, 50}
];, 
{n, Show[gtemp, ImageSize-> 600]}
]
]

Este gráfico muestra como cambia la solución a medida que se aumenta la precisión -- 
Notar que a medida que aumenta la precisión el comportamiento a largos tiempos se va estabilizando.

Manipulate[Show[tplot[[i]], PlotRange-> {-25, 25}], {i, 1, Length@tplot, 1, Appearance-> "Open"}]

Todos los gráficos en uno solo.

Show[tplot, PlotRange-> All]