Главная
Магия тензорной алгебры: Часть 9 — Вывод тензора угловой скорости через параметры конечного поворота. Применяем голову и Maxima

Магия тензорной алгебры: Часть 9 — Вывод тензора угловой скорости через параметры конечного поворота. Применяем голову и Maxima +25

16.07.2015 20:06

maisvendoo 7 6241 Источник

Содержание

Введение

Утекло уже порядком времени, как я обещал получить тензор угловой скорости твердого тела, выразив его через параметры конечного поворота. Если взглянуть на КДПВ, то станет понятно, почему я так долго думал — стопка бумаги на столе, это ход моих мыслей.

Преобразование тензорных выражений то ещё удовольствие...

Жестокие тензоры не хотели упрощаться. Вернее, они то хотели, но при преобразованиях, раскрытии скобок, в силу невнимательности возникали мелкие ошибки, которые не позволяли взглянуть на картину в целом. В итоге результат таки был получен. Не последнюю роль в этом сыграла СКА Maxima, которой я обратился, во многом благодаря статье пользователя EugeneKalentev. Акцент упомянутой статьи смещался в сторону вычислительной работы с тензорами, компоненты которых представлены конкретными структурами данных. Меня же интересовал вопрос работы с абстрактными тензорами. Оказалось, что Maxima может с ними работать, хоть и нет так, как может быть хотелось, но всё же она серьезно упростила мне жизнь.

Итак, мы возвращаемся к механике твердого тела, а заодно посмотрим, как работать с тензорами в Maxima.

1. Немного о конечном выражении тензора поворота

В прошлой статье мы, изучив принципы свертки произведения тензоров Леви-Чивиты, получили вот такое выражение для тензора поворота

$B_k^{\,m} = u^{\,m} \, g_{jk} \, u^{\,j} - \cos\varphi \, \left( \delta_{j}^{\,m} \, \delta_{k}^{\,q} - \delta_{k}^{\,m} \, \delta_{j}^{\,q} \right) \, u_{q} \, u^j + \sin\varphi \, g^{mi} \, U_{ik} \right \quad (1)$

Его можно упростить ещё, для этого поработаем с выражением в скобках во втором слагаемом

$\left( \delta_{j}^{\,m} \, \delta_{k}^{\,q} - \delta_{k}^{\,m} \, \delta_{j}^{\,q} \right) \, u_{q} \, u^j = \delta_{j}^{\,m} \, \delta_{k}^{\,q} \, u_{q} \, u^j - \delta_{k}^{\,m} \, \delta_{j}^{\,q} \, u_{q} \, u^j =$

$= u^m \, u_{k} - \delta_{k}^{\,m} \, u^q \, u_{q} = u^m \, u_{k} - \delta_{k}^{\,m} \quad (2)$

А все потому, что свертка тензора с дельтой Кронекера приводит к замене немого индекса тензора на свободный индекс дельты Кронекера, например $\delta_{k}^{\,q} \, u_{q} = u_{k}$ . Действуя таким образом, и учтя, что вектор $\vec{u}$ имеет длину равную единице, а значит свертка $u^{q} \, u_{q} = 1$ , мы и получаем (2). Тогда выражение тензора поворота выглядит ещё проще

$B_k^{\,m} = \left(1 - \cos\varphi \right )u^{\,m} \, u_{\,k} + \cos\varphi \, \delta_{k}^{\,m} + \sin\varphi \, g^{mi} \, U_{ik} \right \quad (3)$

как раз то выражение, которое я долго и нудно пытался вывести в седьмой статье. Что ещё раз подтверждает правило, что ничего нельзя изучать наполовину…

Для чего я это всё проделал. Во-первых, окончательно реабилитироваться в глазах читателей. Во-вторых — чтобы нам было с чем сравнивать результаты, которые мы чуть ниже получим с помощью Maxima. Ну и в третьих, чтобы сказать, что на выражении (3) наша сладкая жизнь и заканчивается, а начинается ад монструозных преобразований.

2. Применяем Maxima для вывода тензора угловой скорости

За работу с абстрактными тензорами в Maxima отвечает модуль itensor, документация по которому представлена в оригинале, в переводе на русский язык, и ещё одной книгой, основанной на документации.

Запускаем Maxima, хоть в консоли, хоть используя один из её графических фронтэндов, и пишем

kill(all);
load(itensor);
imetric(g);
idim(3);

Здесь мы чистим память (функция kill()), удаляя из неё все определения и загруженные модули, загружаем пакет itensor (функция load()), говорим, что метрический тензор будет именоваться g (функция imetric()), а так же, обязательно указываем размерность пространства (функция idim()), ибо по-умолчанию СКА считает, что работает в 4-мерном пространстве-времени с метрикой Римана. Мы с вами работаем в трехмерном пространстве классической механики, а метрика у нас любая невырожденная без кручения.

Вводим тензор поворота

B:ishow( u([],[l])*g([j,k],[])*u([],[j]) - 
cos(phi)*'levi_civita([],[l,q,i])*u([q],[])*'levi_civita([i,j,k],[])*u([],[j]) + 
sin(phi)*g([],[l,i])*'levi_civita([i,j,k])*u([],[j]))$

В itensor тензоры декларируются идентификатором, после которого в скобках указываются через запятую списки ковариантных и контравариантных индексов. Функция levi_civita() задает тензор Леви-Чивиты, а штрих перед ней означает, что данный тензор не надо вычислять. Maxima известна своей манерой вычислять и упрощать выражения если это возможно, и если штрих не поставить, то тензор Леви-Чивиты превратится в «тыкву», а конкретно, будет предпринята попытка определить его через обобщенную дельту Кронекера, что в наши планы не входит.

Символ $ является альтернативой точки с запятой и запрещает вывод на экран результатов преобразований. Для отображения информации будем использовать ishow(). Подставляемое в неё выражение выводится на экрана в привычной индексной нотации записи тензоров. На экране это выглядит так

Причем вводим мы не упрощенную его форму, с двойным векторным произведением, не вводя антисимметричного тензора $\mathbf{U}$ , с тем чтобы на первом этапе получить (3) и проверить, как машина упрощает тензорные выражения. В человеческом виде это выражение выглядит так

$g_{j\,k}\,u^{l+j}-\varepsilon_{i\,j\,k}\,\cos \Phi\,\varepsilon^{i \,l\,q}\,u_{q}\,u^{j}+g^{i\,l}\,\varepsilon_{i\,j\,k}\,\sin \Phi\,u ^{j}$

правда сумма индексов в первом слагаемом наводит на подозрения в какой-то особой магии, которая нам неведома. Но данный вывод сгенерирован непосредственно Maxima, с помощью заклинания

load(tentex);
tentex(B);

что позволяет получить вывод в виде кода LaTeX. Из недостатков — шаманство с индексами и заглавная буква, обозначающая угол поворота, но и то и другое в принципе исправимо, а порадовала меня Maxima тем, что генерирует TeX-вывод не такой избыточный, как например Maple.

Идем дальше. Дальше вспоминаем выражение, определяющее нам угловую скорость через тензор поворота

$\Omega_{\,mk} = g_{\,mp} \, B_{l}^{'\,p} \, \dot{B}_{k}^{\,l} \quad (4)$

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

Прежде всего упростим введенный тензор поворота

B:ishow(expand(lc2kdt(B)))$

Функция lc2kdt() предназначена специально для упрощения выражений содержащих тензор Леви-Чивиты. Она старается свернуть этот тензор там где это возможно, давая на выходе комбинацию сумм и произведений дельт Кронекера. Так выглядит результат

К результату lc2kdt() применяется так же функция expand(), раскрывающая скобки. Без этого Maxima выполняет свертку очень неохотно.

Теперь попытаемся вычислить полученное выражение, выполнив свертку. Функция contract() один из наиболее надежных способов выполнения свертки

B01:ishow(contract(B))$

На выходе получаем (там где это возможно, буду приводить LaTeX-вывод, генерируемый Maxima с помощью функции tentex())

$\delta_{k}^{l}\,\cos \Phi\,u_{q}\,u^{q}-u_{k}\,\cos \Phi\,u^{l}+u_{ k}\,u^{l}+g^{i\,l}\,\varepsilon_{i\,j\,k}\,\sin \Phi\,u^{j}$

что аналогично выражению (3). Это убедило меня в корректности работы программы и правильности собственных действий. Значит можно двигаться дальше.

Единственное, чего не учитывает Maxima, что вектор, вокруг которого производится поворот имеет длину равную единице. Поэтому он не сворачивает выражение $u^{q} \, u_{q}$ в скаляр. Как ему об этом сказать я не придумал и в документации этого не вычитал, хотя искал тщательно. В дальнейшем это даст серьезные трудности, о которых мы поговорим, а пока я выкрутился с помощью костыля, заменив неудобные тензоры единицей

B01:ishow(subst(1, u([],[q]), B01))$
B01:ishow(subst(1, u([q],[]), B01))$

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

B10:ishow(subst(-phi, phi, B01))$

Раз мы собираемся брать производную, то Maxima должна знать, какие тензоры и числовые параметры зависят от времени. Указываем, что от времени зависит направление оси поворота и угол поворота

depends([u,phi], t);

Чтобы наши тензоры соответствовали формуле (4), выполним переименование индексов

B10:ishow(subst(p,l,B10))$
B10:ishow(subst(l,k,B10))$
B10:ishow(subst(i1,i,B10))$
B10:ishow(subst(j1,j,B10))$
B10:ishow(subst(s,q,B10))$

Тем самым мы поменяли имена свободных индексов, для корректного умножения в (4), и переопределили имена немых индексов. По правилам положено, чтобы имена немых индексов в перемножаемых тензорах не совпадали.

Теперь берем производную по времени от тензора поворота

dBdt:ishow(diff(B01,t))$

получая на выходе

Убеждаемся, что дифференцирование выполнилось корректно. Ну и наконец вводим формулу (4), применяя последовательно раскрытие скобок и упрощение с учетом тензора Леви-Чивиты

exp1:ishow(lc2kdt(expand(g([m,p],[])*B10*dBdt)))$

То, что получается в итоге, страшно показывать, но извольте

Понять что-либо в этой каше трудно. Но видно, что вылезли комбинации дельт Кронекера (функция kdelta() задает дельту Кронекера любого требуемого ранга) со странным обозначением немых индексов, содержащих номер после %. Дело в том, что Maxima, при преобразованиях нумерует немые индексы. Попробуем упростить всё это

Во-первых, ещё раз прогоним выражение через упрощение с учетом тензора Леви-Чивиты (lc2kdt()). Затем выполним свертку (contract()). После чего попытаемся упростить нашего «крокодила», применив к нему функцию canform(), выполняющую нумерацию немых индексов и упрощение тензора. Эта функция рекомендуется разработчиками для выполнения упрощений

exp2:ishow(canform(contract(lc2kdt(exp1))))$

В итоге наблюдаем серьезное похудение «крокодила»

Но! В первом же слагаемом мы видим векторное произведение орта оси поворота самого на себя, а оно должно быть равно нулю. Maxima этого пока не поняла, надо ей указать на возможность подобного упрощения. Делаем это конструкцией

exp3:ishow(canform(contract(expand(applyb1(exp2,lc_l,lc_u)))))$

Функция applyb1() задает правила упрощения подвыражений, входящих в упрощаемое выражение. В качестве аргументов передаются выражение и список правил. Правил у нас два: lc_l и lc_u — правила преобразования подвыражений с символом Леви-Чивиты с нижними (lc_l) и верхними (lc_u) индексами. При этом мы снова раскрываем скобки, которые могут появится после преобразования подвыражений, выполняем свертку и упрощение. В итоге наблюдаем ещё одно сокращение крокодиловой массы

Наверное это ещё не предел. Но меня удивило, вероятно по незнанию и недостатку опыта следующее

Обратим внимание на переменную kdelta. Эксперименты с Maxima позволили выяснить что это — след дельты Кронекера 2 ранга, равный размерности пространства (след единичной матрицы). В нашем случае это число «3». Именно тройка должна стоять на месте kdelta. Но почему-то не стоит, возможно как-то неправильно настроены переменные конфигурации тензорного пакета. Если принять kdelta равной трем, то образуется куча разнознаковых подобных слагаемых, которые в сумме дают ноль.
Все свертки вида
$u^{i} \, u_{i}$
это модуль орта поворота, а он не меняется и равен единице. Как сказать об этом Maxima для меня пока не выяснено.
Вытекает из предыдущей проблемы, ибо

$u^{i} \, u_{i} = 1$

а значит

$\dot u^{i} \, u_{i} + u^{i} \, \dot u_{i} = 0$

откуда, в силу коммутативности операции скалярного умножения имеем

$\dot u^{i} \, u_{i} = 0, \quad u^{i} \, \dot u_{i} = 0$

Последние выражения часто встречаются в вышеприведенном результате.

Я был бы крайне признателен за подсказку знающих по перечисленным проблемам, ибо изучение документации пока не пролило свет на из решение. Функция замены subst() тут не срабатывает.

В связи с этим я снова взялся за ручку и бумагу, чтобы «доупрощать» тензор угловой скорости. Но Maxima существенно облегчила мне задачу, за что ей спасибо.

3. С помощью ручки, бумаги, напильника и бубна...

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

$\Omega_{\,km} = \left(1 - \cos\varphi \right )\left(\dot u_{\,k} \, u_{\,m} - u_{\,k} \, \dot u_{\,m}\right ) + \sin\varphi \, \left(1 - \cos\varphi \right ) \, u^{\,i} \left( \varepsilon_{\,ijk} \, \dot u^{\,j} \, u_{\,m} - \varepsilon_{\,ilm} \, \dot u^{\,l} \, u_{\,k} \right ) +$

$+ \sin\varphi \cos\varphi \, \dot u^{\,j} \varepsilon_{\,jkm} + \dot\varphi \, \varepsilon_{\,ikm} \, u^{i} \quad (5)$

Это выражение до боли напоминает, аналогичное, полученное для ортогонального базиса и выраженное в матричном виде

Погорелов Д. Ю. Введение в моделирование динамики систем тел. Стр. 31.

Тензор (5) антисимметричный, переставим индексы, меняя знаки соответствующих выражений. Первое слагаемое — разность ковариантных тензорных произведений, второе есть транспонированное первое, а такая сумма, дает на выходе антисимметричный тензор. В остальных слагаемых содержится тензор Леви-Чивиты, меняющий знак при перестановке индексов

$\Omega_{\,mk} = \left(1 - \cos\varphi \right )\left(\dot u_{\,m} \, u_{\,k} - u_{\,m} \, \dot u_{\,k}\right ) + \sin\varphi \, \left(1 - \cos\varphi \right ) \, u^{\,i} \left( \varepsilon_{\,ilk} \, \dot u^{\,l} \, u_{\,m} - \varepsilon_{\,ijm} \, \dot u^{\,j} \, u_{\,k} \right ) +$

$- \sin\varphi \cos\varphi \, \dot u^{\,j} \varepsilon_{\,jmk} - \dot\varphi \, \varepsilon_{\,imk} \, u^{i}$

И еще раз переставим индексы в последних двух слагаемых

$\Omega_{\,mk} = \left(1 - \cos\varphi \right )\left(\dot u_{\,m} \, u_{\,k} - u_{\,m} \, \dot u_{\,k}\right ) + \sin\varphi \, \left(1 - \cos\varphi \right ) \, u^{\,i} \left( \varepsilon_{\,ilk} \, \dot u^{\,l} \, u_{\,m} - \varepsilon_{\,ijm} \, \dot u^{\,j} \, u_{\,k} \right ) +$

$+ \sin\varphi \cos\varphi \, \varepsilon_{\,mjk} \, \dot u^{\,j} + \dot\varphi \, \varepsilon_{\,mik} \, u^{i} \quad (6)$

чтобы их свертка давала кососимметричные матрицы, привычные для представления векторного произведения в матричном виде.

Дальнейшее упрощение (6) возможно, если принять базис декартовым. Тогда ковариантные и контравариантные компоненты совпадут, ненулевые элементы тензора Леви-Чевиты по модулю станут равны единице, и (6) можно еще немного упростить. Мы же получили это выражение для произвольного базиса.

Заключение

Вот так, благодаря тензорным средствам Maxima я наконец разобрался с задачей выражения тензора угловой скорости через параметры конечного поворота. А заодно показал читателю живой пример работы с тензорами в СКА.

В следующий раз мы получим из (6) псевдовекторы угловой скорости и углового ускорения и вплотную приблизимся к тензорному описанию кинематики твердого тела.

Благодарю за внимание!

Продолжение следует...

Комментарии (7)

PapaBubaDiop
17.07.2015 07:57
#8503275
И все таки хотелось бы картинки формулы тензора угловой скорости, зависящей от углов поворота в декартовой системе координат. Я выписал ее лет 20 назад на бумажку из какой-то статьи, но со временем чернила стерлись, одна надежда на вас.
1. maisvendoo Автор
  17.07.2015 08:24
  #8503295
  зависящей от углов поворота в декартовой системе координат
  
  Конкретный вид матрицы зависит от того какие это углы. Если речь идет о поворотах вокруг координатных осей, то можно глянуть здесь.
  
  А можно получить непосредственно из (3), например для оси Z. Формула (3) в матричном виде
  
  $\mathbf{B} = \left(1 - \cos\varphi \right ) \mathbf{u} \, \mathbf{u}^{\,T} + \cos\varphi \, \mathbf{E} + \sin\varphi \, \mathbf{U}$
  
  Орт оси поворота и порождаемая им кососимметричная матрица
  
  $\mathbf{u} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \quad \mathbf{U} = \begin{bmatrix} 0 && -1 && 0 \\ 1 && 0 && 0 \\ 0 && 0 && 0 \end{bmatrix}$
  
  Матричное произведение столбца на строку
  
  $\mathbf{u} \, \mathbf{u}^{\,T} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \begin{bmatrix} 0 && 0 && 1 \end{bmatrix} = \begin{bmatrix} 0 && 0 && 0 \\ 0 && 0 && 0 \\ 0 && 0 && 1 \end{bmatrix}$
  
  Складываем всё вместе
  
  $\mathbf{B_z} = \begin{bmatrix} 0 && 0 && 0 \\ 0 && 0 && 0 \\ 0 && 0 && 1 - \cos\varphi \end{bmatrix} + \begin{bmatrix} \cos\varphi && 0 && 0 \\ 0 && \cos\varphi && 0 \\ 0 && 0 && \cos\varphi \end{bmatrix} + \begin{bmatrix} 0 && -\sin\varphi && 0 \\ \sin\varphi && 0 && 0 \\ 0 && 0 && 0 \end{bmatrix} =$
  
  $= \begin{bmatrix} \cos\varphi && -\sin\varphi && 0 \\ \cos\varphi && \sin\varphi && 0 \\ 0 && 0 </p></div><br><ol class=$
2. PapaBubaDiop
  17.07.2015 09:11
  #8503331
  -1
  Эту матрицу бе-Зет и второклассник знает наизусть, тем более она у вас с ошибкой написана))
  А вот вокруг вектора на произвольный угол нужна на хабре, чтобы бумажку не искать.
  1. PapaBubaDiop
    17.07.2015 09:13
    #8503335
    -1
    И пожалуйста тензор деформации в углах поворота фи и пси.
  1. maisvendoo Автор
    17.07.2015 09:14
    #8503339
    вокруг вектора на произвольный угол
    
    Вот она
    $\mathbf{B} = \left(1 - \cos\varphi \right ) \mathbf{u} \, \mathbf{u}^{\,T} + \cos\varphi \, \mathbf{E} + \sin\varphi \, \mathbf{U}$
    
    c ошибкой
    
    Да, перепутал при наборе
    $= \begin{bmatrix} \cos\varphi && -\sin\varphi && 0 \\ \sin\varphi && \cos\varphi && 0 \\ 0 && 0 && 1 \end{bmatrix}$

maisvendoo Автор
17.07.2015 08:24
#8503291
del

PavelSandovin
17.07.2015 10:15
#8503423
>> Жестокие тензоры не хотели упрощаться. Вернее, они то хотели, но при преобразованиях, раскрытии скобок, в силу невнимательности возникали мелкие ошибки, которые не позволяли взглянуть на картину в целом.

Вы не пробовали нотацию Пенроуза для тензорных вычислений? Она изложена в дополнении к книге Пенроуз Р., Риндлер В. Спиноры и пространство-время (в 2-х томах). Т.1: Два-спинорное исчисление и релятивистские поля.

Выглядит как пляшущие человечки, но, говорят, очень хороша для ручных вычислений.

Магия тензорной алгебры: Часть 9 — Вывод тензора угловой скорости через параметры конечного поворота. Применяем голову и Maxima +25

Содержание

Введение

1. Немного о конечном выражении тензора поворота

2. Применяем Maxima для вывода тензора угловой скорости

3. С помощью ручки, бумаги, напильника и бубна...

Заключение

Комментарии (7)

PapaBubaDiop

maisvendoo Автор

PapaBubaDiop

PapaBubaDiop

maisvendoo Автор

maisvendoo Автор

PavelSandovin