Animation and physical simulation: from keyframes to cloth and fluids

Why this subject matters now

Physical simulation used to be a specialty confined to film effects houses and a handful of engineering codes. A cloth or fluid shot took a render farm overnight, and the person who set it up was a technical director, not a generalist. Three things changed that. The first is that implicit integration and position-based methods made stiff systems, cloth and hair above all, run stably at large time steps, so a garment simulation that once needed thousands of tiny substeps now runs in a game frame. The second is that the GPU turned the embarrassingly parallel inner loops of particle and grid solvers into real-time workloads: a fluid on a grid, a million-particle sand pile, a character's hair, all now update at interactive rates.

The third shift is the one a practitioner today is expected to understand that they were not five years ago: the simulator became differentiable. If every step of a physics engine is a composition of differentiable operations, then the gradient of a final loss, how far a thrown object landed from its target, how much a cloth deviated from a scan, can be propagated back through the entire trajectory to the initial conditions, material parameters, or a control policy. This turns simulation from a forward oracle into a component of a learning system. Reinforcement learning environments like Brax run thousands of parallel rigid-body simulations on an accelerator; differentiable engines like DiffTaichi and Warp let a gradient-based optimizer discover a control signal or a material stiffness directly. Understanding when a gradient through a simulator is meaningful and when a contact discontinuity makes it garbage is now part of the job.

Underneath all of it is a small, stable body of mathematics: interpolation, the algebra of rotations, the numerical analysis of ordinary differential equations, and a few continuum-mechanics energies. The methods that look most modern, XPBD cloth, neural control in a differentiable engine, are recombinations of these pieces. This page derives the pieces, works the arithmetic, and points at the systems that deploy them at scale.

Interpolation: turning keyframes into motion

Animation begins with sparse data. An artist or a motion-capture rig specifies a value, a position, an angle, a blend weight, at a few instants in time called keyframes, and the job of interpolation is to produce a continuous, plausible curve through or near them. The choice of interpolant is the difference between motion that looks mechanical and motion that looks alive, because the eye is exquisitely sensitive to the second derivative of position: a discontinuity in acceleration reads as a jerk.

Linear, then Hermite

Linear interpolation between two keys \(p_0\) at \(t=0\) and \(p_1\) at \(t=1\) is \(\text{lerp}(t) = (1-t)\,p_0 + t\,p_1\). It is continuous in value but its derivative jumps at every key, so a linearly interpolated path has a visible corner wherever the direction changes. The fix is to interpolate not just values but tangents. A cubic Hermite segment is the unique cubic that matches endpoint values \(p_0, p_1\) and endpoint tangents \(m_0, m_1\). Writing the cubic as \(p(t) = a t^3 + b t^2 + c t + d\) and imposing \(p(0)=p_0\), \(p(1)=p_1\), \(p'(0)=m_0\), \(p'(1)=m_1\) gives four equations in the four coefficients. Solving them and regrouping by the data yields the Hermite basis:

