4. XPBD и SmallSteps

Получаем уравнения

Ограничения действуют силами $f$. У этих ограничений есть потенциальная энергия $U$.

$$ f = -\frac{\partial U(x)}{\partial x} $$

$$ M \ddot{x} = f $$

M – это матрица массы.

$$ M \ddot{x} = -\frac{\partial U(x)}{\partial x} $$

Теперь используют неявную численную схему:

$$ \begin{equation} M \frac{x^{n+1} - 2x^n + x^{n-1}}{{\Delta t}^2} = -\frac{\partial U(x^{n+1})}{\partial x} \end{equation} $$

Эта схема неявная, так как $x^{n+1}$ зависит от U, которое зависит от $x^{n+1}$. Для почти всех систем неявная схема абсолютно устойчива.

Вводим потенциал

Теперь нужно ввести $U$ для ограничений. Мы считаем что уравнение ограничения $C(x) = 0$. зависит только от $x$. А еще мы считаем, что оно дифференцируемо. Вообще можно в качестве потенциала использовать любую функцию, но нам полезно добиться хороших свойств от него.

• $U(x) = 0$ если $C(x) = 0$ Т.е когда уравнение ограничения выполняется, потенциал равен 0 и сира ограничения тоже равна нулю. Что логично.
• При отклонении от $C(x) = 0$ потенциал растет. Ну и сила ограничения тоже будет расти.

Для этого хорошо подходит квадратичный потенциал:

$$ U(x) = \frac{1}{2} С(x)^T * \alpha^{-1} * C(x) $$

Здесь $C(x)$ это вектор ограничений, а $\alpha$ это матрица compliance. Вводится термин compliance $\alpha$: Это параметр обратный жесткости. Оно в пределах $(0; +\inf)$ Чем он больше тем мягче ограничение. Эта матрица может быть не только диагональной. Благодаря этому можно комбинаю жесткости разных ограничений.

Дальше вводим переменную $\lambda$ которая будет множителем Лагранжа.

$$ \lambda = -\tilde{\alpha} * C(x) $$

Это уравнение можно переписать вот так:

$$ C(x) + \tilde{\alpha} * \lambda = 0 $$

Тут $\tilde{\alpha}$ это $\frac{\alpha} {\Delta t^2}$.

Уравнения которые нужно решить

Подставляем уравнение силы и $\lambda$ в уравнение движения (1) и получаем финальную систему, которую нужно решить:

$$ \begin{equation} \begin{aligned} M(x^{n+1} - \tilde{x}) - \Delta t^2 \nabla C(x^{n+1}) *\lambda^{n+1} = 0 \\ C(x^{n+1}) + \tilde{\alpha} * \lambda^{n+1} = 0 \end{aligned} \end{equation} $$

Где $\tilde{x} = x^n + \Delta t * v^n$

Решаем

Получили нашу любимую нелинейную систему уравнений. Нам нужно найти  $x^{n+1}$ и $\lambda^{n+1}$. В общем виде ее решить невозможно. Поэтому в статье предлагают решать ее итерационным методом Ньютона. Обозначим уравнения из (2) как

$$ \begin{aligned} g(x^{n+1}, \lambda^{n+1}) = 0 \\ h(x^{n+1}, \lambda^{n+1}) = 0 \end{aligned} $$

Линеаризовываем их относительно $x^{n+1}$ и $\lambda^{n+1}$:

Раньше мы писали ${n+1}$, и это означало значение на следующем шаге. Мы сейчас будем искать $x^{n+1}$ и $\lambda^{n+1}$ итеративно. И поэтому будем писать $x^{i}$ и $\lambda^{i}$, где $i$ – номер итерации. Т.е. $x^{i}$ это приближение к $x^{n+1}$ на $i$-й итерации.

$$ \begin{bmatrix} K, &\nabla C^{T}(x^{i}) \\ \nabla C{(x^{i})}, &\tilde{\alpha} \end{bmatrix} \begin{bmatrix} \Delta x^{i} \\ \Delta \lambda^{i} \end{bmatrix} = \begin{bmatrix} g(x^{i}, \lambda^{i}) \\ h(x^{i}, \lambda^{i}) \end{bmatrix} $$

Где

$$ \begin{aligned} &K = \partial{g(x^{n+1}, \lambda^{n+1})} / \partial{x^{n+1}} = \\ &= M - \Delta t^2 \nabla C(x^{n+1}) \nabla C^{T}(x^{n+1}) \\ \end{aligned} $$

Из этого уравнения мы получаем $\Delta x^{i}$ и $\Delta \lambda^{i}$, которые мы используем для улучшения наших приближений:

$$ \begin{aligned} x^{i+1} = x^{i} + \Delta x^{i} \\ \lambda^{i+1} = \lambda^{i} + \Delta \lambda^{i} \end{aligned} $$

И так делаем пока не сойдемся. В качестве начального приближения $x^{0} = \tilde{x}$ и $\lambda^{0} = 0$.

Проблема этого метода в том, что его не всегда просто реализовать. Нужно расчитывать матрицу $K$ у которой второе слагаемое – Гессиан функции $C(x)$. Это сложно.

Упрощаем

• Гессиан ограничений дает поправку как $O(\Delta t^2)$. Это мало. Мы можем его просто выкинуть. Тогда $K = M$. Это константная матрица масс, которая известна с самого начала.
• $g(x^{i}, \lambda^{i}) = 0$. Для начальных условий $g(x^{0}, \lambda^{0}) = 0$. Остальные итерации изменяют $g$ на малую величину, которая зависит от изменения градиента ограничений. Если ограничение линейное, то $g(x^{i}, \lambda^{i}) = 0$ всегда.

И вот мы получаем простую систему уравнений:

$$ \begin{bmatrix} M, &\nabla C^{T}(x^{i}) \\ \nabla C{(x^{i})}, &\tilde{\alpha} \end{bmatrix} \begin{bmatrix} \Delta x^{i} \\ \Delta \lambda^{i} \end{bmatrix} = \begin{bmatrix} 0 \\ h(x^{i}, \lambda^{i}) \end{bmatrix} $$

Переписав уравнения в удобной форме получаем:

$$ \begin{equation} \begin{bmatrix} \nabla C(x^{i}) M^{-1} \nabla C^{T}(x^{i}) + \tilde{\alpha} \end{bmatrix} \Delta \lambda = -C(x^{i}) - \tilde{\alpha} \lambda^{i} \end{equation} $$

$$ \begin{equation} \Delta x = M^{-1} \nabla C^{T}(x^{i}) \Delta \lambda \end{equation} $$

Дампинг здесь вводится через потенциал диссипации энергии.

$$ D(x, v) = \frac{1}{2} \dot{C}(x)^T \beta \dot{C}(x) = \frac{1}{2} v^T*\nabla C(x)^T \beta \nabla C(x) * v $$

Где $\beta$ это матрица дампинга.

Сила дампинга равна:

$$ f = -\frac{\partial D(x, v)^T}{\partial v} = -\nabla C(x)^T \beta \nabla C(x) * v $$

Точно также здесь можно ввести переменную $\lambda_{damp}$:

$$ \lambda_{damp} = -\tilde{\beta} \nabla C(x) * v $$

Где $\tilde{\beta} = \frac{\beta}{\Delta t^2}$

Дальше мы делаем очень крутой переход. Сила дампинга должна быть сонаправлена с силой самого ограничения. Поэтому мы можем просто сложить $\lambda_{elastic}$ и $\lambda_{damp}$ и получить общую силу ограничения.

$$ \lambda = \lambda_{elastic} + \lambda_{damp} = -\tilde{\alpha} \nabla C(x) - \tilde{\beta} \nabla C(x) * v $$

$$ h(x, \lambda) = C(x) + \tilde{\alpha} \lambda + \tilde{\alpha}\tilde{\beta} \nabla C(x) * v = 0 $$

Подставляя это в уравнение (3) и линеаризуя получаем:

$$ \left[\mathbf{I} + \frac{\tilde{\alpha} \tilde{\beta}}{\Delta t} \nabla \mathbf{C}(\mathbf{x}_i) \mathbf{M}^{-1} \nabla \mathbf{C}(\mathbf{x}_i)^{\mathrm{T}} + \tilde{\alpha}\right] \Delta \boldsymbol{\lambda} = -\mathbf{h}(\mathbf{x}_i, \boldsymbol{\lambda}_i) $$

И в итоге получаем вот такое решение для $\Delta \lambda$: