Newton's Method for Implicit Integration
Solving such a large, sparse, non-linear system is complex and computationally expensive since we need to evaluate $H(x^{n})$ in equation (12), which is called the ***Hessian***. Therefore, implicit integration is rarely used in real-time applications.
(md-position-based-dynamics)=
## Position-Based Dynamics
### Methodology
The position-based approach directly modifies vertex positions instead of evaluating forces. Hence, it does not need to solve for the *Hessian*. We first revisit the problem from the energy potential and Newton's equations of motion. Then we introduce how the position-based approach works.
The system has an energy potential $U(x)$, where $x=[x_1, x_2, ..., x_n]^T$ represents the system state using vertex coordinates. The elastic energy potential $U$ is:
```{math}
U = \frac{1}{2} k \Delta x^2 \tag{13}
```
where $k$ is the spring constant representing stiffness; $\Delta x$ represents the extension, i.e., how far the vertex moves beyond the rest position. For a single constraint $C_1$, its potential energy has the following form:
```{math}
U(x) = \frac{1}{2} k C_{1} (x)^2 \tag{14}
```
For a system with $m$ constraints, it can be written in matrix form:
```{math}
U(x) = \frac{1}{2} C^T(x) \alpha^{-1} C(x) \tag{15}
```
where $C(x) = [C_{1}(x), ..., C_{m}(x)]^T$ is the constraint column vector; $\alpha$ is a block diagonal compliance matrix corresponding to the inverse of stiffness $k_1, ..., k_m$.
According to Newton's equations of motion:
```{math}
f_{elastic}^{T} = M \ddot{x} \tag{16}
```
Substituting equation (15) into (16), we have:
```{math}
M \ddot{x} = - \nabla U^T = -\nabla C^T(x) \alpha^{-1} C(x) \tag{17}
```
in which $M$ is a diagonal mass matrix of cloth vertices.
We use the implicit Euler method to perform time discretization on equation (17):
```{math}
M(\frac{x^{n+1} - 2x^n + x^{n-1} }{\Delta t^2}) = - \nabla U^T(x^{n+1}) \tag{18}
```
where $\alpha$ is a block diagonal matrix representing inverse stiffness. For more details, please refer to [Appendix: Implicit Euler Time Discretization](implicit-euler-time-discretization). The equation can be simplified as:
```{math}
M(\frac{x^{n+1} - \tilde{x} }{\Delta t^2}) = - \nabla C^T(x^{n+1}) \alpha^{-1} C(x^{n+1}) \tag{19}
```
where $n$ represents the iteration count; $\tilde{x} = 2 x^n - x^{n-1} = x^n + v^n \Delta t$ is the **predicted** position. The elastic force $f_{elastic}$ is evaluated by the gradient of the energy potential $U(x)$:
```{math}
f_{elastic} = - \nabla U(x) = - \nabla C^T(x) \alpha^{-1} C(x) \tag{20}
```
XPBD converts it into a constrained problem by introducing a Lagrange multiplier:
```{math}
\lambda_{elastic} = - \tilde{\alpha}^{-1} C(x) \tag{21}
```
where $\lambda_{elastic} = [\lambda_{1}, ..., \lambda_{m}]^{T}$ represents the vector of constraint multipliers; $\tilde{\alpha} = \frac{\alpha}{\Delta{t}^{2}}$ is the compliance matrix.
Substituting $\lambda$ into equation (19, 20), we have the discrete constrained equations of motion and label them as $g$ and $h$, respectively:
```{math}
\begin{equation}
\begin{cases}
g(x, \lambda) = M(x^{n+1} - \tilde{x}) - \nabla C(x^{n+1})^T \lambda^{n+1} = 0 \\
h(x, \lambda) = C(x^{n+1}) + \tilde{\alpha} \lambda^{n+1} = 0
\end{cases}
\tag{22}
\end{equation}
```
Our goal is to find $x$ and $\lambda$ that satisfy:
```{math}
\begin{equation}
\begin{cases}
g(x, \lambda) = 0 \\
h(x, \lambda) = 0
\end{cases}
\tag{23}
\end{equation}
```
Linearizing equations (23), we obtain the following linear Newton subproblem. Notice that we omit the superscript $(n+1)$ that represents the time step; to emphasize per-iteration unknowns we use the subscript $(i)$. For more information on the linearization, please refer to [Appendix: Linearization of g and h](linearization-of-g-and-h). The linearized equation becomes:
```{math}
\begin{equation}
\begin{bmatrix}
\mathbf{K} & -\nabla C^{T}(x_i) \\
\nabla C(x_i) & \tilde{\alpha}
\end{bmatrix}
\begin{bmatrix}
\Delta x \\
\Delta \lambda
\end{bmatrix}
= -
\begin{bmatrix}
g(x_i, \lambda_i) \\
h(x_i, \lambda_i)
\end{bmatrix}
\tag{24}
\end{equation}
```
where $K = \frac{\partial g}{\partial x}$ is related to the second derivative of energy, i.e., it has a *Hessian* term as in equation (12). The system can be solved for $\Delta x$ and the multiplier $\Delta \lambda$ accordingly:
```{math}
\begin{equation}
\begin{cases}
\lambda_{i+1} = \lambda_i + \Delta \lambda \\
x_{i+1} = x_i + \Delta x
\end{cases}
\tag{25}
\end{equation}
```
This iteratively updates $x$ and $\lambda$ that satisfy equation (23). However, as mentioned before, the evaluation of the *Hessian* is too expensive, and the system needs a line search algorithm to be robust. Therefore, XPBD makes two approximations to simplify the problem:
1. $K \approx M$, thus omitting geometric stiffness and Hessian terms. This approximation slows down the local convergence rate but does not change the global error.
2. $g(x_i, \lambda_i)=0$, which is trivially true for the first Newton iteration when initialized with $x_0 = \tilde x$ and $\lambda_0 = 0$.
The simplified problem is given by:
```{math}
\begin{equation}
\begin{bmatrix}
M & -\nabla C^{T}(x_i) \\
\nabla C(x_i) & \tilde{\alpha}
\end{bmatrix}
\begin{bmatrix}
\Delta x \\
\Delta \lambda
\end{bmatrix}
=
-
\begin{bmatrix}
0 \\
h(x_i, \lambda_i)
\end{bmatrix}
\tag{26}
\end{equation}
```
The problem can be written as a reduced system using a partition matrix (see [Appendix: Reduced System via Partition Matrix](reduced-system-via-partition-matrix)) or the Schur complement:
```{math}
[\nabla C(x_i)M^{-1} \nabla C^{T}(x_i) + \tilde{\alpha}] \Delta \lambda = - C(x_i) - \tilde{\alpha} \lambda_i \tag{27}
```
The positional update is given by:
```{math}
\Delta x = M^{-1} \nabla C^{T}(x_i) \Delta \lambda \tag{28}
```
Regarding the update for the Lagrange multiplier, if we take a single constraint with index $j$, $\Delta \lambda_{j}$ is given by:
```{math}
\Delta \lambda_{j} = \frac{ - C_j(x_i) - \tilde{\alpha_j} \lambda_{ij} }{ \nabla C_j(x_i)M^{-1} \nabla C_{j}^{T}(x_i) + \tilde{\alpha_j} } \tag{29}
```
If we set $\alpha_j = 0$, it degenerates to the scaling factor $s$ in classic PBD:
```{math}
s_j = \frac{- C_j(x_i) }{ \nabla C_j(x_i)M^{-1} \nabla C_{j}^{T}(x_i) } \tag{30}
```
and the $\Delta x$ comes to:
```{math}
\Delta x = k_j s_j M^{-1} \nabla C^{T}(x_i) \tag{31}
```
### Step Integration
In the PBD method, the cloth stiffness is related to the time step size and iteration count, which has limited use. We employed the XPBD scheme to update vertex positions in *XRTailor*. The overall algorithm that integrates clothing dynamics is:
