Описывается свободное вращение твердого тела. Для вектора угловой скорости, который определяет движение относительно инерциальной системы координат, решаются уравнения Эйлера, но при этом вектор представляется через компоненты в системе координат рассматриваемого тела (в ней моменты инерции сохраняются). В численном решении показывается, что движения вокруг главной оси неустойчивы, что связано с промежуточным значением момента инерции.

Для вектора угловой скорости используем [w1, w2, w3], а для трех главных моментов инерции – I1, I2, I3. N1, N2, N3 – это составляющие вращающего момента тела в собственной системе координат, которые для свободного вращения равны нулю:

• E1:=I1*diff(w1(t),t)-(I2-I3)*w2(t)*w3(t)=N1;

• E2:=I2*diff(w2(t),t)-(I3-I1)*w3(t)*w1(t)=N2;

• E3:=I3*diff(w3(t),t)-(I1-I2)*w1(t)*w2(t)=N3;

• E1:=subs(N1=0,E1): E2:=subs(N2=0,E2): E3:=subs(N3=0,E3):

В единицах СИ:

• M:=1.63;

• a:=235*10^(-3);

• b:=154*10^(-3);

• c:=17*10^(-3);

Чтобы найти три главных момента инерции, используем результат для прямоугольной пластины, вращающейся вокруг центра масс:

• I1:=M/12*(b^2+c^2);

• I2:=M/12*(a^2+c^2);

• I3:=M/12*(a^2+b^2);

Уравнения готовы для численного решения, и мы можем записать начальные условия. Движение начинается вокруг оси 1 (наименьший момент инерции должен быть устойчивым), и добавляем небольшие вклады по остальным осям (при бросании тела, например кирпича, никогда не удастся идеально совпасть с осью симметрии).

Для главной оси вращения зададим начальное значение ≈ 1 оборот в секунду, а остальные ≈ 1 % от этой величины.

• E1;

• IC:=w1(0)=6.28,w2(0)=5/100,w3(0)=5/100;

• sol:=dsolve({E1,E2,E3,IC},{w1(t),w2(t),w3(t)},
  numeric, output=listprocedure):
• W1:=subs(sol,w1(t)): W2:=subs(sol,w2(t)):
  W3:=subs(sol,w3(t)):
• plot(['W1(t)','W2(t)','W3(t)'],t=0..10,
  color=[red,blue,green],axes=boxed,thickness=3);

Результат показан на рис. 4.13.1.

Рис. 4.13.1

Теперь сделаем то же самое для второй оси (в нашем случае – той, которая имеет промежуточное значение момента инерции: I1 < I2 < I3):

• IC:=w1(0)=5/100,w2(0)=6.28,w3(0)=5/100;

• sol:=dsolve({E1,E2,E3,IC},{w1(t),w2(t),w3(t)},
  numeric, output=listprocedure):
• W1:=subs(sol,w1(t)): W2:=subs(sol,w2(t)):
  W3:=subs(sol,w3(t)):
• plot(['W1(t)','W2(t)','W3(t)'],t=0..10,
  color=[red,blue,green],axes=boxed,thickness=3);

По синей кривой на рис. 4.13.2 видно, что, если продолжить изменение ориентации вращения в диапазоне между крайними значениями, то движение «перевернется».

Рис. 4.13.2

Кстати, этот вопрос интересен для теории устойчивости дифференциальных уравнений: как следует анализировать причину происходящего с рассматриваемой системой с помощью выкладок Maple.