Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

How can I simulate the rotation of a pot lid about an axis that is quickly rotating as in the following image? I'm particular interested in the motion of the red dot.

rotating pot lid

Here is my attempt:

φ = π/6; r = 1;
line = Rotate[{Line[
     Table[{Cos[θ], Sin[θ], 0}, {θ, 0, 
       2 π, π/20}]], {Red, Pink, PointSize[Large], 
     Point[{1, 0, 0}]}}, φ, {0, 1, 0}];
    Rotate[line, θ, {0, 0, 
      1}], θ Cos[φ], {Sin[φ] Cos[θ], 
     Sin[φ] Sin[θ], Cos[φ]}, {0, 0, 0}], {Line[
     Table[{r Cos[φ] Cos[θ1], 
       r Cos[φ] Sin[θ1], -r Sin[φ]}, {θ1, 
       0, θ, π/20}]]}, {Blue, PointSize[Large], 
    Point[{0, 0, 0}]}}, Axes -> True, PlotRange -> 2, 
  AxesLabel -> {"x", "y", "z"}], {θ, 0, 2 π}]

visualization attempt

share|improve this question
This might be what you're looking for. – cormullion Apr 14 '13 at 22:16
Everything happens so fast! To the down-voters, I'm surprised that the "effort" is so much important than the question itself on this sites. Actually I did some work about the problem before I post it here. Anyway, I added some code to the question if that can make you happier. – xslittlegrass Apr 14 '13 at 23:58
Following your edit, I've undeleted the question, removed my downvote, edited your post and reopened it. Please note that emphasizing on some "basic effort" is not something strange. To you or someone who only asks questions, it might seem like we're being annoying pricks, but if you actually hang out on the site answering questions, you'll soon find out that there are quite a number of users who just want us to do all their coding for them. Those are all shot down and closed/deleted quickly for good reason — if we let it become the norm, then soon people will lose interest in sticking around. – R. M. Apr 15 '13 at 5:50
@rm-rf Thanks for doing that. I now understand the situation here. I will pay more attention to the form of the questions in the future to make your work easier. I really appreciate what you do to this community.Initially I didn't put any code in the question was because I'm more interested in learning the ideas and concepts from you creative experts tackling a real physical problem, than mealy code techniques, and I was a little worried that some codes or words pointing to specific directions would harms the creative diverse ideas on the problem,maybe that worry is not necessary.Thanks again. – xslittlegrass Apr 15 '13 at 18:36
@MichaelE2 I'm interested in making it better in the sense of visualization and experimentable. For example, if I could experiment with the initial angule, initial angular momentum, etc., that would be great. Thanks. – xslittlegrass Apr 16 '13 at 15:23
up vote 12 down vote accepted

In my opinion, this problem is the same as a rolling ellipsoid on a plane without slipping, which can solved by Euler equations for rigid body.

The Euler equations in body axes are

$$ I_1\dot\omega_1-\omega_2\omega_3(I_2-I_3)=0\\ I_2\dot\omega_2-\omega_3\omega_1(I_3-I_1)=0\\ I_3\dot\omega_3-\omega_1\omega_2(I_1-I_2)=0 $$

with initial conditions $$ \omega_1(0)=0,\omega_2(0)=0,\omega_3(0)=\omega_{30} $$

where $\omega_i$ and $I_i$ are the angular frequency and momentum of inertia in the three principle axes.

Then the three Euler angles $\theta, \psi, \phi$ are solved using the Euler angle equations:

\begin{eqnarray} \dot\phi\sin\theta\sin\psi+\dot\theta\cos\psi&=&\omega_1\\ \dot\phi\sin\theta\cos\psi-\dot\theta\sin\psi&=&\omega_2\\ \dot\phi\cos\theta+\dot\psi&=&\omega_3 \end{eqnarray}

After get the Euler angles function of time, we can easily make a movie. The line on the ellipsoid is called polhode (black curve) and the line on the plane is called herpolhode (red curve).

a = 3; b = 4; c = 1;(*The Length of the cubic*)
I1 = M/12 (b^2 + c^2); I2 = M/12 (a^2 + c^2); I3 = M/12 (a^2 + b^2); M = 12;(*Three momentum of inertia in principle axes*)
w10= 1; w20 = 2; w30 = 1;(*The initial conditions*)
L = Sqrt[w10^2 I1^2 + w20^2 I2^2 + w30^2 I3^2];(*The magnitude of the angular momentum*)  
T = 1/2 w10^2 I1 + 1/2 w20^2 I2 + 1/2 w30^2 I3 ;(*The kenetic engergy*)
d = Sqrt[2 T]/L;(*The distance of the origin to the plane*)

solw = NDSolve[{I1 w1'[t] == (I2 - I3) w2[t] w3[t], 
     I2 w2'[t] == (I3 - I1) w1[t] w3[t], 
     I3 w3'[t] == (I1 - I2) w2[t] w1[t], w1[0] == w10, w2[0] == w20, 
     w3[0] == w30}, {w1, w2, w3}, {t, 0, 20}] // 
   Flatten;(*Solve Euler equations, in the body axes*)

solang = NDSolve[{\[Phi]'[
          t] Sin[\[Theta][t]] Sin[\[Psi][t]] + \[Theta]'[
          t] Cos[\[Psi][t]] == w1[t],
      \[Phi]'[t] Sin[\[Theta][t]] Cos[\[Psi][t]] - \[Theta]'[
          t] Sin[\[Psi][t]] == w2[t],
      \[Phi]'[t] Cos[\[Theta][t]] + \[Psi]'[t] == w3[t], \[Phi][0] == 
       0.01, \[Theta][0] == 0.01, \[Psi][0] == 0.01} /. 
     solw, {\[Phi], \[Theta], \[Psi]}, {t, 0, 20}] // 
   Flatten;(*Solve the equations of Euler angles*)

w1[t_] = w1[t] /. solw; w2[t_] = w2[t] /. solw; w3[t_] = w3[t] /. solw;
\[Phi][t_] = \[Phi][t] /. solang; \[Theta][t_] = \[Theta][t] /. 
  solang; \[Psi][t_] = \[Psi][t] /. solang;
A[t_] = {{Cos[\[Psi][t]], Sin[\[Psi][t]], 0}, {-Sin[\[Psi][t]], 
     Cos[\[Psi][t]], 0}, {0, 0, 1}}.{{1, 0, 0}, {0, Cos[\[Theta][t]], 
     Sin[\[Theta][t]]}, {0, -Sin[\[Theta][t]], 
     Cos[\[Theta][t]]}}.{{Cos[\[Phi][t]], Sin[\[Phi][t]], 
     0}, {-Sin[\[Phi][t]], Cos[\[Phi][t]], 0}, {0, 0, 
     1}};(*Transformation matrix as a function of time*)

cpilherpolhode = 
  Compile[{{t, _Real}}, 
   Transpose[{{Cos[\[Psi][t]], Sin[\[Psi][t]], 0}, {-Sin[\[Psi][t]], 
        Cos[\[Psi][t]], 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 
        Cos[\[Theta][t]], Sin[\[Theta][t]]}, {0, -Sin[\[Theta][t]], 
        Cos[\[Theta][t]]}}.{{Cos[\[Phi][t]], Sin[\[Phi][t]], 
        0}, {-Sin[\[Phi][t]], Cos[\[Phi][t]], 0}, {0, 0, 1}}].(1/Sqrt[
      2 T] {w1[t], w2[t], w3[t]})];(*Compile the herpolhode equations*)

polhode = 
  Table[Point[(1/Sqrt[2 T] {w1[t], w2[t], w3[t]})], {t, 0, 10, 
     0.1}] /. solw;(*polhode*) 
elps = Scale[{Sphere[], 
      Line@Table[{1.01 Cos[u], 1.01 Sin[u], 0}, {u, 0., 2 Pi, 
         2 Pi/30}], i, {0, 1, 0}], {i, 0, \[Pi], \[Pi]/3}], 
      Table[{1.01 Cos[u], 1.01 Sin[u], 0}, {u, 0., 2 Pi, 
        2 Pi/30}], \[Pi]/2, {1, 0, 0}]}, {1/Sqrt[I1], 1/Sqrt[I2] , 1/
    Sqrt[I3] }, {0, 0, 0}];(*ellipsoid of inertia with lines on it*)

elpspolhode = {elps, polhode};(*combine the ellipsoid and the polhode*)

plane = Plot3D[
   L/(w30 I3) (d - (w10 I1)/L x - (w20 I2)/L y), {x, -0.4, 
    0.4}, {y, -0.5, 0.5}, Mesh -> None, 
   PlotStyle -> Opacity[0.5]];(*plane*)

      Rotate[elpspolhode, \[Phi][t], {0, 0, 1}], \[Theta][
       t], {Cos[\[Phi][t]], Sin[\[Phi][t]], 0}], \[Psi][
      t], {Sin[\[Theta][t]] Sin[\[Phi][
         t]], -Sin[\[Theta][t]] Cos[\[Phi][t]], Cos[\[Theta][t]]}]], 
   plane}, ParametricPlot3D[cpilherpolhode[tt], {tt, 0, t}, 
   PlotStyle -> Red], AspectRatio -> 1, PlotRange -> 0.5, 
  Boxed -> False, ViewAngle -> 0.24252597899424388`, 
  ViewPoint -> {0.09410273512100081`, -1.389522619330309`, 
  ViewVertical -> {0.4688415188450498`, 0.3712202655129461`, 
    0.8014880814332579`}], {t, 0.001, 20, 0.1}]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.