Войти
Физика для игрСтатьи

Итеративные LCP-солверы. Прочёл Xperiens на gamedev_lecture.

Автор:

Disclaimer: лог почищен, некоторые опечатки поправлены. полный лог

[18:04] <XperienS> значит, пару слов на введение
[18:05] <XperienS> тема для тех, кто уже определился что он хочет делать физику "на ограничениях"
[18:05] <XperienS> пару доков в тему:
[18:06] <XperienS> Kenny Erleben "Stable, Robust, and Versatile Multibody Dynamics Animation"
[18:06] <XperienS> здесь есть много инфы для размышления
[18:06] <XperienS> более сокращенный вариант без оснований -
[18:07] <XperienS> Erin Catto "Iterative Dynamics with Temporal Coherence"
[18:07] <XperienS> Гугл поможет если что
[18:07] <XperienS> значит постановка задачи
[18:07] <XperienS> у нас есть матрица A, составленная из якобиан и так называемых матриц масс
[18:08] <XperienS> еще есть вектор RHS ("right hand side") который составляется в зависимости от ограничений действующих на тело (пару тел)
[18:09] <XperienS> нам нада узнать лямбду (для краткости это икс), чтобы вычислить силу, которая это ограничение и осуществляет
[18:10] <XperienS> формула щас не так важна потому как лекция про солверы )
[18:10] <XperienS> получается уравнение вида A*x = RHS
[18:10] <XperienS> что нужно получить, это итеративное решение для х
[18:11] <XperienS> разложем матрицу A
[18:11] <XperienS> A = L + D + U
[18:11] <XperienS> где
[18:11] <XperienS> L - строго нижнетреугольная матрица
[18:11] <XperienS> D - диагональная
[18:12] <XperienS> а U, соответственно, строго верхнетреугольная матрица
[18:12] <XperienS> теперь имеем выражение (L + D + U)*x = RHS
[18:12] <XperienS> замечу что пока мы решаем простое СЛАУ, налогаемые ограничения будут разрешены позже
[18:13] <XperienS> из вышеописанной формулы легко получить D*x = RHS - (L + U)*x
[18:14] <XperienS> и наконец финальная формула x' = (D^-1)*(RHS - (L + U)*x)
[18:14] <XperienS> именно такая последовательность обусловлена тем, что матрицу D легче всего инвертировать, т.к. она у нас диагональная
[18:15] <XperienS> [18:14] <Zeux> матрица D необратима, если хоть один диагональный элемент 0. Что тогда?
[18:16] <XperienS> там получаются матрицы масс в диагонали
[18:16] <XperienS> а недвижимые тела характеризуются БЕСКОНЕЧНОЙ массой
[18:16] <XperienS> вобщем, тела с нулевой массой не поддерживаются )
[18:17] <XperienS> так вот, напомню, что инвертированная диагональная матрица это тоже диагональная матрица, элементы которой являются инверсией элементов изначальной матрицы
[18:17] <XperienS> так же следует заметить что матрица состоит из блоков и она диагональная в том случае если принимать за элемент один блок
[18:18] <XperienS> т.е. инвертированая D это просто инвертирование всех элементов-блоков
[18:19] <XperienS> полученная формула есть формула метода Якоби
[18:19] <XperienS> ща запишу более удобную для имплементации и дальнейших рассуждений формулу
[18:20] <XperienS> x'i = (RHS_i - S(Li,j * xj)[j=0,i-1] - S(Ui,j * xj)[j=i+1,n-1]) / Di,i;
[18:20] <XperienS> S - суммирование, в кв. скобках показан параметр суммирования
[18:21] <XperienS> заметно, что вычисляя x'i, мы уже вычислили все элементы до инджекса i, как бы странно это ни звучало ))
[18:22] <XperienS> поэтому для того, чтобы увеличить скорость решения (т.е. время, за которое мы подбираемся до искомого ответа с некоторой точностью)
[18:22] <XperienS> следует в S(Li,j * xj)[j=0,i-1] использовать уже вычисленные значения x
[18:23] <XperienS> формула примет вид
[18:23] <XperienS> x'i = (RHS_i - S(Li,j * x'j)[j=0,i-1] - S(Ui,j * xj)[j=i+1,n-1]) / Di,i;
[18:23] <XperienS> и название метода сменится на метод Гаусса-Сейделя
[18:24] <XperienS> (вообще-то по немецки читается Зайделя, но это отдельная тема для дискуссии - в некоторой литературе встречается смежный вариант Зейделя)
[18:25] <XperienS> забыл упомянуть что x' это новый икс на текущем шаге а х - икс, вычисленый на предыдущем шаге
[18:25] <XperienS> должен возникнуть вопрос - если это первый шаг, то чему равно это самое х?
[18:25] <XperienS> используется такое понятие как "начальное предположение"
[18:26] <XperienS> т.е. первоначальный икс мы устанавливаем в некоторое известное значение
[18:26] <XperienS> есть несколько методов
[18:26] <XperienS> первый это, естественно, установить изначальный икс в 0-вектор
[18:27] <XperienS> это самый простой способ, но он увеличивает время решения системы (увеличивает время сходимости)
[18:28] <XperienS> второй - это установить икс в конечные значения, полученные в предыдущем кадре (не итерации, а именно кадре, или временном шаге)
[18:28] <XperienS> этот метод значительно ближе к истине и не требует особых затрат умственных ресурсов
[18:29] <XperienS> последний это кэширование предыдущих результатов для определенных типов соединений
[18:29] <XperienS> возьмем например контактное соединение
[18:29] <XperienS> в прошлом кадре мы уже нашли нужную лямбду
[18:30] <XperienS> и, используя когерентность симуляции, можно использовать предыдущие значения как инишл гесс
[18:31] <XperienS> для этого мы создаем "кучу" этих лямбд а потом ищем в ней лямбду для текущего контакта например по индексам тел участвующих в нем
[18:32] <XperienS> начальное предположение позволяет уменьшить время сходимости до нужного значения с определенной точностью, как я уже говорил
[18:32] <XperienS> теперь время поговорить об этой самой сходимости
[18:32] <XperienS> формула
[18:33] <XperienS> ||x' - x|| < EPS
[18:33] <XperienS> где ||a|| - норма вектора a
[18:34] <XperienS> если модуль разницы векторов х' и x меньше чем некоторая точность EPS
[18:34] <XperienS> то можно останавливать вычисления
[18:35] <XperienS> это позволяет убрать ненужные итерации и снизить ресурсоемкость
[18:35] <XperienS> на практике достаточно вычислять |x'i - xi| < EPS для всех элементов векторов х' и x
[18:37] <XperienS> теперь о том, как можно уменьшить время сходимости
[18:37] <XperienS> можно использовать метод SOR - successive overrelaxation
[18:37] <XperienS> последовательная верхняя релаксация
[18:37] <XperienS> для того, чтобы понять че это за
[18:38] <XperienS> надо ввести понятие остаточного вектора
[18:39] <XperienS> опять же, на пальцах, остаточный вектор это вектор, являющийся разницей между x' и x
[18:39] <XperienS> т.е. x' = x + r
[18:39] <XperienS> r - остаточный вектор (residual vector)
[18:40] <XperienS> для того, чтобы подобраться быстрее к нужному решению (т.е. чтоб выполнить тест на сходимость)
[18:40] <XperienS> нужно увеличить остаточный вектор (т.е. увеличить "скачки" с которыми мы подбираемся к нужному решению)
[18:40] <XperienS> x' = x + w*r
[18:41] <XperienS> параметр w это параметр релаксации
[18:41] <XperienS> если омега больше 1, то это метод SOR
[18:42] <XperienS> если меньше 1 то это метод SUR (successive under-relaxation, увеличивает время просчета, и значительно увеличивает точность)
[18:42] <XperienS> т.е. для того, чтобы алгоритм быстрее сходился до нужной нам точности нада поставить w больше чем 1, но главное не переборщить
[18:43] <XperienS> а то алгоритм будет скакать вокруг да около нужного нам решения
[18:44] <XperienS> отмечу что без теста на сходимость апгрейд до Гаусса-Сейделя с последовательной верхней релаксации не имеет смысла
[18:44] <XperienS> т.е. прироста в скорости мы не получим если будет просто фиксированное количество итераций
[18:46] <XperienS> надеюсь по итеративным методам пока понятно )
[18:47] <XperienS> значит что мы щас имеем так это метод для решения СЛАУ
[18:47] <XperienS> итеративный
[18:47] <XperienS> называется SOR - GS
[18:47] <XperienS> теперь нам надо установить некоторые ограничения на ограничивающие силы, во как
[18:48] <XperienS> например
[18:49] <XperienS> у нас есть контактное ограничение
[18:49] <XperienS> а это значит что это ограничение может отталкивать тела
[18:49] <XperienS> а притягивать не может
[18:49] <XperienS> значит, ограничивающая сила должна быть больше нуля
[18:50] <XperienS> а из этого следует что лямбда должна быть тоже больше нуля

Страницы: 1 2 Следующая »

30 марта 2006 (Обновление: 21 янв. 2010)