Войти
ПрограммированиеФорумФизика

Аккуратный вывод JV = c

Страницы: 1 2 3 Следующая »
#0
2:00, 18 апр 2014

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

В двумерном случае всё понятно, интересует только трёхмерный. Пусть уравнение связи на движение двух тел имеет вид С(p1, R1, p2, R2) = 0; здесь p1, R1 - позиция центра масс и матрица поворота(нет, не кватернион) первого тела, p2, R2 - второго. Интересует именно общий метод, как из уравнения на позиции получить уравнение на скорости:
J*(v1, w1, v2, w2)T = 0;
Давайте посмотрим, что уже есть:
SoRRoW в статье: http://www.gamedev.ru/community/gd_physcomm/articles/?id=4889 рассматривает только ball-socket соединение, дифференцирует его по времени и замечает, что вау, получилось как раз то, что нужно. интересует общий случай, а не только ball-socket.
Suslik в статье: http://www.gamedev.ru/community/gd_physcomm/articles/phys_engine_development напрочь забил за уравнение связи и сразу пляшет от якобиана. засранец.
XperienS в статье: http://www.gamedev.ru/community/gd_physcomm/articles/LCP_derivation положил болт на матрицу поворота и дифференцирует уравнение связи по обобщённым координатам:
Изображение
Всё здорово, но у нас положение тел x описывается не какими-то волшебными обобщёнными координатами, а вполне конкретной матрицей поворота, по которой дифференцировать, ээ.. нелегко.

Более того, в предыдущей формуле используется равенство
dx/dt=v, что верно опять же лишь для обобщённых координат, а для случая
x=(p, R)T, выполняется dx/dt = (v, w*R)T, а вовсе не dx/dt = (v, w)T. более того, под произведением w*R понимается постолбцовое векторное произведение w и R, что тоже не так-то очевидно.

Kenny Erleben устанавливает связь между
dx/dt и (v, w) посредством дополнительной матрицы S:
dx/dt = S*(v, w), в данном случае S = ((1, 0), (0, R)) и якобианом на самом-то деле является вовсе не
dC/dx, а dC/dx*S, но операция дифференцирования d/dx, содержащая дифференцирование по матрице, от этого вовсе не выглядит менее магической.

#1
3:40, 18 апр 2014

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

\( 0=\partial_tC(p^a,p^b,R^a,R^b)= \)

\( =\partial_{p^a}C_i\dot{p^a_i}+\partial_{p^b}C_i\dot{p^b_i}+\partial_{R^a}C_{ij}\dot{R^a_{ij}}+\partial_{R^b}C_{ij}\dot{R^b_{ij}}= \)

\( =\partial_{p^a}C_iv^a_i+\partial_{p^b}C_iv^b_i+\partial_{R^a}C_{ij}R^a_{ik}\varepsilon_{kjl}\omega^a_l+\partial_{R^b}C_{ij}R^b_{ik}\varepsilon_{kjl}\omega^b_l. \)

Я тут мог напортачить с производными от матрицы поворота, но вроде как-то так.

P. S. Встроенный LaTeX тут какой-то кривой, точки сверху еле видно.

#2
6:14, 18 апр 2014

э? ну так что с этим делать-то? во-первых, что ты назвал символом эпсилон в dR/dt? во-вторых, для меня по-прежнему не очевидно, каким образом это сведётся к какому-нибудь известному якобиану вроде ball-socket соединения.

#3
15:44, 18 апр 2014

Эпсилон — это абсолютно антисимметричный тензор, из него делаются векторные произведения: \([a\times b]_i=\varepsilon_{ijk}a_jb_k\).

В общем, я тут проверил, более правильным будет другой порядок множителей:

\( \partial_{p^a}C_iv^a_i+\partial_{p^b}C_iv^b_i+\partial_{R^a}C_{ij}\varepsilon_{ikl}\omega^a_lR^a_{kj}+\partial_{R^b}C_{ij}\varepsilon_{ikl}\omega^b_lR^b_{kj}=0. \)

Теперь для ball-socket.
Исходная связь:

\( C=D_i[p^a_i-p^b_i+R^a_{ij}r^a_j-R^b_{ij}r^b_j], \)

где D — произвольный вектор, а r — смещение точки соединения от центра масс.
Подставляем в решение:

\( D_iv^a_i-D_iv^b_i+D_ir^a_j\varepsilon_{ikl}\omega^a_lR^a_{kj}-D_ir^b_j\varepsilon_{ikl}\omega^b_lR^b_{kj}=0. \)

Выкидываем D и перегруппировываем:

\( v^a_i-v^b_i+\varepsilon_{ikl}R^a_{kj}r^a_j\omega^a_l-\varepsilon_{ikl}R^b_{kj}r^b_j\omega^b_l=0. \)

В векторных обозначениях:

