Group Abstract Group Abstract

Message Boards Message Boards

Faster ways to solve symmetric top problem?

Ox Clouding
Posted 7 years ago

I'm currently trying to solve the symmetric top problem using Mathematica. It's essentially solving this differential equation:

Sec[\[Theta][t]/
   2] (1 - Cos[\[Theta][t]] + 2 Cos[2 \[Theta][t]] + 
    2 (1 + Cos[\[Theta][t]]) Derivative[1][\[Theta]][t]^2) == 0

(There are parameters, but I have set them to some numerical values in this equation).

I found this implementation Nutation of a Symmetric Top - Wolfram Demonstrations Project but it's quite complicated. So I try to do it with NDSolve:

eqN={Sec[\[Theta][t]/
    2] (1 - Cos[\[Theta][t]] + 2 Cos[2 \[Theta][t]] + 
     2 (1 + Cos[\[Theta][t]]) Derivative[1][\[Theta]][t]^2) == 
  0, \[Theta][0] == 1};

solN = NDSolve[eqN, \[Theta][t], {t, 0, 4}]

On my computer, it gives warnings like: "Maximum number of xxx steps reached at the point t == xxx". And I plot out the solutions using:

plotf = \[Theta][t] /. solN; Plot[plotf, {t, 0, 4}]

The result is like this:

Plot of theta

It looks like the solver won't proceed past when it reach the max/min point of theta, but before that everything is fine.

So here is my question: Without actually manually doing integration like the one in Demonstrations project did, are there any ways to solve it using NDSolve capability of Mathematica?

Thanks in advance!
POSTED BY: Ox Clouding
17 Replies
Sort By:
Ox Clouding
Ox Clouding
Posted 7 years ago
POSTED BY: Ox Clouding
POSTED BY: Neil Singer
POSTED BY: Neil Singer
Ox Clouding
Ox Clouding
Posted 7 years ago
POSTED BY: Ox Clouding
POSTED BY: Neil Singer
Ox Clouding
Ox Clouding
Posted 7 years ago
POSTED BY: Ox Clouding
Ox Clouding
Ox Clouding
Posted 7 years ago
POSTED BY: Ox Clouding
POSTED BY: Neil Singer
Ox Clouding
Ox Clouding
Posted 7 years ago
POSTED BY: Ox Clouding

That should teach me not to post without actually looking at the equations in detail! Sorry for the misinformation. Good catch Mariusz!

from what I see on the Systemmodeler webpage, it's don't have a physics interface built-in that one can draw, say, a ellipsoid and specify the gravity and fix point. So how does one simulate this? Do we have to manually put in the differential equation?

As to SystemModeler, you would model the top as a rigid body using the Multibody library elements -- specify the mass and inertia properties, the location of the center of mass and the contact point. Constrain the contact point on a surface and give it some initial conditions. SystemModeler derives the equations of motion and solves them, including an animation. There is a gyroscope example in the WSM examples that shows a rotating, swinging body. Unfortunately I do not see a closer example. If I get a chance, I will knock up a quick model and post it for you.

Regards,

Neil
POSTED BY: Neil Singer

Try:

eqN = {Sec[?[t]/2] (1 - Cos[?[t]] + 2 Cos[2 ?[t]] + 2 (1 + Cos[?[t]]) Derivative[1][?][t]^2) == 
      0, ?[0] == 1};
solN = NDSolve[eqN, ?[t], {t, 0, 4*Pi}, Method -> {"EquationSimplification" -> "Residual"}]
Plot[Evaluate[?[t] /. solN], {t, 0, 4*Pi}]

enter image description here

1:

 eq = {Sec[?[t]/2] (1 - Cos[?[t]] + 2 Cos[2 ?[t]] + 2 (1 + Cos[?[t]]) Derivative[1][?][t]^2) == 0};
 INIT = Derivative[1][?][0] /. Solve[eq /. t -> 0 /. ?[0] -> 1, Derivative[1][?][0]](*A second initial conditions *)

Then we have: $\theta '(0)=-\sqrt{\frac{-1+\cos (1)-2 \cos (2)}{2 (1+\cos (1))}}$ and $\theta '(0)=\sqrt{\frac{-1+\cos (1)-2 \cos (2)}{2 (1+\cos (1))}}$

 sol = ParametricNDSolve[{eq, ?[0] == 1, ?'[0] == a}, ?, {t, 0, 4*Pi}, {a}, Method -> {"EquationSimplification" -> "Residual"}];
 Plot[{?[(INIT[[1]])][t] /. sol, ?[(INIT[[2]])][t] /. sol}, {t, 0, 4*Pi}]

enter image description here

2: 

 eq = {Sec[?[t]/2] (1 - Cos[?[t]] + 2 Cos[2 ?[t]] + 2 (1 + Cos[?[t]]) Derivative[1][?][t]^2) == 0};
 INIT = Derivative[1][?][0] /. Solve[eq /. t -> 0 /. ?[0] -> 1, Derivative[1][?][0]](*A second initial conditions *)
  sol = ParametricNDSolve[{D[eq,t], ?[0] == 1, ?'[0] == a}, ?, {t, 0, 4*Pi}, {a}];
  Plot[{?[(INIT[[1]])][t] /. sol, ?[(INIT[[2]])][t] /. sol}, {t, 0, 4*Pi}]
POSTED BY: Mariusz Iwaniuk
Ox Clouding
Ox Clouding
Posted 7 years ago

It works!

But I am more confused. Why this solving method avoids the problem? I couldn't find any description of this "Residual" thing in Mathematica documentation. Is it more robust than default method or what?

Thank you!
POSTED BY: Ox Clouding

Hello

I'm ordinary user not Jedi Master in Mathematica. I only know how much it is written in the documentation. About "Residual" find here: (NDSolve->Options->Method->EquationSimplification)

Regards MI

POSTED BY: Mariusz Iwaniuk

Ox,

NDSolve is solving the equation, however, your solution is becoming complex (you are integrating into a region in which the derivative of Theta is the square root of a negative number.

This suggests that either your equations are wrong (for the physics you are modeling) or the top enters a region where those equations no longer apply (ie it fall over).

I would start by checking the physics.

One way to try this is to use Wolfram SystemModeler to model the top.

I hope this helps.

Regards,

Neil
POSTED BY: Neil Singer
Ox Clouding
Ox Clouding
Posted 7 years ago

Well, I'm certain the equations are physical since I obtained it directly from Goldstein's book, and the parameters are physical too. It's would be strange if it leads to complex numbers.

Regardless, Mariusz Iwaniuk posted a solution to the exactly same equation, so I don't think the equations are wrong. If so, then, why Mathematica would pick up complex numbers with default solver?

Thank you!
POSTED BY: Ox Clouding
Ox Clouding
Ox Clouding
Posted 7 years ago

This maybe off topic but, from what I see on the Systemmodeler webpage, it's don't have a physics interface built-in that one can draw, say, a ellipsoid and specify the gravity and fix point. So how does one simulate this? Do we have to manually put in the differential equation?
POSTED BY: Ox Clouding
Ox Clouding
Ox Clouding
Posted 7 years ago

[Double posted, please delete]
POSTED BY: Ox Clouding
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard