Problema 15 de la guía de cuerpo rígido. Es un giróscopo 'rápido' pesado con un punto fijo.

La solución numérica:

Module[{ϕ, θ, ψ, t, I1 = 1, I3 = 2, ω3, ψp0, ϕp0, θp0, θ0, ϕ0, ψ0, g =10.5293},

(*Condiciones iniciales*)

(*Las velocidades*)
ψp0 = 5;
ϕp0 = 0;
θp0 = 0;

(*Las posiciones*)
ϕ0 = 0;
ψ0 = 0;
θ0 = π/2;     (*Varién esta condición inicial. En la guía esta puesta en π/2 únicamente para hacer las cuentas más simples*)

ω3 = ψp0 + ϕp0 Cos[θ0];

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

sol = NDSolve[
{
(*Las ecuaciones diferenciales*)

I1 θ''[t] == -I3 ω3 ϕ'[t] Sin[θ[t]] + I1 ((ϕ'[t])^2) Sin[θ[t]] Cos[θ[t]] + g Sin[θ[t]],
I1 ϕ''[t] ((Sin[θ[t]])^2)+2 I1 ϕ'[t] θ'[t] Sin[θ[t]] Cos[θ[t]] - I3 ω3 θ'[t] Sin[θ[t]] == 0,
ψ'[t] + ϕ'[t] Cos[θ[t]] == ω3, 

(*Condiciones iniciales. La ecuación para ψ es de primer orden, así que no es necesario dar ψ' inicial. *)

θ[0] == θ0, 
ϕ[0] == ϕ0,
ψ[0] == ψ0,
θ'[0] == θp0,
ϕ'[0] == ϕp0,

(*La integración se detiene cuando ϕ da una vuelta completa*)

WhenEvent[ϕ[t] > 2π, "StopIntegration"]
},
{ϕ, θ, ψ},
{t, 0, tmax = ∞}
][[1]]
];

(*La solución numéricca se usa ahora para definir las funciones ϕ, θ, y ψ. Por otro lado se extrae el período Tϕ.*)

{ϕ = sol[[1]][[2]], θ = sol[[2]][[2]], ψ = sol[[3]][[2]]}

Tϕ = sol[[1]][[2]][[1]][[1]][[2]]

(* ¿A dónde fue a parar θ luego de que ϕ variase en 2π? *)
θ[Tϕ]

Auxiliar: número de períodos en θ por cada vuelta de ϕ según la solución del problema linealizado: ((((I3 ω3))^2)/g)

((((2 5))^2)/10.5293)1.

Animación

(* Curva que sigue el extremo de la barra*)
ez = {0,0,1};
curva = ParametricPlot3D[(largo=2)(Cos[θ[t]] ez - Sin[θ[t]] {-Sin[ϕ[t]], Cos[ϕ[t]], 0}), {t, 0, Tϕ}, PlotStyle->Orange, PlotPoints->100, Background->Black]

(*Grilla de fondo*)
lmax = 4;
bajo = 4;
grilla =Show[Graphics3D[{Opacity[0.05],White, Table[Line[{{-lmax, i, -bajo}, {lmax, i, -bajo}}], {i, -lmax, lmax, 0.5}], Table[Line[{{i, -lmax, -bajo}, {i, lmax, -bajo}}], {i, -lmax, lmax, 0.5}]}, Boxed-> False], Background-> Black]

La imagen que sirve de textura para el disco

im = ;

Un gráfico de prueba para t fijo

tpr = 1;

x = ϕ[tpr]; y = θ[tpr]; z = ψ[tpr];
ex = {1,0,0};
ey = {0,1,0};

eρ = {Cos[x], Sin[x], 0}; 
eϕ = {-Sin[x], Cos[x], 0};
eϕm = {-Sin[ϕm = 2(π/(Tϕ))tpr], Cos[ϕm],0};
ez = {0,0,1}; 
eψ = Cos[y] eϕ + Sin[y] ez;

e1 = Cos[z] eρ + Sin[z] eψ;
e2 = -Sin[z] eρ + Cos[z] eψ;
e3 = Cos[y] ez - Sin[y] eϕ;

Show[

curva,
grilla,
Graphics3D[{Red,Specularity[White, 100], Tube[{{0,0,0},largo e3}, (rt = 0.05)0.5],White, Tube[{largo e3, (largo+rt/2) e3},1.5rt]}],

(*El disco se grafica como si fuera medio elipsoide muy aplanado*)

ParametricPlot3D[(apr=1/2.)( Cos[ϕ] Sin[θ] e1+  Sin[ϕ] Sin[θ] e2) +0.1 Cos[θ] e3, {ϕ, 0, 2π}, {θ, 0, π/2},
AspectRatio->Automatic,
PlotPoints->20,
Mesh-> False,
Boxed-> False,
Axes-> False,
PlotStyle-> {Directive[Texture[im]]},
TextureCoordinateFunction->({((Cos[2π#4] Sin[#5 π/2]+1)/2), ((Sin[2π #4] Sin[#5 π/2]+1)/2)}&)
],

PlotRange-> 10 {{-1, 1}, {-1, 1}, {-1, 1}},
Background-> Black,
Boxed-> False,
Ticks-> False,
Axes-> False,
AxesOrigin-> {0,0,0},
ViewVector-> {-6eϕm+0ex+1 ez, {0,0,0}}, 
ViewAngle-> 45°,
ViewVertical-> {0,0,1}
]

La animación en sí misma

DynamicModule[{tpr},
Manipulate[

x = ϕ[tpr]; y = θ[tpr]; z = ψ[tpr];

ex = {1,0,0};
ey = {0,1,0};
ez = {0,0,1}; 

eρ = {Cos[x], Sin[x], 0}; 
eϕ = {-Sin[x], Cos[x], 0};
eψ = Cos[y] eϕ + Sin[y] ez;

(*Este vector avanza con la veloccidad media de precesión*)
eϕm = {-Sin[ϕm = 2(π/(Tϕ))tpr], Cos[ϕm],0}; 

e1 = Cos[z] eρ + Sin[z] eψ;
e2 = -Sin[z] eρ + Cos[z] eψ;
e3 = Cos[y] ez - Sin[y] eϕ;

Show[

curva,
grilla,
Graphics3D[{Red,Specularity[White, 100], Tube[{{0,0,0},largo e3}, (rt = 0.05)0.5],White, Tube[{largo e3, (largo+rt/2) e3},1.5rt]}],

ParametricPlot3D[(apr=1/2.)( Cos[ϕ] Sin[θ] e1+  Sin[ϕ] Sin[θ] e2) +0.1 Cos[θ] e3, {ϕ, 0, 2π}, {θ, 0, π/2},
AspectRatio->Automatic,
 
(* Bajar el número de puntos si la animación se hace muy lenta. De todas formas la animación es una etapa previa a la salida de los archivos con las imágenes; es más una referencia que otra cosa*)
PlotPoints->10,

Mesh-> False,
Boxed-> False,
Axes-> False,
PlotStyle-> {Directive[Texture[im]]},
TextureCoordinateFunction->({((Cos[2π#4] Sin[#5 π/2]+1)/2), ((Sin[2π #4] Sin[#5 π/2]+1)/2)}&)
],

PlotRange-> 10 {{-1, 1}, {-1, 1}, {-1, 1}},
Background-> Black,
Boxed-> False,
Ticks-> False,
Axes-> False,
AxesOrigin-> {0,0,0},

(*La posición de la cámara. Aquí se mueve con el movimiento medio de ϕ *)
ViewVector-> {-6eϕm+0ex+1 ez, {0,0,0}},

ViewAngle-> 45°,
ViewVertical-> {0,0,1}
],

{tpr, 0, Tϕ, 0.01, Appearance->"Open", AnimationRate-> .5},
TrackedSymbols:> tpr
]
]

Salida de archivos

$HistoryLength = 0;

(*Esto rasteriza la imágenes a una buena resolución *)
antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]];

frame = 0;
nframes = 400;
Monitor[
tsal = Table[
++frame;
x = ϕ[tpr]; y = θ[tpr]; z = ψ[tpr];
ex = {1,0,0};
ey = {0,1,0};

eρ = {Cos[x], Sin[x], 0}; 
eϕ = {-Sin[x], Cos[x], 0};
eϕm = {-Sin[ϕm = 2(π/(Tϕ))tpr], Cos[ϕm],0};
ez = {0,0,1}; 
eψ = Cos[y] eϕ + Sin[y] ez;

e1 = Cos[z] eρ + Sin[z] eψ;
e2 = -Sin[z] eρ + Cos[z] eψ;
e3 = Cos[y] ez - Sin[y] eϕ;

imtemp = Show[
curva,
grilla,
Graphics3D[{Red,Specularity[White, 100], Tube[{{0,0,0},largo e3}, (rt = 0.05)0.5],White, Tube[{largo e3, (largo+rt/4) e3},rt]}],

ParametricPlot3D[(apr=1/2.)( Cos[ϕ] Sin[θ] e1+  Sin[ϕ] Sin[θ] e2) +0.1 Cos[θ] e3, {ϕ, 0, 2π}, {θ, 0, π/2},
AspectRatio->Automatic,
PlotPoints->50,
Mesh-> False,
Boxed-> False,
Axes-> False,
PlotStyle-> {Directive[Texture[im]]},
TextureCoordinateFunction->({((Cos[2π#4] Sin[#5 π/2]+1)/2), ((Sin[2π #4] Sin[#5 π/2]+1)/2)}&)
],

PlotRange->  {10{-1, 1},10 {-1, 1}, {-10, 1}},
Background-> Black,
Boxed-> False,
Ticks-> False,
ImageSize-> 600,
Axes-> False,

(*En este caso la cámara está fija*)
ViewVector-> {-4ey+2ex+2 ez, {0,0,0}}, 

ViewAngle-> 60°,
ViewVertical-> {0,0,1}
];

Export[StringJoin["C:\\users\\juan zanella\\desktop\\sal 1\\frame_", If[frame < 10, "000", If[frame < 100, "00", "0"]], ToString[frame], ".png"], imsal = antialias[imtemp, 2],"PNG"];,

{tpr, 0, Tϕ-Tϕ/nframes, Tϕ/nframes}
];, {frame, Show[imsal, ImageSize-> 400]}]

Pruebas auxiliares para el ploteo del disco

im = ;
apr = 1;
bpr = 1/10;

ParametricPlot3D[{apr Cos[ϕ] Sin[θ], apr Sin[ϕ] Sin[θ], bpr Cos[θ]}, {ϕ, 0, 2π}, {θ, 0, π/2},
AspectRatio->Automatic,
Mesh-> False,
Boxed-> False,
Axes-> False,
PlotStyle-> Texture[im],
TextureCoordinateFunction->({((Cos[2π#4] Sin[#5 π/2]+1)/2), ((Sin[2π #4] Sin[#5 π/2]+1)/2)}&)
]

tpr = 1;

x = ϕ[tpr]; y = θ[tpr]; z = ψ[tpr];
eρ = {Cos[x], Sin[x], 0}; eϕ = {-Sin[x], Cos[x], 0}; ez = {0,0,1}; eψ = Cos[y] eϕ + Sin[y] ez;

e1 = Cos[z] eρ + Sin[z] eψ;
e2 = -Sin[z] eρ + Cos[z] eψ;
e3 = Cos[y] ez - Sin[y] eϕ;


Show[

ParametricPlot3D[apr Cos[ϕ] Sin[θ] e1+ apr Sin[ϕ] Sin[θ] e2 +bpr Cos[θ] e3, {ϕ, 0, 2π}, {θ, 0, π/2},
AspectRatio->Automatic,
Mesh-> False,
Boxed-> False,
Axes-> False,
PlotStyle-> {Directive[Texture[im]]},
TextureCoordinateFunction->({((Cos[2π#4] Sin[#5 π/2]+1)/2), ((Sin[2π #4] Sin[#5 π/2]+1)/2)}&)
],

PlotRange-> 2 {{-1, 1}, {-1, 1}, {-1, 1}},
Background-> Black,
ImageSize-> 600
]

Búsqueda auxiliar de g para órbita cuasi-cerrada:

(*Se construye una tabla con soluciones para todo un intervalo de valores de g y se busca el valor que hace que luego de un período de movimiento en ϕ el valor de θ haya regresado lo más cerca posible de su valor inicial de π/2*)

tθ = Table[
Module[{ϕ, θ, ψ, t, I1 = 1, I3 = 2, ω3, ψp0, ϕp0, θp0, θ0, ϕ0, ψ0, g},

g = gpr;
ψp0 = 5;
ϕp0 = 0;
θp0 = 0;
θ0 = π/2;
ϕ0 = 0;
ψ0 = π/3;

ω3 = ψp0 + ϕp0 Cos[θ0];

sol = NDSolve[
{
I1 θ''[t] == -I3 ω3 ϕ'[t] Sin[θ[t]] + I1 ((ϕ'[t])^2) Sin[θ[t]] Cos[θ[t]] + g Sin[θ[t]],
I1 ϕ''[t] ((Sin[θ[t]])^2)+2 I1 ϕ'[t] θ'[t] Sin[θ[t]] Cos[θ[t]] - I3 ω3 θ'[t] Sin[θ[t]] == 0,
ψ'[t] + ϕ'[t] Cos[θ[t]] == ω3, 

θ[0] == θ0, 
ϕ[0] == ϕ0,
ψ[0] == ψ0,
θ'[0] == θp0,
ϕ'[0] == ϕp0,

WhenEvent[ϕ[t] > 2π, "StopIntegration"]
},
{ϕ, θ, ψ},
{t, 0, tmax = 100}
][[1]]
];
θpr = sol[[2]][[2]];
Tϕ = sol[[1]][[2]][[1]][[1]][[2]];
{gpr,  θpr[Tϕ]-π/2}, {gpr,10.525,10.535, 0.0001}];

ListPlot[tθ,  PlotRange-> All]

Position[tθ, Min[tθ[[All, 2]]]]

tθ[[44]]