\( \mathbf{v}^a-\mathbf{v}^b+[R^a\mathbf{r}^a\times\mathbf{\omega}^a]-[R^b\mathbf{r}^b\times\mathbf{\omega}^b]=0. \)

Вроде то что надо.

#4
18:46, 18 апр 2014

}:+()___ [Smile]
гмм.. да, похоже на то. в тензорном виде действительно проще получилось.

Прошло более 10 месяцев
#5
4:17, 18 фев 2015

пишу диссер, оформляю формулы. мне кажется, что длинную формулу из #3 должно быть как-то возможно оформить в матричном виде или просто более наглядном.
итак, мои рассуждения следующие. я просто опустил все линейные члены и член с поворотом для второго тела и положил \(C(p^a,p^b,R^a,R^b)=C(R)\), просто для краткости. далее расписываем:
\(\frac{dC(R)}{dt}= \frac{dC}{dR}\frac{dR}{dt}= \frac{dC}{dR}(\omega \otimes R) \)
цель состоит в том, чтобы привести это к виду:
\(\frac{dC(R)}{dt}=J\omega\)
очевидно, что это возможно, так как это возможно в тензорном виде, но хочется выразить \(J\)именно как матрицу. мои действия:
\((\omega \otimes R)_{ij}=\eps_{ipq}\omega_p R_{qj}\) смайл, обрати внимание, что индексы у \(\eps\) наоборот. у тебя в моих обозначениях вместо \(\eps_{ipq}\) стоит \(\eps_{iqp}\), дальше ты ещё раз теряешь знак, когда снова из тензорного вида получаешь тензорное произведение. это я что-то не заметил или у тебя баг? продолжаем:
\(\left(\frac{dC(R)}{dt}\right)_{ij}=\)
\(\left(\frac{dC}{dR}\right)_{ij}\left(\frac{dR}{dt}\right)_{ij}=\)
\(\left(\frac{dC}{dR}\right)_{ij}\eps_{ipq}\omega_p R_{qj}=J_p\omega_p\)
то есть в принципе
\(J_p=\eps_{ipq}\left(\frac{dC}{dR}\right)_{ij} R_{qj}=\eps_{ipq}\left(\left(\frac{dC}{dR}\right)R^T\right)_{iq}\)
но что это такое? что означает произведение \(\eps_{ijk}A_{jk}\), как это записывается в матричном виде?

насколько я понимаю, запись вида \(\eps_{ipq}A_{pj}B_{jq}\) можно трактовать как сумму векторных произведений - 1-й столбец \(A\) на 1-ю строку \(B\) + 2-й столбец \(A\) н 2-й столбец \(B\) итд. но это какое-то отстойное объяснение и скорее следствие из того, чем этот объект реально является.

#6
6:37, 18 фев 2015

Suslik
> смайл, обрати внимание, что индексы у [cht]\eps[/cht] наоборот. у тебя в моих обозначениях вместо [cht]\eps_{ipq}[/cht] стоит [cht]\eps_{iqp}[/cht], дальше ты ещё раз теряешь знак, когда снова из тензорного вида получаешь тензорное произведение. это я что-то не заметил или у тебя баг?
Это все с точностью до определения векторных произведений (которые меняют знак при отражениях пространства), там надо все честно выписывать и проверять.

> что означает произведение [cht]\eps_{ijk}A_{jk}[/cht], как это записывается в матричном виде?
Это вектор, дуальный антисимметричной части A. По-моему, нормальных матричных обозначений для этого нету.
В принципе, можно ввести оператор дуальности и выражать через него (особенно полезно в 2D, ибо там это оператор поворота вектора на 90°).
Я поэтому тензорные обозначения и люблю, что они позволяют разруливать сложные случаи.

#7
6:44, 18 фев 2015

}:+()___ [Smile]
> Я поэтому тензорные обозначения и люблю, что они позволяют разруливать сложные
> случаи.
просто программируешь-то ты ведь не в тензорном виде, а в матричном, да и понимаешь тоже. поэтому всякие тензора хороши для промежуточных вычислений, но от тензоров выше второго ранга обычно нужно и возможно избавляться. но во всяких неклассических случаях действительно тензорное выражение вроде того, что у меня получилось для якобиана, вероятно, проще никак не объяснять, а принять как гештальдт.

> Это все с точностью до определения векторных произведений (которые меняют знак
> при отражениях пространства), там надо все честно выписывать и проверять.
просто ты сам выписал в начале #3 формулу для правосторонней системы координат, а остальные выкладки - насколько я понял, в левосторонней

#8
14:53, 18 фев 2015

Suslik
> но от тензоров выше второго ранга обычно нужно и возможно избавляться.
Дак обычно таких тензоров всего один (\(\varepsilon_{ijk}\)), у которого 6 ненулевых компонентов, причем единичных по модулю.
Например,
\(\varepsilon_{ijk}A_{jk}=\{A_{yz}-A_{zy},\quad A_{zx}-A_{xz},\quad A_{xy}-A_{yx}\}_i.\)
Для некоторых выражений с абсолютно антисимметричным тензором ввели векторные обозначения (векторное произведение, детерминант), для остальных — нет.

> просто ты сам выписал в начале #3 формулу для правосторонней системы координат, а остальные выкладки - насколько я понял, в левосторонней
Вполне может быть, ибо выводил на коленке. Я часто такие вещи подгоняю по факту, до получения правдоподобного результата.

#9
1:38, 20 фев 2015

Ещё вопрос. Формализм, будь он неладен. Пусть мы ввели функцию ограничения \({C}(X)=0\), продифференцировали её, получили(опустим вышерассмотренную проблему со вращением, которая, на самом деле, показывает, что всё немного не так просто):
\(\frac{dC}{dt}=\frac{dC}{dX}\frac{dX}{dt}=JV\)
И теперь нам нужно перейти к обобщённым силам(или, если проинтегрировать по времени, импульсам), возникающим в связях, которые по идее равны
\(F=J^T\lambda\), вопрос в том, как это сделать проще всего.

У XperienS'а есть статья, в которой этот вопрос рассматривается: http://www.gamedev.ru/community/gd_physcomm/articles/LCP_derivation , но, как мне кажется, у него вводится слишком много дополнительных сущностей и необходимое соотношение можно получить из более прозрачных принципов. Например, в той статье доказывается и одновременно фактически неявно используется тот факт, что \(\frac{dC}{dX}\) является направлением ограничивающей силы. Но если это - направление, то помножив на \(\lambda\), мы ведь фактически и получим саму силу, разве нет? 

Откуда по берётся принцип минимизации работы ограничивающих сил? откуда вообще берётся формула?:
Изображение
В комментах изрядный срачик по этому вопросу уже возник, ответа я так внятного и не нашёл. Если рассматриваются идеальные связи, работа которых по определению равна нулю, то откуда инфа, что минимум функции работы вообще равен 0, а не минус бесконечности, например?
Изображение

В википедии постулируется(хоть и с оговоркой), что возможное перемещение всегда ортогонально вектору ограничивающей силы: https://en.wikipedia.org/wiki/D%27Alembert%27s_principle#Derivati… special_cases то есть тоже условие, которое я пытаюсь откуда-то получить, постулируется.

Короче нужна помощь - опять ответ известен, надо его как можно аккуратнее и проще(чтобы не плодить сущности) получить.

#10
2:23, 20 фев 2015

Если бы я еще че-то помнил из механики... По идее, то, что обзывается "принципы", — это аксиомы.
В общем, нужны точные определения и утверждения, из которых можно исходить, а иначе я не понимаю постановку задачи.

#11
2:32, 20 фев 2015

}:+()___ [Smile]
собственно, в этом и вопрос - что принять за аксиомы, чтобы требуемые соотношения было проще всего получить? XperienS в статье утверждает, что необходимо принять за аксиому вариационный принцип Д'Аламбера, он, как мне кажется, тут как ослу пятая нога и вполне можно ограничиться просто вторым законом Ньютона и определением идеальной связи(работа возникающих в ней элементарных сил на любом возможном перемещении равна нулю).

#12
3:17, 20 фев 2015

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

#13
4:15, 20 фев 2015

Я более-менее разобрался со всем, кроме принципа минимизации работы идеальных связей:
Изображение
Откуда он берётся? Почему из того, что идеальные связи не совершают работу, следует, что работа, совершаемая ими, минимальна?

#14
8:39, 20 фев 2015

Suslik
> Почему из того, что идеальные связи не совершают работу, следует, что работа, совершаемая ими, минимальна?
Это, вообще, походу, какой-то бред. Что как обозначено, непонятно, и, в конечном счете, все равно приходит к \(F_c=J^T\lambda\).
К тому же самому можно сразу прийти из принципа наименьшего дейтсвия: при нахождении минимума добавляются стандартным способом множители Лагранжа, которые при варьировании дают силоподобный член, обзываемый силами реакции.

\(\min\int L(r,\dot r)dt,\)

\(L\to \tilde L=L+\lambda C(r),\)

\(\delta\tilde L=\frac{\partial L}{\partial r}\delta r-\frac d{dt}\left(\frac{\partial L}{\partial\dot r}\right)\delta r+\lambda\frac{\partial C}{\partial r}\delta r,\)

\(\frac{\partial L}{\partial r}=F_e,\qquad\qquad\qquad\qquad\lambda\frac{\partial C}{\partial r}=\lambda J=F_c.\)

Т. е., можно сказать, что выражение через якобиан — это, практически, определение.

Страницы: 1 2 3 Следующая »
ПрограммированиеФорумФизика

Тема в архиве.