Рассчитаем матрицу момента инерции в некоторой системе координат и покажем, что можно найти такую систему, в которой матрица диагональна:

• restart; with(LinearAlgebra):
• with(plots): with(plottools):

Начнем с примера, в котором интегралы рассчитываются напрямую.

Имеем однородную среду с постоянной плотностью массы, задаваемой как отношение полной массы к полному объему.

• rho:=M/V;

Пусть тело имеет форму сигары – эллипсоида вращения, описываемого квадратичной формой:

• A:=Matrix([[1,1/2,1/3],[1/2,4,1],[1/3,1,4]]);

• A-Transpose(A);

Полученный нулевой результат показывает, что А выбрана действительной и симметричной:

• EVals:=evalf(Eigenvalues(A));

Числовой шум в собственных значениях (ненулевая комплексная часть) происходит при применении для расчета evalf с 10 знаками.

Матрица имеет положительные собственные значения, которые относятся к неглавным осям эллипсоида. Если бы один (или два) из них были отрицательны, то квадратичная поверхность представляла бы собой листы гиперболоида.

Результат можно предвидеть. В главных осях квадратичная форма имеет вид:

QF = EVals[1]*x^2 + EVals[2]*y^2 + EVals[3]*z^2.

Для выбора шкалы можно положить QF = 1 и получить, что собственные значения равны

EV[i] = 1/a_i2,

где a_i – неглавные оси.

Следовательно, объем эллипсоида рассчитывается как:

• V0:=evalf(4/3*Pi*mul(sqrt(1/EVals[i]),i=1..3));

• seq(Re(sqrt(1/EVals[i])),i=1..3);

Теперь вычислим квадратичную форму (рис. 4.7.1):

• Rv:=Vector([x,y,z]);

• QF:=Transpose(Rv) . A . Rv;

• QF:=simplify(QF);

• surf:=solve(QF-1,z);

• plot3d(surf[1],x=-1.1..1.1,y=-1.1..1.1,axes=boxed,
  grid=[30,30],style=patchcontour,shading=zhue,
  scaling=constrained);

Рис. 4.7.1

Получившийся график выглядит не так, как хотелось бы, поэтому применим implictplot (рис. 4.7.2).

• implicitplot3d(QF-1,x=-1.1..1.1,y=-1.1..1.1,z=-1.1..1.1,
  axes=boxed,numpoints=2000,scaling=constrained,
  style=patchcontour,shading=zhue);

Рис. 4.7.2

• EVecs:=evalf( Eigenvectors(A) );

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

• E1v:=convert(SubMatrix(EVecs[2],1..3,1..1),Vector);

• E1v:=Normalize(map(Re,E1v),Euclidean);

• E2v:=Normalize(map(Re,convert(SubMatrix(EVecs[2],1..3,
  2..2),Vector)),Euclidean);

• E3v:=Normalize(map(Re,convert(SubMatrix(EVecs[2],1..3,
  3..3),Vector)),Euclidean);

• E1 := plottools[arrow]([0, 0, 0], E1v, .2, .4, .1, color=red):
• E2 := plottools[arrow]([0, 0, 0], E2v, .2, .4, .1, color=blue):
• E3 := plottools[arrow]([0, 0, 0], E3v, .2, .4, .1, color=green):
• QFp:=implicitplot3d(QF-1,x=-1.1..1.1,y=-1.1..1.1,
  z=-1.1..1.1,axes=boxed,numpoints=2000,scaling=constrained,
  style=wireframe,color=gold):
• display(E1,E2,E3,QFp,axes=boxed);

Рис. 4.7.3

• E1v . E2v;

• E1v . E1v;

• E1v . E3v;

• E2v . E3v;

Для лучшего понимания геометрии повращаем картинку.

Изменим матрицу в квадратичной форме и повторим расчеты.

Продолжим конструирование ортогональной матрицы вне нормированных собственных векторов. Преобразуем матрицу, чтобы показать, что она диагонализована этим преобразованием:

• R:=Matrix([E1v,E2v,E3v]);
• Ri:=MatrixInverse(R);

Результаты применения этих команд (выдачи Maple) не приводятся. Их можно посмотреть самостоятельно.

Сделанная из собственных векторов (столбцы) матрица преобразования является ортогональной, т. е. обратная ей получается в результате транспонирования. Теперь рассчитаем преобразование матрицы А матрицей R. Правильный вид получается так (ниже покажем, что другой способ не работает!):

• Ri . A . R;

На диагонали преобразованной матрицы – собственные значения, а недиагональные элементы – нули (приближенно).

• R . A . Ri;

Результаты применения (выдачи Maple) этих команд не приводятся, чтобы можно было в этом убедиться самостоятельно.

Рассчитайте квадратичную форму с помощью диагонализованной матрицы. Сделайте недиагональные элементы нулями (точно), т. е. создайте диагональную форму из собственных значений, и рассмотрите уравнения для квадратичной формы (нет пересекающихся членов).

Выполните преобразование главных осей (как сделано выше) для некоторой другой матрицы А с тремя положительными собственными значениями.

Расчет объема. Можно ли рассчитать объем V в декартовых координатах? Выбираем квадратичную форму и масштабный множитель, равный 1.

Выполним интегрирование по очереди. Вначале интегрирование по z:

• I1:=int(1,z=fz[1]..fz[2]);

• fy:=solve(I1,y);

• plot([fy[1],fy[2]],x=-1.1..1.1,color=[red,blue],
  numpoints=800,scaling=constrained);

Результат представлен на рис. 4.7.4.

Рис. 4.7.4

Следуя интегрированию по z, нужно интегрировать I1 по площади, ограниченной кривыми fy[1] и fy[2]. Будем интегрировать по y между нижней (синей) и верхней (красной) кривыми. Затем будем интегрировать по x в пределах. Каковы они? Это точки, в которых производная y(x) становится бесконечной:

• x0:=solve(1/diff(fy[1],x),x);

Прямая попытка взять определенный интеграл неудачна:

• I2:=int(I1,y=fy[1]..fy[2]);
• I3idf:=int(expand(I2),x);

(Получившиеся в результате действия команд громоздкие выдачи не приводятся.)

В конце концов, можно бы взять этот интеграл численно. Для этого нужно знать диапазон x:

• evalf(x0);

• plot(I2,x=x0[1]..x0[2]);

Надо численно интегрировать получившуюся функцию (рис. 4.7.5) и исправить знак.

Рис. 4.7.5

• V:=evalf(Int(-I2,x=x0[1]..x0[2]));

• V0:=evalf(4/3*Pi*mul(sqrt(1/Eigenvalues(A)[i]),i=1..3),15);

Объемное интегрирование работает! Из трех интегралов два берутся аналитически, третий – численно.

Теперь рассчитаем матрицу момента инерции:

• M:=1; rho;

• MI:=Matrix(3,3):
• for i from 1 to 3 do: for j from 1 to i do:
• if i=1 and j=1 then w:=Rv[2]^2+Rv[3]^2; elif i=2 and j=2 then w:=Rv[1]^2+Rv[3]^2; elif i=3 and j=3 then w:=Rv[1]^2+Rv[2]^2; else w:=-Rv[i]*Rv[j]; fi:
• I1:=int(w*rho,z=fz[1]..fz[2]);
• I2idf:=unapply(int(I1,y),y);
• I2:=expand(I2idf(fy[2])-I2idf(fy[1]));
• I2:=simplify(I2);
• MI[i,j]:=evalf(Int(I2,x=x0[1]..x0[2]));
• od: od:
• print(MI);

Теперь заполним матрицу целиком с помощью симметрии:

• for i from 1 to 3 do: for j from i to 3 do: MI[i,j]:=MI[j,i]: od: od:
• print(MI);

• Eigenvalues(MI);

Это главные моменты инерции.

• EVecsMI:=Eigenvectors(MI);

(Получившаяся в результате действия команды громоздкая выдача не приводится.)

