Сделаем анимацию движения простого гироскопа. Цель в том, чтобы объяснить, как описывающие движение гироскопа уравнения Эйлера для вектора угловой скорости дают его компоненты в связанной с телом «модельной» системе отсчета (в ней моменты инерции не изменяются, и в ней легко задать внешний гравитационный вращающий момент). Задача состоит в том, чтобы соединить эти компоненты с лабораторной системой отсчета. Это делается с помощью двух углов поворота θ и . Интересно ввести трение, чтобы объяснить, например, спящий режим волчка:

• restart: with(linalg):

Начнем с вектора полного углового момента:

• Lv:=vector([-diff(theta(t),t)*I1,
  -diff(phi(t),t)*sin(theta(t))*I2,
  (diff(phi(t),t)*cos(theta(t))+w0)*I3]);

Вектор угловой скорости («большое» Ω, вращения нет): он определяет модельную систему отсчета, связанную с телом.

• Wv:=vector([-diff(theta(t),t),
  -diff(phi(t),t)*sin(theta(t)),
  diff(phi(t),t)*cos(theta(t))]);

Скорость изменения вектора углового момента задается векторным произведением Ω и L плюс скорость изменения «модельной» системы отсчета тела. Получим векторное произведение:

• WcrL:=map(simplify,crossprod(Wv,Lv));

Определим вектор с производной по времени от углового момента в «модельной» неинерциальной системе отсчета тела:

• Ldotv:=map(expand,map(diff,op(Lv),t));

Получили общую скорость изменения L (dL/dt_lab, которая равна крутящему моменту), выраженную в единичных векторах «модельной» системы:

• LdotS:=map(simplify,evalm(Ldotv+WcrL));

Вместо уравнений Эйлера применим теорему об угловом моменте: гравитационный момент всегда около оси x (ось, обозначаемая 1):

• eq1:=LdotS[1]=-M*g*l*sin(theta(t));

• eq2:=LdotS[2]=0;

• eq3:=LdotS[3]=0;

Теперь подставим значения для моментов инерции:

• l:=1: g:=0: M:=1: I1:=1/4+1: I2:=I1: I3:=1/2: w0:=20:

Начальное условие: вначале горизонтальное направление гироскопа по оси X в лабораторной системе, и затем он отпущен.

• th0:=Pi/2: ph0:=0: thdot0:=0: phidot0:=0:

Пара eq1/eq3 не работает! Это потому, что eq3 представляет собой только сохранение трехкомпонентного углового момента импульса:

• eq3;

• t:='t': sol:=dsolve({eq1,eq2,theta(0)=th0,
  D(theta)(0)=thdot0, phi(0)=ph0,D(phi)(0)=phidot0},
  {theta(t),phi(t)},numeric,method=gear,output=listprocedure);

• Theta:=subs(sol,theta(t)): Phi:=subs(sol,phi(t)):
• plot(['Theta(t)','Phi(t)'],t=0..20,color=[blue,green],
  numpoints=1000,thickness=3);

Зеленая кривая на рис. 4.15.1 () описывает прецессию, в то время как колебание в Ω (синяя кривая) соответствует нутации около некоторого среднего значения (которое определяется начальными условиями; см. также рис. 4.15.2).

Рис. 4.15.1

• plot('Theta(t)',t=0..4,numpoints=1000,thickness=3);

Рис. 4.15.2

Чтобы нарисовать движение гироскопа, надо перейти в лабораторную систему координат. Определяем декартовы компоненты вектора положения в единицах длины плеча l, полярного угла (угла нутации) Ω и азимутального угла (угла прецессии) :

• with(plottools):

Пишем анимацию, т. е. определяем последовательность во времени 3D-графиков и используем команду «стрелка» для отображения Ω -вектора в лабораторной системе (рис. 4.15.3). Чтобы построить изменение во времени вектора оси гироскопа, берутся численные решения для θ(t) и (t):

• N:=300: for i from 1 to N do: t:=0.03*i:
  Prec:=vector([l*sin(Theta(t))*cos(Phi(t)),
  l*sin(Theta(t))*sin(Phi(t)),l*cos(Theta(t))]) ;
  ll[i] := arrow(vector([0,0,0]), Prec, .2, .4, .1):
  od:
• i:='i': plots[display3d]([seq(ll[i],i=1..N)],
  color=red,orientation=[66,48],axes=boxed,insequence=true);

Рис. 4.15.3