Each basis function is a cubic, and one checks directly that \(h_{00}(0)=1\) with every other basis and every basis derivative vanishing at \(t=0\) except \(h_{10}'(0)=1\), which is exactly the statement that \(p_0\) controls the start value and \(m_0\) the start tangent. A whole animation curve is a chain of Hermite segments, and if adjacent segments share endpoint values and tangents the curve is \(C^1\): continuous position and velocity, no corners.

Catmull-Rom: tangents from the data

Hermite splines need tangents, and asking an artist to set a tangent at every key is tedious. The Catmull-Rom spline computes each tangent automatically from the neighboring keys. For a uniform Catmull-Rom through points \(p_{i-1}, p_i, p_{i+1}, p_{i+2}\), the tangent at \(p_i\) is the centered difference \(m_i = \tfrac{1}{2}(p_{i+1} - p_{i-1})\). Substituting these tangents into the Hermite form gives the segment between \(p_i\) and \(p_{i+1}\) as a cubic in the four surrounding points, which in matrix form is

The curve passes exactly through every control point, which is why it is the default for camera paths and motion trajectories: the artist places points and the curve interpolates them with a smooth, tension-free feel. The generalization that lets the artist dial the tension and bias of the tangents is the Kochanek-Bartels spline, standard in animation packages. A caution that matters in practice: Catmull-Rom can overshoot, sending the curve outside the convex hull of the keys between two of them, so for values that must stay bounded (an opacity, a normalized weight) a monotone cubic or a clamped tangent is safer.

Easing: the timing of motion

Interpolation places a value in space; easing places it in time. An easing function is a reparameterization \(u = f(t)\) with \(f(0)=0\), \(f(1)=1\), applied before the spatial interpolant, so that motion accelerates out of a key and decelerates into the next instead of moving at constant speed. The smoothstep polynomial \(f(t) = 3t^2 - 2t^3\) is the cubic Hermite that eases with zero velocity at both ends; its derivative \(6t(1-t)\) vanishes at \(0\) and \(1\), so there is no velocity discontinuity where a still key meets a moving one. Smootherstep, \(6t^5 - 15t^4 + 10t^3\), additionally zeroes the acceleration at the ends and is the quintic used where even the acceleration must be continuous. In interface and game animation the timing is usually authored as a cubic Bezier easing curve with two handle points, the same construction as a Bezier segment, evaluated by solving for the parameter that gives a queried \(x\).

Rotation: why quaternions, and slerp

Interpolating a position is interpolating in a vector space, where the average of two points is a point. Interpolating an orientation is not, because the set of rotations is a curved manifold, the group \(SO(3)\), and the naive average of two rotation matrices is not a rotation. The representation chosen for orientation determines whether interpolation is well behaved, and the near-universal choice in animation and simulation is the unit quaternion.

Gimbal lock: why not Euler angles

Euler angles encode a rotation as three successive turns about coordinate axes, for example yaw then pitch then roll. They are compact and human-readable, but they have a structural defect. When the middle rotation reaches \(\pm 90^\circ\) it aligns the first and third axes, so the first and third rotations now turn about the same line and one degree of freedom is lost: this is gimbal lock. Near that configuration a small change in orientation demands a large, ill-conditioned change in the angles, and interpolating angles across the singularity produces a wild, unphysical tumble. The defect is not a bad choice of axes; it is topological. There is no continuous, singularity-free map from three numbers onto \(SO(3)\), because \(SO(3)\) is not homeomorphic to any open subset of \(\R^3\) in a way that avoids this. Quaternions dodge the problem by using four numbers with one constraint, living on the unit sphere \(S^3\), which double-covers \(SO(3)\) smoothly and without singularities.

The quaternion algebra

A quaternion is \(q = w + x\,i + y\,j + z\,k\), a real part \(w\) and a vector part \((x,y,z)\), with the multiplication rules \(i^2 = j^2 = k^2 = ijk = -1\). These force \(ij = k\), \(jk = i\), \(ki = j\) and the reversed products with a sign flip, so quaternion multiplication is noncommutative. Writing \(q = (w, \mathbf{v})\) with \(\mathbf{v} = (x,y,z)\), the product of \(q_1 = (w_1, \mathbf{v}_1)\) and \(q_2 = (w_2, \mathbf{v}_2)\) collects into

The conjugate is \(\bar q = (w, -\mathbf{v})\), the norm is \(\lVert q\rVert^2 = w^2 + x^2 + y^2 + z^2 = q\bar q\), and for a unit quaternion the inverse is just the conjugate. A rotation by angle \(\theta\) about a unit axis \(\mathbf{n}\) is encoded by the unit quaternion \(q = \big(\cos\tfrac{\theta}{2},\; \sin\tfrac{\theta}{2}\,\mathbf{n}\big)\), and a point \(\mathbf{p}\), embedded as the pure quaternion \((0, \mathbf{p})\), is rotated by the sandwich product \(q\,(0,\mathbf{p})\,\bar q\). The half-angle is the reason \(S^3\) double-covers \(SO(3)\): \(q\) and \(-q\) rotate identically, because rotating the axis-angle by \(2\pi\) negates the quaternion but returns the same orientation. That sign ambiguity is not a nuisance to be ignored; it is exactly the thing to handle before interpolating, as the next section shows.

Deriving slerp

Spherical linear interpolation, slerp, moves along the shortest great-circle arc between two unit quaternions at constant angular speed. Treat \(q_0\) and \(q_1\) as unit vectors in \(\R^4\) separated by angle \(\Omega\), where \(\cos\Omega = q_0 \cdot q_1\). We want a curve \(q(t)\) on the unit sphere, starting at \(q_0\), ending at \(q_1\), that sweeps the arc uniformly, so \(q(t)\) subtends angle \(t\,\Omega\) from \(q_0\). Any point on the plane spanned by \(q_0\) and \(q_1\) is \(q(t) = \alpha(t)\,q_0 + \beta(t)\,q_1\). Two conditions fix \(\alpha\) and \(\beta\): the angle from \(q_0\) to \(q(t)\) is \(t\Omega\), and the angle from \(q(t)\) to \(q_1\) is \((1-t)\Omega\). Taking dot products,

This is a \(2\times 2\) linear system in \(\alpha, \beta\). Solving it, using the identity \(\cos(t\Omega) - \cos((1-t)\Omega)\cos\Omega = \sin((1-t)\Omega)\sin\Omega\) that follows from the angle-subtraction formula, gives the standard result:

The coefficients reduce to \((1-t)\) and \(t\) as \(\Omega \to 0\), so slerp degrades gracefully to lerp for nearby orientations, where one should switch to normalized lerp to avoid dividing by a tiny \(\sin\Omega\). Two practical points. First, before interpolating one flips the sign of \(q_1\) if \(q_0 \cdot q_1 < 0\), choosing the representative on the same hemisphere so the interpolation takes the short way around rather than the long \(> 180^\circ\) arc. Second, slerp between two keys gives constant angular velocity, but chaining slerps across several keys is only \(C^0\) in angular velocity; smooth orientation splines use the squad construction, a spherical Bezier built from slerp, to get \(C^1\) rotation.

Skeletal animation and skinning

A character is animated by a skeleton: a tree of bones, each with a local transform relative to its parent, so that walking the tree from the root composes a world transform for every bone. Animation drives the joint rotations, usually stored as quaternions and interpolated by slerp between keyframes. The skeleton is invisible; what the viewer sees is a skin, a triangle mesh, and the question of skinning is how to move each mesh vertex when the bones move.

Linear blend skinning

Linear blend skinning, also called smooth skinning or the skeletal subspace deformation, binds each vertex to a small set of bones with weights that sum to one. Let vertex \(\mathbf{v}\) in the rest pose be influenced by bones \(j\) with weights \(w_j\), \(\sum_j w_j = 1\). Each bone \(j\) has a rest-pose world transform \(B_j\) and a current animated world transform \(T_j\); the transform that carries a point from the bone's rest frame to its animated frame is \(M_j = T_j B_j^{-1}\). Linear blend skinning takes the weighted average of the transformed vertex over its bones:

The blend is linear in the matrices, which is what makes it fast, one small weighted sum of matrix-vector products per vertex, trivially parallel, and the reason it has run in the vertex shader of essentially every real-time character for two decades. Its defect is also in that linearity. A weighted average of two rigid transforms is not a rigid transform: if a vertex is influenced equally by two bones that have twisted apart, the averaged matrix has a smaller rotational part and the mesh collapses inward. This is the candy-wrapper artifact, the pinch that appears at a twisting elbow or wrist, and it is intrinsic to averaging matrices rather than rotations.

Dual-quaternion skinning

Dual-quaternion skinning fixes the candy-wrapper by blending in a representation where the average of two rigid motions stays close to rigid. A dual quaternion \(\hat q = q_0 + \varepsilon\, q_\varepsilon\), with \(\varepsilon^2 = 0\), packs a rotation into the non-dual part \(q_0\) and the translation into the dual part; unit dual quaternions represent rigid transforms exactly as unit quaternions represent rotations. Skinning then blends the per-bone dual quaternions with the skin weights and normalizes the result, so the rotational component is interpolated on the sphere rather than through it. The twist no longer shrinks the limb, and Kavan and colleagues showed in 2008 that the method costs only modestly more than linear blend skinning and fits the same shader pipeline. Its own artifact is a slight bulge, the joint-bulging that comes from the rotation being preserved too rigidly, but for most characters it is the better default. Production rigs frequently layer corrective blend shapes or a learned corrective on top of either skinning method to capture muscle and skin detail that no bone-weighted blend can.

Numerical integration of the equations of motion

Every physical simulation reduces to the same task: integrate a system of ordinary differential equations forward in time. Newton's second law for a system of particles is second order, \(\mathbf{M}\ddot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \dot{\mathbf{x}})\), and is rewritten as a first-order system in the state \(\mathbf{y} = (\mathbf{x}, \mathbf{v})\):

An integrator advances \(\mathbf{y}\) by a step \(h\). The whole art is in the choice of integrator, because the same physics can be stable or explode depending on it, and the trade is always between accuracy, stability, and cost per step.

The methods

Explicit (forward) Euler evaluates the derivative at the current state: \(\mathbf{y}_{n+1} = \mathbf{y}_n + h\,f(\mathbf{y}_n)\). It is first-order accurate, the local truncation error per step is \(O(h^2)\) and the global error \(O(h)\), and it is the cheapest possible method, one force evaluation per step. It is also the least stable.

Semi-implicit (symplectic) Euler uses the freshly updated velocity to advance position: \(\mathbf{v}_{n+1} = \mathbf{v}_n + h\,\mathbf{M}^{-1}\mathbf{f}(\mathbf{x}_n)\), then \(\mathbf{x}_{n+1} = \mathbf{x}_n + h\,\mathbf{v}_{n+1}\). It costs exactly the same as explicit Euler and is still first order, but it is symplectic for conservative systems: it conserves a slightly perturbed energy exactly, so the energy oscillates in a bounded band instead of drifting. This is why it, not explicit Euler, is the default in games.

Midpoint and RK4. The classical fourth-order Runge-Kutta method takes four derivative evaluations per step, sampling the slope at the start, twice at the middle, and at the end, and combines them so the local error is \(O(h^5)\) and the global error \(O(h^4)\). Halving the step cuts the error by sixteen, which lets RK4 take far larger steps than Euler for the same accuracy, at four times the per-step cost. It is the standard high-accuracy explicit method for smooth, non-stiff systems such as a swinging pendulum or an orbit.

Implicit (backward) Euler evaluates the derivative at the unknown future state: \(\mathbf{y}_{n+1} = \mathbf{y}_n + h\,f(\mathbf{y}_{n+1})\). For a nonlinear \(f\) this is an equation to be solved each step, usually by one or a few Newton iterations, which is expensive. What it buys is unconditional stability on stiff systems: the step size is limited by accuracy, not by blow-up. It is also dissipative, it removes energy, which for a stiff cloth is a feature (the cloth settles) and for a pendulum is a defect (it damps to rest).

Stability regions from the test equation

The stability of a linear multistep or one-step method is analyzed on the scalar test equation \(y' = \lambda y\), whose exact solution is \(y(t) = y_0 e^{\lambda t}\). For \(\mathrm{Re}(\lambda) < 0\) the true solution decays, and a good integrator should produce a numerical sequence that also decays, or at least does not grow. A general oscillatory or damped physical mode looks locally like this equation with \(\lambda\) the eigenvalue of the system's Jacobian, so the scalar analysis governs the real multivariable case mode by mode.

Applying explicit Euler to \(y' = \lambda y\) gives \(y_{n+1} = y_n + h\lambda y_n = (1 + h\lambda)\,y_n\), so \(y_n = (1 + h\lambda)^n y_0\). The numerical solution stays bounded exactly when the amplification factor satisfies \(\lvert 1 + h\lambda\rvert \le 1\). For a real negative \(\lambda\) this is \(-1 \le 1 + h\lambda \le 1\), i.e. \(0 \le h \le 2/\lvert\lambda\rvert\): a hard upper bound on the step. The larger the \(\lvert\lambda\rvert\), the stiffer the mode and the tinier the permitted step, no matter how much the physics itself has settled down. Implicit Euler gives \(y_{n+1} = y_n + h\lambda y_{n+1}\), so \(y_{n+1} = y_n / (1 - h\lambda)\) and the amplification is \(1/\lvert 1 - h\lambda\rvert\). For any \(\mathrm{Re}(\lambda) \le 0\) and any \(h > 0\), \(\lvert 1 - h\lambda\rvert \ge 1\), so the factor is at most one: the method is stable for every step size. This is the whole reason implicit integration exists.

y' = λ y,  λ real < 0        exact:   y_n = e^(λ h n) y_0  → 0

explicit Euler   amp = |1 + hλ|        stable only if  h ≤ 2/|λ|
symplectic Euler (oscillator) bounded if   hω < 2
implicit Euler   amp = 1/|1 - hλ|      stable for ALL h > 0  (A-stable)

The symplectic-Euler bound comes from the same analysis applied to the undamped oscillator \(\ddot x = -\omega^2 x\). Its update matrix maps \((x_n, v_n)\) to \((x_{n+1}, v_{n+1})\) with eigenvalues on the unit circle exactly when the trace condition \(h\omega < 2\) holds, so a symplectic step is stable up to a step of roughly one over the highest frequency, and within that band it neither gains nor loses energy on average.

Mass-spring systems

The simplest deformable model is a network of point masses connected by springs. Each spring between particles \(i\) and \(j\) with rest length \(\ell_0\) and stiffness \(k\) exerts, by Hooke's law, a force along the line between the particles proportional to how far the current length \(\ell = \lVert \mathbf{x}_i - \mathbf{x}_j\rVert\) departs from rest:

A damping term \(-k_d\,(\mathbf{v}_i - \mathbf{v}_j)\) projected onto the spring direction removes oscillation energy. Assembling every spring force on every particle and stepping the state gives a cloth, a rope, or a soft blob, depending on the connectivity. The model is beloved because it is trivial to implement and understand, and distrusted because it is not a faithful discretization of any continuum: the effective material stiffness depends on the mesh resolution and topology, so refining the mesh changes the physics, which the finite element method fixes.

The reason mass-spring cloth was hard before the late 1990s is exactly the stability analysis above. To look like fabric rather than jelly a cloth needs stiff springs, large \(k\), which means a large \(\omega = \sqrt{k/m}\), which forces explicit integrators to a step \(h \lesssim 2/\omega\) that shrinks as the square root of stiffness. A believable garment could need thousands of substeps per frame, and any overshoot detonated the mesh. The field was stuck between cloth that was too soft and cloth that exploded, until implicit integration removed the step-size limit.

The finite element method for deformable solids

Where mass-spring is a heuristic, the finite element method (FEM) is a principled discretization of continuum elasticity, and it is what film-quality flesh, muscle, and soft-body simulation are built on. The object is a solid whose rest shape is meshed into elements, usually tetrahedra. Deformation is described by the map \(\phi\) from a rest point \(\mathbf{X}\) to its deformed position \(\mathbf{x} = \phi(\mathbf{X})\), and the local stretching is captured by the deformation gradient \(\mathbf{F} = \partial\mathbf{x}/\partial\mathbf{X}\), a \(3\times 3\) matrix that is constant within a linear tetrahedron.

Elastic energy and the stiffness matrix

A hyperelastic material stores an energy density \(\Psi(\mathbf{F})\) that penalizes departure from the rest shape, and the total elastic energy is its integral over the body, \(E = \sum_e V_e\,\Psi(\mathbf{F}_e)\), summed over elements of rest volume \(V_e\). The elastic force on a node is the negative gradient of this energy with respect to the node's position, \(\mathbf{f} = -\partial E/\partial\mathbf{x}\), and its derivative, the second derivative of the energy, is the stiffness matrix \(\mathbf{K} = \partial^2 E/\partial\mathbf{x}^2\). For small strains and a linear material the energy is quadratic, \(E = \tfrac12\,\mathbf{u}^\top \mathbf{K}\, \mathbf{u}\) in the displacement \(\mathbf{u}\), and \(\mathbf{K}\) is a constant sparse symmetric matrix assembled from per-element blocks that depend on the material's Lame parameters and the element geometry. Implicit integration of an FEM body reduces, each step, to solving a large sparse linear system built from \(\mathbf{M}\) and \(\mathbf{K}\), which is where most of the compute goes.

Corotational FEM for large rotations

Linear (small-strain) elasticity has a notorious failure: it is not invariant to rotation. Rotating an undeformed element rigidly produces a nonzero linear strain, so a linearly elastic beam that merely swings visibly inflates, the ghost forces grow with the rotation angle. The corotational method fixes this cheaply. Each element's current deformation gradient \(\mathbf{F}\) is factored by a polar decomposition \(\mathbf{F} = \mathbf{R}\,\mathbf{S}\) into a rotation \(\mathbf{R}\) and a symmetric stretch \(\mathbf{S}\); the linear elastic law is then applied in the unrotated frame, so only the genuine stretch \(\mathbf{S} - \mathbf{I}\), not the rotation, generates force. Corotational FEM keeps the speed and simplicity of a linear model while behaving correctly under the large rotations that any animated character undergoes, and it is the workhorse of interactive soft-body simulation. Fully nonlinear hyperelastic models, Neo-Hookean and its variants, go further and remain valid under extreme deformation and inversion, at higher cost per element, and are the standard for offline film-quality flesh.

Rigid-body dynamics

A rigid body cannot deform, so its state is only a position \(\mathbf{x}\) of the center of mass and an orientation, together with their rates, the linear velocity \(\mathbf{v}\) and the angular velocity \(\boldsymbol\omega\). The dynamics split into a translational part governed by \(\mathbf{f} = m\mathbf{a}\) and a rotational part that is considerably subtler.

The inertia tensor and the Newton-Euler equations

Rotational inertia is not a scalar but a \(3\times 3\) symmetric matrix, the inertia tensor \(\mathbf{I}\), which relates angular velocity to angular momentum by \(\mathbf{L} = \mathbf{I}\,\boldsymbol\omega\). It encodes how mass is distributed about the axes; a long rod is easy to spin about its length and hard about a transverse axis, and the tensor's eigenvectors are the principal axes along which spin and momentum align. In the body frame \(\mathbf{I}\) is constant, but in the world frame it rotates with the body: \(\mathbf{I}_\text{world} = \mathbf{R}\,\mathbf{I}_\text{body}\,\mathbf{R}^\top\). The equations of motion are the Newton-Euler equations,

where \(\boldsymbol\tau\) is the applied torque. The cross-product term is the gyroscopic term, present even with no torque, and it is the source of tumbling: a body spun about an intermediate principal axis flips periodically, the tennis-racket effect, purely from \(\boldsymbol\omega \times \mathbf{I}\boldsymbol\omega\). Orientation is integrated as a quaternion, whose rate is \(\dot q = \tfrac12\,(0, \boldsymbol\omega)\,q\); after each step the quaternion is renormalized to stay on the unit sphere, which is cheaper and more robust than integrating a rotation matrix and re-orthonormalizing it.

Collision detection: broad and narrow phase

Bodies collide, and finding contacts is split into two phases for efficiency. The broad phase quickly rejects pairs that cannot possibly touch, using cheap bounding volumes, axis-aligned boxes sorted along a sweep axis, or a spatial hash, to reduce the \(O(n^2)\) pairs to a short list of candidates. The narrow phase then tests each candidate pair exactly, computing whether and where they intersect. For convex shapes the standard narrow-phase algorithm is GJK, the Gilbert-Johnson-Keerthi algorithm, which decides intersection by asking whether the Minkowski difference of the two shapes contains the origin. It never forms that difference explicitly; instead it walks a sequence of simplices, at most a tetrahedron, toward the origin using only a support function that returns the farthest point of a shape in a given direction, so it works for any convex shape defined by such a function and converges in a handful of iterations. When GJK reports penetration, the companion EPA algorithm expands the final simplex to recover the penetration depth and contact normal.

Collision response by impulses

Once a contact with normal \(\mathbf{n}\) is found, the response must instantaneously change the velocities so the bodies separate rather than interpenetrate. The impulse-based method applies an instantaneous impulse \(j\,\mathbf{n}\) at the contact point. Let \(\mathbf{r}_A, \mathbf{r}_B\) be the contact point relative to each body's center of mass, and let the relative velocity at the contact along the normal be \(v_\text{rel} = \mathbf{n}\cdot(\mathbf{v}_A + \boldsymbol\omega_A\times\mathbf{r}_A - \mathbf{v}_B - \boldsymbol\omega_B\times\mathbf{r}_B)\). A coefficient of restitution \(e\in[0,1]\) prescribes that the outgoing normal speed be \(-e\,v_\text{rel}\), perfectly elastic at \(e=1\) and perfectly inelastic at \(e=0\). Solving for the impulse that achieves this, accounting for how the impulse changes both the linear and angular velocities of each body, gives

The denominator is the effective inverse mass at the contact along \(\mathbf{n}\): it blends the two bodies' translational inverse masses with their rotational responses to an impulse at that lever arm. The impulse then updates \(\mathbf{v}_A \mathrel{+}= j\,\mathbf{n}/m_A\) and \(\boldsymbol\omega_A \mathrel{+}= \mathbf{I}_A^{-1}(\mathbf{r}_A\times j\,\mathbf{n})\), with opposite signs for \(B\). Friction is a second impulse in the tangent plane, capped by the Coulomb limit \(\lvert j_t\rvert \le \mu\,j\).

Constraint-based contact and the LCP

Sequential impulses handle one contact at a time and iterate, which is fast and robust and is what most game engines ship. A more principled view treats all contacts at once as a complementarity problem. At each contact the normal force \(\lambda \ge 0\) can only push, never pull, and the normal separation velocity \(a \ge 0\) must be non-negative, and the two are complementary: either the bodies are separating with zero force, or they are in contact with positive force, so \(\lambda\,a = 0\). Collecting these across all contacts, with \(a = \mathbf{A}\lambda + \mathbf{b}\) the linear relation between forces and resulting separation velocities, yields the linear complementarity problem \(a = \mathbf{A}\lambda + \mathbf{b},\; a\ge 0,\; \lambda\ge 0,\; \lambda^\top a = 0\). Solving the LCP gives a set of contact impulses that are simultaneously consistent, which matters for stacks and joints where treating contacts independently drifts. Modern engines solve it approximately with a projected Gauss-Seidel iteration, which is exactly the sequential-impulse loop seen from the optimization side.

Cloth and hair: implicit integration made them practical

Cloth is a stiff mass-spring or thin-shell system, and the stability argument above is the reason it resisted real-time simulation for so long. The turning point was the 1998 result of Baraff and Witkin, who argued that the right answer to stiff cloth is not smaller steps but implicit integration, and showed how to take large steps stably. Their scheme is implicit (backward) Euler on the cloth's mass-spring-bending system. Writing the update in terms of the velocity change \(\Delta\mathbf{v}\) and linearizing the forces about the current state, the implicit step becomes a single large sparse symmetric linear system

with \(\partial\mathbf{f}/\partial\mathbf{x}\) the stiffness matrix. The matrix on the left is symmetric positive definite, so the system is solved with a modified conjugate-gradient iteration, and the modification is how Baraff and Witkin enforce constraints such as a pinned vertex directly inside the solver by filtering the search directions. The payoff is that the step size is no longer governed by the spring stiffness; a garment that would have needed thousands of explicit substeps runs at a handful of implicit steps per frame. The same implicit-solve template, mass matrix minus scaled stiffness, conjugate gradient, underlies most production cloth to this day, augmented with strain limiting to stop over-stretch, bending energies for a realistic drape, and robust collision handling to keep the cloth from passing through itself and the body. Hair is the same physics on many thin strands, typically modeled as inextensible rods with bending and twisting energies, and it shares the stiffness problem and the implicit or position-based cure.

Fluids: from Navier-Stokes to stable fluids

Fluid motion is governed by the incompressible Navier-Stokes equations, a statement of momentum conservation together with the constraint that the fluid neither compresses nor expands:

The velocity field \(\mathbf{u}\) is pushed around by four effects: advection \((\mathbf{u}\cdot\nabla)\mathbf{u}\), the fluid carrying its own momentum; the pressure gradient \(\nabla p\), which enforces incompressibility; viscous diffusion \(\nu\nabla^2\mathbf{u}\); and body forces \(\mathbf{g}\) like gravity. The advection term is nonlinear and is what makes the equations hard, and a naive explicit discretization of it is unstable at any interesting resolution, which kept fluids out of interactive graphics until 1999.

Stam's stable fluids

Stam's 1999 stable fluids method made real-time fluid possible by splitting the step into substeps that are each unconditionally stable, at the cost of accuracy. The key move is the semi-Lagrangian advection. Rather than pushing quantities forward, which can overshoot and blow up, it asks, for each grid cell, where the fluid now in this cell came from a step ago: trace the velocity field backward from the cell center by \(-h\,\mathbf{u}\), and interpolate the old field at that departure point. Because the new value is an interpolated average of old values, it can never exceed the range of the old field, so semi-Lagrangian advection is unconditionally stable regardless of step size, the fluid analog of the same averaging argument that stabilizes semi-Lagrangian cloth. The trade is numerical diffusion: the interpolation smears the field, so vorticity is lost and the fluid looks too viscous, which later methods counter with vorticity confinement or higher-order advection.

After advection and force integration the velocity field is generally not divergence-free, so the final substep is pressure projection, which subtracts off the compressible part. By the Helmholtz decomposition any vector field splits uniquely into a divergence-free part and the gradient of a scalar. Setting \(\mathbf{u} = \mathbf{u}^* - \nabla p\) and imposing \(\nabla\cdot\mathbf{u} = 0\) gives a Poisson equation for the pressure, \(\nabla^2 p = \nabla\cdot\mathbf{u}^*\), a large sparse symmetric linear system solved each step; subtracting its gradient leaves a divergence-free field. Projection is where the incompressibility constraint is actually enforced, and the Poisson solve is the most expensive part of a grid fluid, which is why multigrid and preconditioned conjugate gradient are the standard solvers.

Grid, particle, and hybrid

There are two families of fluid discretization and a hybrid that dominates production. The Eulerian (grid) view, above, stores velocity and pressure on a fixed grid; it enforces incompressibility crisply through the pressure Poisson solve but suffers numerical diffusion in advection. The Lagrangian view, smoothed-particle hydrodynamics (SPH), carries the fluid as particles that move with the flow and interact through smoothing kernels; it handles free surfaces and splashes naturally and conserves mass exactly, but enforcing incompressibility on scattered particles is awkward and it can look clumpy. The hybrid FLIP and PIC methods, and the affine variant APIC, get the best of both: they carry quantities on particles for advection, so there is no grid diffusion, but transfer to a grid each step to do the pressure projection, so incompressibility is enforced well. FLIP or APIC on a grid is the standard for high-end liquid effects, and the material point method (MPM) extends the same particle-grid transfer to snow, sand, and elastoplastic solids.

Position-based dynamics

Position-based dynamics (PBD), introduced by Muller and colleagues in 2007, is a different answer to stability than implicit integration, and it has taken over real-time simulation because it is simple, fast, and unconditionally stable by construction. Instead of computing forces, integrating to velocities, and integrating to positions, PBD works directly on positions. Each step first does an unconstrained prediction, moving every particle by gravity and its velocity, and then repeatedly projects the predicted positions onto a set of constraints until they are satisfied.

A constraint is a function \(C(\mathbf{x})\) that should equal zero, a distance constraint \(C = \lVert\mathbf{x}_i - \mathbf{x}_j\rVert - \ell_0\) for a spring, or a volume or bending constraint. PBD linearizes \(C\) and moves the involved particles along the gradient of \(C\), weighted by inverse mass, by exactly the amount that zeroes the linearized constraint:

Iterating this projection over all constraints, in Gauss-Seidel fashion, drives the system toward a configuration that satisfies them, and the velocity is recovered afterward as the change in position over the step. The method is stable because it never applies a force that could overshoot; it moves positions directly and a projection cannot make a constraint more violated. Its defect is that the effective stiffness depends on the number of solver iterations and the time step, so the material gets stiffer as one iterates more, which is physically arbitrary.

XPBD

Extended position-based dynamics (XPBD), from Macklin and colleagues in 2016, removes that defect by giving each constraint a physical compliance \(\alpha\), the inverse of a stiffness, and carrying a Lagrange multiplier for the constraint across the iteration. The per-iteration position update becomes

with the position moved by \(\Delta\lambda\,\nabla C\) weighted by inverse mass. Because \(\tilde\alpha\) carries the true compliance divided by \(h^2\), the converged stiffness is now independent of the iteration count and time step, so a cloth authored at a given stiffness behaves the same regardless of solver budget. XPBD is the backbone of modern real-time cloth, soft bodies, and granular effects, and it is what runs when a game character's cape simulates in a millisecond.

Differentiable simulation

A simulator is a function that maps initial conditions, parameters, and controls to a trajectory, and if every operation in it is differentiable then automatic differentiation gives the gradient of any scalar loss on the trajectory with respect to those inputs. This reframes simulation as a layer in a learning system. The gradient of "the thrown ball missed the target by this much" with respect to the release velocity is available directly, so gradient descent finds the throw, and the gradient with respect to a neural control policy's weights trains the policy through the physics rather than around it. DiffTaichi, from Hu and colleagues in 2019 and 2020, built a differentiable programming language for physics in which a handful of lines specify a mass-spring, MPM, or fluid simulator whose gradient is generated automatically, and demonstrated learning locomotion controllers and material parameters end to end. Brax, from Freeman and colleagues at Google in 2021, put a rigid-body physics engine entirely in JAX so that thousands of environments run in parallel on an accelerator and the whole simulation is differentiable and vmap-able, collapsing the reinforcement-learning wall-clock for locomotion tasks. NVIDIA's Warp brings the same differentiable, GPU-kernel approach to a Python simulation framework.

The caveat is essential and is where careful practitioners earn their keep: contact and collision are discontinuous, and the gradient through a hard contact is either zero or undefined. A ball that bounces has a trajectory whose dependence on initial height is piecewise smooth with a kink at each bounce, so a naive gradient can point the wrong way or vanish. Differentiable engines address this with softened or randomized contact models, with gradient estimators that blend the analytic gradient with a stochastic one, and by choosing losses that avoid stepping across the discontinuity. Knowing when the gradient of a simulator is trustworthy, smooth elastic deformation, yes; a stiff impulsive collision, treat with suspicion, is a core competence in this area.

Implementation

The integrators first, because they are the foundation everything else rests on. The block below implements explicit Euler, symplectic Euler, and RK4, then reproduces the pendulum energy-drift numbers from Problem 4. It is runnable as written and prints the drifts quoted in the table.

import numpy as np

g, L = 9.81, 1.0                       # pendulum: gravity, length

def deriv(s):                          # state s = [theta, omega]
    theta, omega = s
    return np.array([omega, -g / L * np.sin(theta)])

def energy(s):                         # specific energy, conserved exactly
    theta, omega = s
    return 0.5 * omega**2 + g / L * (1.0 - np.cos(theta))

def step_euler(s, h):                  # explicit forward Euler, O(h) global
    return s + h * deriv(s)

def step_symplectic(s, h):             # semi-implicit: new v drives new x
    theta, omega = s
    omega = omega + h * (-g / L * np.sin(theta))
    theta = theta + h * omega
    return np.array([theta, omega])

def step_rk4(s, h):                    # classical 4th order, 4 evals
    k1 = deriv(s)
    k2 = deriv(s + 0.5 * h * k1)
    k3 = deriv(s + 0.5 * h * k2)
    k4 = deriv(s + h * k3)
    return s + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)

def run(step, h=0.01, T=10.0):
    s = np.array([1.0, 0.0])           # released from 1 rad at rest
    E0 = energy(s)
    for _ in range(int(T / h)):
        s = step(s, h)
    return (energy(s) - E0) / E0 * 100.0

for name, step in [("euler", step_euler),
                   ("symplectic", step_symplectic),
                   ("rk4", step_rk4)]:
    print(f"{name:11s} energy drift over 10 s: {run(step):+.4f} %")
# euler       energy drift over 10 s: +113.3634 %
# symplectic  energy drift over 10 s:  -1.2489 %
# rk4         energy drift over 10 s:  -0.0000 %

import jax, jax.numpy as jnp
from jax import lax

g, L = 9.81, 1.0

def deriv(s):
    theta, omega = s
    return jnp.array([omega, -g / L * jnp.sin(theta)])

def energy(s):
    theta, omega = s
    return 0.5 * omega**2 + g / L * (1.0 - jnp.cos(theta))

def step_rk4(s, h):
    k1 = deriv(s)
    k2 = deriv(s + 0.5 * h * k1)
    k3 = deriv(s + 0.5 * h * k2)
    k4 = deriv(s + h * k3)
    return s + h / 6.0 * (k1 + 2 * k2 + 2 * k3 + k4)

@jax.jit
def run(h=0.01, T=10.0):
    s0 = jnp.array([1.0, 0.0])
    n = int(T / h)
    # lax.scan keeps the whole rollout on the accelerator
    s_final, _ = lax.scan(lambda s, _: (step_rk4(s, h), None), s0, None, length=n)
    return (energy(s_final) - energy(s0)) / energy(s0) * 100.0

print("rk4 energy drift over 10 s:", float(run()), "%")

# Because run is a pure JAX function, jax.grad(run) differentiates the entire
# 1000-step rollout: the sensitivity of the final energy drift to the step h,
# release angle, or gravity is one autodiff call away. That is the seed of a
# differentiable simulator.

#include <array>
#include <cmath>
#include <cstdio>

using State = std::array<double, 2>;   // {theta, omega}
constexpr double g = 9.81, L = 1.0;

State deriv(const State& s) {
    return { s[1], -g / L * std::sin(s[0]) };
}
double energy(const State& s) {
    return 0.5 * s[1] * s[1] + g / L * (1.0 - std::cos(s[0]));
}
State axpy(const State& a, double h, const State& b) {   // a + h*b
    return { a[0] + h * b[0], a[1] + h * b[1] };
}
State step_rk4(const State& s, double h) {
    State k1 = deriv(s);
    State k2 = deriv(axpy(s, 0.5 * h, k1));
    State k3 = deriv(axpy(s, 0.5 * h, k2));
    State k4 = deriv(axpy(s, h, k3));
    return { s[0] + h / 6.0 * (k1[0] + 2*k2[0] + 2*k3[0] + k4[0]),
             s[1] + h / 6.0 * (k1[1] + 2*k2[1] + 2*k3[1] + k4[1]) };
}
int main() {
    State s = {1.0, 0.0};
    double E0 = energy(s), h = 0.01;
    for (int i = 0; i < 1000; ++i) s = step_rk4(s, h);
    std::printf("rk4 drift over 10 s: %+.6f %%\n", (energy(s) - E0) / E0 * 100.0);
}

Next, one substep of a position-based-dynamics cloth: predict positions under gravity, then project distance constraints with inverse-mass weighting, exactly the update derived in Problem 6. Pinned particles are given zero inverse mass so the projection leaves them fixed. This is the shape of the inner loop that runs in a real-time cloth engine.

import numpy as np

def pbd_step(x, v, w, edges, rest, h=1/60, iters=20, g=np.array([0.0, -9.81, 0.0])):
    """One PBD substep for a mass-spring cloth.
    x    (N,3) positions      v (N,3) velocities
    w    (N,)  inverse masses (0 == pinned)
    edges (E,2) index pairs   rest (E,) rest lengths
    """
    p = x + h * v + (h * h) * g * (w[:, None] > 0)   # predict (skip pinned gravity)
    for _ in range(iters):                            # Gauss-Seidel projection
        for e, (i, j) in enumerate(edges):
            d = p[i] - p[j]
            L = np.linalg.norm(d)
            if L < 1e-9:
                continue
            C = L - rest[e]                           # constraint value
            n = d / L                                 # gradient direction
            denom = w[i] + w[j]
            if denom == 0:
                continue
            dp = (C / denom) * n                      # scaled correction
            p[i] -= w[i] * dp
            p[j] += w[j] * dp
    v = (p - x) / h                                   # velocity from motion
    return p, v

# 3-particle strip pinned at one end, released under gravity
x = np.array([[0.,0.,0.], [1.,0.,0.], [2.,0.,0.]])
v = np.zeros((3, 3))
w = np.array([0.0, 1.0, 1.0])                          # particle 0 pinned
edges = np.array([[0, 1], [1, 2]])
rest = np.array([1.0, 1.0])
for _ in range(120):                                  # 2 seconds at 60 Hz
    x, v = pbd_step(x, v, w, edges, rest)
print("settled positions:\n", np.round(x, 3))         # hangs down from the pin

#include <vector>
#include <array>
#include <cmath>

using V3 = std::array<double, 3>;
struct Edge { int i, j; double rest; };

static double norm(const V3& a) {
    return std::sqrt(a[0]*a[0] + a[1]*a[1] + a[2]*a[2]);
}

// One PBD substep. w[k] == 0 pins particle k.
void pbd_step(std::vector<V3>& x, std::vector<V3>& v,
              const std::vector<double>& w,
              const std::vector<Edge>& edges,
              double h, int iters) {
    const V3 g = {0.0, -9.81, 0.0};
    std::vector<V3> p = x;
    for (size_t k = 0; k < x.size(); ++k)
        for (int d = 0; d < 3; ++d)
            p[k][d] += h * v[k][d] + (w[k] > 0 ? h*h*g[d] : 0.0);

    for (int it = 0; it < iters; ++it)
        for (const auto& e : edges) {
            V3 dv = { p[e.i][0]-p[e.j][0], p[e.i][1]-p[e.j][1], p[e.i][2]-p[e.j][2] };
            double Lc = norm(dv);
            if (Lc < 1e-9) continue;
            double C = Lc - e.rest, denom = w[e.i] + w[e.j];
            if (denom == 0.0) continue;
            double s = C / (denom * Lc);
            for (int d = 0; d < 3; ++d) {
                p[e.i][d] -= w[e.i] * s * dv[d];
                p[e.j][d] += w[e.j] * s * dv[d];
            }
        }
    for (size_t k = 0; k < x.size(); ++k)
        for (int d = 0; d < 3; ++d) {
            v[k][d] = (p[k][d] - x[k][d]) / h;
            x[k][d] = p[k][d];
        }
}

Finally, the quaternion slerp of Problem 1, with the sign-flip that keeps the interpolation on the short arc and the near-parallel fallback to normalized lerp.

import numpy as np

def slerp(q0, q1, t):
    q0 = q0 / np.linalg.norm(q0)
    q1 = q1 / np.linalg.norm(q1)
    dot = float(q0 @ q1)
    if dot < 0.0:                       # pick the short way around SO(3)
        q1, dot = -q1, -dot
    if dot > 0.9995:                    # nearly parallel: lerp then renormalize
        q = q0 + t * (q1 - q0)
        return q / np.linalg.norm(q)
    omega = np.arccos(dot)              # angle between the quaternions
    so = np.sin(omega)
    a = np.sin((1 - t) * omega) / so
    b = np.sin(t * omega) / so
    return a * q0 + b * q1

c = np.cos(np.pi / 4)
q0 = np.array([c, c, 0, 0])            # 90 deg about x
q1 = np.array([c, 0, c, 0])            # 90 deg about y
print(np.round(slerp(q0, q1, 0.5), 6))
# [0.816497 0.408248 0.408248 0.      ]  -> 70.53 deg about (1,1,0)/sqrt(2)

#include <array>
#include <cmath>

using Quat = std::array<double, 4>;    // {w, x, y, z}

static double dot(const Quat& a, const Quat& b) {
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + a[3]*b[3];
}

Quat slerp(Quat q0, Quat q1, double t) {
    double d = dot(q0, q1);
    if (d < 0.0) { for (auto& c : q1) c = -c; d = -d; }   // short arc
    if (d > 0.9995) {                                       // lerp fallback
        Quat q; double n = 0.0;
        for (int i = 0; i < 4; ++i) { q[i] = q0[i] + t*(q1[i]-q0[i]); n += q[i]*q[i]; }
        n = std::sqrt(n);
        for (auto& c : q) c /= n;
        return q;
    }
    double omega = std::acos(d), so = std::sin(omega);
    double a = std::sin((1.0 - t) * omega) / so;
    double b = std::sin(t * omega) / so;
    Quat q;
    for (int i = 0; i < 4; ++i) q[i] = a * q0[i] + b * q1[i];
    return q;
}

How it is done in practice

The gap between these derivations and a shipped system is mostly robustness and scale. A game physics engine such as Bullet, PhysX, or Jolt runs a fixed small time step, typically \(1/60\) or \(1/120\) of a second with substepping, and spends its budget on the broad phase (a sweep-and-prune or a dynamic bounding-volume tree), the narrow phase (GJK and EPA for convex shapes, specialized tests for boxes and spheres), and a constraint solver that is a projected Gauss-Seidel sweep over contacts and joints, warm-started from the previous frame's impulses so that a stack of boxes converges in a few iterations rather than hundreds. Numerical robustness dominates the engineering: contact manifolds must be reduced to a stable set of points, penetration must be resolved without adding energy (Baumgarte stabilization or a split-impulse scheme), and rounding must never let a resting object jitter or sink.

Real-time cloth and soft bodies in games are almost universally XPBD now, running dozens of solver iterations over tens of thousands of constraints per frame on the GPU, with strain limiting and self-collision handled by a spatial hash. Film production runs the heavier machinery: implicit FEM flesh with Neo-Hookean energies solved by Newton with a sparse Cholesky or conjugate-gradient inner solve, FLIP or APIC liquids on adaptive grids with millions of particles, and the material point method for snow and sand, all offline at minutes to hours per frame. The differentiable and massively-parallel engines are the newest production category: Brax and Warp and Isaac Gym run thousands of rigid-body environments at once on a single accelerator to train robot controllers, where the value is not photorealism but throughput, tens of thousands of simulated steps per second per environment, so a policy sees billions of transitions in hours.

The current research frontier

Differentiable simulation is the most active front. DiffTaichi (Hu et al., MIT and others) and Warp (NVIDIA) make gradients through elastic, fluid, and MPM simulators routine, and the open question is contact: how to define a useful gradient through the discontinuity of an impact. Groups at MIT, Berkeley, Toronto, and Google DeepMind have proposed smoothed contact, randomized-smoothing gradient estimators, and implicit-function-theorem gradients through the LCP, and there is no settled answer. A parallel line pushes learned surrogates: graph-network simulators from DeepMind (Sanchez-Gonzalez, Pfaff, and colleagues) learn to predict particle or mesh dynamics directly, trading physical guarantees for speed and for differentiability by construction, and are competitive on fluids and deformables where a fast approximate answer beats a slow exact one.

On the classical side, the projective-dynamics and its generalizations (Bouaziz and colleagues at EPFL, later work at IST Austria and elsewhere) reformulate implicit integration as an optimization with a fixed system matrix that can be prefactored, giving fast, stable elastic simulation, and the family of solvers around it, ADMM and Anderson acceleration, is where much of the deformable-body speed comes from. In fluids, the push is toward better advection with less diffusion (BFECC, the reflection method, and monotonic higher-order schemes) and toward the material point method as a unifying framework across snow, sand, foam, and cloth, with contributions from UCLA, the University of Pennsylvania, UBC, and others. Robotics has pulled rigid-body simulation toward differentiable, GPU-parallel, contact-rich engines, MuJoCo (now at Google DeepMind), Brax, Isaac, and analytically-differentiable variants, because sim-to-real transfer of learned control depends on both throughput and gradient quality. The unifying theme across all of it is that the boundary between simulation and learning has dissolved: the simulator is now something one optimizes through, not merely runs.

Open source to read

• taichi-dev/taichi: a differentiable programming language for high-performance physics; read the MPM and mass-spring examples in python/taichi/examples to see a full simulator and its autodiff in a page of code.
• google/brax: a JAX rigid-body engine for massively parallel, differentiable RL; open brax/positional and brax/generalized to compare the position-based and reduced-coordinate dynamics backends.
• NVIDIA/warp: a Python framework for differentiable GPU kernels; the warp/sim module has cloth, FEM, and rigid-body integrators with gradients, a good read for how autodiff and CUDA kernels coexist.
• bulletphysics/bullet3: the production game engine; start in src/BulletDynamics/ConstraintSolver for the sequential-impulse solver and in src/BulletCollision for GJK and the broad phase.
• InteractiveComputerGraphics/PositionBasedDynamics: the reference PBD and XPBD implementation from the group that invented the method; the constraint-projection files are the clearest statement of the update in this page.
• google-deepmind/mujoco: the contact-rich robotics simulator; the engine core in src/engine shows the soft-constraint contact model and the reduced-coordinate formulation used for control.
• NVIDIAGameWorks/FleX: the unified particle solver behind cloth, fluids, and rigids in games, a case study in one position-based framework spanning every material.

Common misconceptions

"Smaller time steps make explicit Euler stable on a spring." On an undamped oscillator explicit Euler's energy grows by \(\sqrt{1 + (h\omega)^2} > 1\) every step, at any step size; shrinking \(h\) only slows the blow-up, it never stops it. Symplectic Euler, at the same cost, is bounded up to \(h\omega < 2\). The instability is qualitative, not a matter of resolution.

"Implicit integration is more accurate." It is more stable, which is a different property. Implicit Euler is only first-order accurate and is strongly dissipative; it lets a stiff cloth take large steps without exploding, but it damps energy and rounds off high-frequency motion. Accuracy and stability are separate axes, and the reason to go implicit is stability on stiff systems, not fidelity.

"Quaternions are just a compact rotation format." The compactness is incidental; the reason they dominate is that they interpolate on \(SO(3)\) without singularities and compose without drift, where Euler angles gimbal-lock and rotation matrices lose orthonormality under repeated multiplication. A quaternion is the double cover of the rotation group, and slerp is a geodesic on it.

"Linear blend skinning is wrong and should be replaced everywhere." It is the right tool for most vertices most of the time; the candy-wrapper artifact only bites at large twists near a joint. Dual-quaternion skinning fixes the twist but introduces a bulge, so production rigs use linear blend skinning by default and switch or add correctives only where the twist is severe.

"Mass-spring cloth is a discretization of a real material." It is a heuristic whose effective stiffness and anisotropy depend on the mesh topology, so refining the mesh changes the physics. The finite element method is the convergent discretization; mass-spring is chosen for simplicity, not fidelity.

"You can always backpropagate through a physics engine." Smooth elastic and fluid dynamics differentiate cleanly, but the gradient through a hard contact or a collision is zero or undefined, and a naive autodiff there yields a useless or misleading gradient. Differentiable engines need softened or randomized contact models precisely because the exact gradient across an impact carries no information.

Self-check

References

1. R. Parent. Computer Animation: Algorithms and Techniques, 3rd ed. Morgan Kaufmann, 2012.
2. R. Bridson. Fluid Simulation for Computer Graphics, 2nd ed. A K Peters/CRC Press, 2015.
3. D. Baraff and A. Witkin. Physically Based Modeling: Principles and Practice. SIGGRAPH course notes, 1997 (updated 2001). Pixar / Carnegie Mellon. graphics.pixar.com/pbm2001
4. D. Baraff and A. Witkin. Large Steps in Cloth Simulation. SIGGRAPH, 1998. doi:10.1145/280814.280821
5. J. Stam. Stable Fluids. SIGGRAPH, 1999. Alias|wavefront. doi:10.1145/311535.311548
6. K. Shoemake. Animating Rotation with Quaternion Curves. SIGGRAPH, 1985. doi:10.1145/325334.325242
7. M. Muller, B. Heidelberger, M. Hennix, and J. Ratcliff. Position Based Dynamics. Journal of Visual Communication and Image Representation, 18(2), 2007. doi:10.1016/j.jvcir.2007.01.005
8. M. Macklin, M. Muller, and N. Chentanez. XPBD: Position-Based Simulation of Compliant Constrained Dynamics. Motion in Games, 2016. doi:10.1145/2994258.2994272
9. L. Kavan, S. Collins, J. Zara, and C. O'Sullivan. Geometric Skinning with Approximate Dual Quaternion Blending. ACM Transactions on Graphics, 27(4), 2008. doi:10.1145/1409625.1409627
10. E. Sifakis and J. Barbic. FEM Simulation of 3D Deformable Solids: A Practitioner's Guide to Theory, Discretization and Model Reduction. SIGGRAPH course notes, 2012. viterbi-web.usc.edu/~jbarbic/femdefo
11. E. G. Gilbert, D. W. Johnson, and S. S. Keerthi. A Fast Procedure for Computing the Distance Between Complex Objects in Three-Dimensional Space. IEEE Journal of Robotics and Automation, 4(2), 1988. doi:10.1109/56.2083
12. Y. Hu, L. Anderson, T.-M. Li, Q. Sun, N. Carr, J. Ragan-Kelley, and F. Durand. DiffTaichi: Differentiable Programming for Physical Simulation. ICLR, 2020. arXiv:1910.00935
13. Y. Hu, T.-M. Li, L. Anderson, J. Ragan-Kelley, and F. Durand. Taichi: A Language for High-Performance Computation on Spatially Sparse Data Structures. ACM Transactions on Graphics, 38(6), 2019. doi:10.1145/3355089.3356506
14. C. D. Freeman, E. Frey, A. Raichuk, S. Girgin, I. Mordatch, and O. Bachem. Brax: A Differentiable Physics Engine for Large Scale Rigid Body Simulation. 2021. arXiv:2106.13281
15. J. Bender, M. Muller, and M. Macklin. A Survey on Position Based Dynamics. Eurographics Tutorials, 2017. doi:10.2312/egt.20171034
16. M. Macklin, M. Muller, N. Chentanez, and T.-Y. Kim. Unified Particle Physics for Real-Time Applications. ACM Transactions on Graphics, 33(4), 2014. doi:10.1145/2601097.2601152
17. M. Muller, B. Heidelberger, M. Teschner, and M. Gross. Meshless Deformations Based on Shape Matching. SIGGRAPH, 2005. doi:10.1145/1073204.1073216
18. Y. Zhu and R. Bridson. Animating Sand as a Fluid. SIGGRAPH, 2005. doi:10.1145/1073204.1073298
19. C. Jiang, C. Schroeder, A. Selle, J. Teran, and A. Stomakhin. The Affine Particle-In-Cell Method. ACM Transactions on Graphics, 34(4), 2015. doi:10.1145/2766996
20. S. Bouaziz, S. Martin, T. Liu, L. Kavan, and M. Pauly. Projective Dynamics: Fusing Constraint Projections for Fast Simulation. ACM Transactions on Graphics, 33(4), 2014. doi:10.1145/2601097.2601116
21. A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia. Learning to Simulate Complex Physics with Graph Networks. ICML, 2020. arXiv:2002.09405
22. E. Todorov, T. Erez, and Y. Tassa. MuJoCo: A Physics Engine for Model-Based Control. IEEE/RSJ IROS, 2012. doi:10.1109/IROS.2012.6386109