Нам нужны графики трех собственных векторов. Далее приведем только команды Maple. Результат их выполнения опустим. Его легко увидеть самому.

• E1vMI:=convert(SubMatrix(EVecsMI[2],1..3,1..1),Vector);
• E1vMI:=Normalize(map(Re,E1vMI),Euclidean);
• E2vMI:=Normalize(map(Re,convert(SubMatrix(EVecsMI[2],
  1..3,2..2),Vector)),Euclidean);
• E3vMI:=Normalize(map(Re,convert(SubMatrix(EVecsMI[2],
  1..3,3..3),Vector)),Euclidean);
• RMI:=Matrix([E1vMI,E2vMI,E3vMI]);

• R;
• E1MI := plottools[arrow]([0, 0, 0], E1vMI, .2, .4, .1, color=red):
• E2MI := plottools[arrow]([0, 0, 0], E2vMI, .2, .4, .1, color=blue):
• E3MI := plottools[arrow]([0, 0, 0], E3vMI, .2, .4, .1, color=green):

Собственные векторы матрицы MI дают оси, вращение тела вокруг которых возможно без смещений.

Заметьте, что если движение относительно осей с экстремальными собственными значениями стабильно, то для оси, обусловленной промежуточными собственными значениями, вращение нестабильно (есть колебания; нестабильность в уравнении Эйлера).

Поразительно, что в механике вращения для тела любой формы (например, помидора) можно найти три оси, вокруг которых возможно чистое вращение без смещающих импульсов. Наиболее близкая форма для описания томата – эллипсоид. Рисунок 4.7.6 показывает, что оси стабильного вращения совпадают с осями, перпендикулярными сечению эллипсоида.

• display(E1MI,E2MI,E3MI,QFp,axes=boxed);

Рис. 4.7.6

Собственные векторы матрицы момента инерции выглядят так же, как собственные векторы, диагонализующие квадратичную форму. Ниже дается пример двух матриц с одним набором собственных векторов, но с разными собственными значениями.

Можно найти векторы, которые одновременно являются собственными векторами матрицы А (что определяет геометрическую форму) и моментами инерции матрицы MI. Эти матрицы различны и имеют разные собственные значения. Согласно теореме линейной алгебры, если две матрицы коммутируют (т. е. для них произведения A.B и B.A одинаковы), то для них можно найти общий набор собственных векторов.

Воспользуемся свойствами.

Вначале ортогональность собственных векторов матрицы MI:

• E1vMI.E2vMI;

• E1vMI.E1vMI;

Конечно, собственные векторы матрицы MI перпендикулярны. Попробуем получить результат, редактируя предыдущие команды.

Следующие три строки показывают, какие собственные векторы матрицы MI соответствуют главным осям эллипсоида:

• E2v.E1vMI;

• E3v.E1vMI;

Первый вектор матрицы MI соответствует второму собственному вектору квадратичной формы:

• E1v . E2vMI;

• E2v . E2vMI;

• E3v . E2vMI;

Второй вектор матрицы MI соответствует первому собственному вектору квадратичной формы:

• E1v . E3vMI;

• E2v . E3vMI;

• E3v . E3vMI;

Здесь меняется знак. Собственные векторы определяются с точностью до постоянного множителя. На рис. 4.7.7 показаны все векторы и QF.

• display(E1,E2,E3,E1MI,E2MI,E3MI,QFp,axes=boxed);

Рис. 4.7.7

Рассчитаем коммутатор двух матриц – той, что для QF, и матрицы MI.

• A . MI - MI . A;

Понятно, что произведения A.MI и MI.A – это одно и то же с точностью численного расчета.

Воспользуйтесь соотношениями между главными осями квадратичной формы и главными осями матрицы момента инерции. Рассмотрите матрицы А, для которых форма близка к круговому сечению, т. е. найдите такие матрицы A, для которых два собственных вектора очень близки.

Что произойдет с собственными векторами, если два собственных значения очень близки? Какая форма дает три одинаковых собственных значения?