Constraint satisfaction in physics: the method of Lagrange multipliers

1. Why not just use springs?

Найпростіший спосіб «прив'язати» дві точки — з'єднати їх дуже жорсткою пружиною: F = −k·(|r| − L)·r̂. Проблема в тому, що жорсткість k, потрібна для «нерозтяжної» мотузки, робить систему чисельно жорсткою (stiff) — явний інтегратор вимагає надзвичайно малого dt, інакше симуляція вибухає осциляціями.

Метод множників Лагранжа дозволяє знайти точну силу, яка утримує обмеження — без будь-якої жорсткості чи власної частоти коливань. Це та сама математика, яка лежить в основі узагальненої механіки Лагранжа: мінімізація дії за умови обмежень.

2. Constrained Lagrangian mechanics

Обмеження записують як неявну функцію C(q) = 0, де q — узагальнені координати системи. Наприклад, для жорсткого стрижня довжини L між точками p₁ і p₂:

Рівняння руху з обмеженнями (constrained Euler-Lagrange) додають до звичайних сил доданок λ · ∇C — силу, спрямовану вздовж градієнта обмеження, з невідомим масштабом λ (множник Лагранжа):

Множник λ обчислюється так, щоб обмеження точно виконувалось на кожному кроці — не наближено, як з пружиною, а точно (з точністю до похибки чисельного інтегрування).

3. What the multiplier λ actually is

λ — це величина сили обмеження, необхідної, щоб утримати систему на многовиді C(q) = 0. Фізично це натяг мотузки, реакція опори, сила зчеплення в шарнірі. Знак λ несе фізичний зміст: для мотузки λ ≥ 0 (мотузка тягне, але не штовхає) — це однобічне обмеження (unilateral constraint), на відміну від жорсткого стрижня, де λ може бути будь-якого знаку (двобічне, bilateral).

4. Velocity-level constraint form

У фізичних рушіях реального часу зручніше працювати не з позиційним обмеженням C(q) = 0, а з його похідною за часом — обмеженням на швидкість:

Підставляючи рівняння руху, отримуємо лінійну систему відносно λ:

Матрицю J · M⁻¹ · Jᵀ називають ефективною масою (effective mass) обмеження — вона показує, наскільки «легко» система реагує на одиничний імпульс уздовж напрямку обмеження.

5. Baumgarte stabilization

Робота лише з Ċ = 0 утримує обмеження на рівні швидкості, але не виправляє накопичену позиційну похибку (drift) через похибки чисельного інтегрування — мотузка поступово «розтягується» на десяті частки відсотка за кроками. Baumgarte-стабілізація додає корекційний доданок, пропорційний позиційній похибці:

Параметр β ∈ [0.1, 0.3] контролює швидкість корекції дрейфу: занадто велике β «висмикує» систему занадто різко (додає штучну енергію), занадто мале — дрейф накопичується видимо. Це компроміс між точністю і стабільністю — типова інженерна проблема числових методів.

6. Sequential impulse solvers

Коли обмежень багато (десятки контактів, шарнірів у ragdoll), розв'язання повної системи J·M⁻¹·Jᵀ·λ = b напряму — дорого. Практичні рушії (Box2D, Bullet, PhysX) використовують Sequential Impulse (Projected Gauss-Seidel): розв'язують кожне обмеження окремо, ітеративно, приймаючи інші за фіксовані, поки система не збіжиться.

Зазвичай достатньо 4–10 ітерацій для правдоподібної збіжності в реальному часі — точний розв'язок не потрібен, важлива стабільність кадр за кадром.

7. Constraint methods compared

8. Псевдокод

Distance constraint with Baumgarte (JS)

// Розв'язати λ для обмеження відстані з Baumgarte-корекцією
function solveDistanceConstraint(a, b, restLength, dt, beta = 0.2) {
  const delta = sub(b.pos, a.pos);
  const dist = length(delta);
  const n = scale(delta, 1 / dist);        // unit direction
  const C = dist − restLength;               // position error

  // Relative velocity along constraint direction
  const relVel = dot(sub(b.vel, a.vel), n);

  // Effective mass along n
  const invMassSum = a.invMass + b.invMass;
  const bias = (beta / dt) * C;             // Baumgarte term

  const lambda = −(relVel + bias) / invMassSum;

  // Apply impulse Jᵀ·λ to both bodies
  a.vel = sub(a.vel, scale(n, lambda * a.invMass));
  b.vel = add(b.vel, scale(n, lambda * b.invMass));
}