Verlet, Leapfrog and RK4: integrators for physics

1. The Euler problem

Метод Ейлера — найпростіший чисельний інтегратор: взяти поточне значення похідної і перемножити на крок часу.

Проблема в тому, що Ейлер не зберігає повну механічну енергію системи. У Гамільтонових системах (маятник, орбіти, пружини) метод Ейлера систематично додає або забирає енергію на кожному кроці. Результат: пружина нескінченно розтягується, а планета закручується в спіраль замість того, щоб залишатися на орбіті.

Метод Ейлера має першосортну локальну похибку O(dt²) і глобальну похибку O(dt). Щоб зменшити похибку вдвічі, треба зменшити крок вдвічі — і вдвічі більше роботи. Це не вирішує проблему нестабільності.

2. Verlet method

Люк Верле запропонував метод у 1967 році для молекулярної динаміки. Ключова ідея: обчислювати позицію через дві попередні позиції, уникаючи прямого зберігання швидкості.

Виводимо з розкладу Тейлора: запишемо x(t+dt) і x(t-dt), потім складемо:

Де a(t) = F(t)/m — прискорення у поточний момент. Швидкість у цій формулі явно не фігурує, але може бути отримана апроксимацією:

Verlet method є симплектичним — він точно зберігає симплектичну форму фазового простору, що для Гамільтонових систем означає збереження енергії на довгих горизонтах. Його глобальна похибка — O(dt²), як і в Ейлера, але без зростаючого зміщення енергії.

3. Velocity Verlet

Velocity Verlet — модифікований метод Верле, де швидкість зберігається явно. Алгоритм двокроковий:

Зверніть увагу: між кроками потрібно перерахувати a(t+dt) з нової позиції. Тобто сила обчислюється двічі за ітерацію. Але Velocity Verlet точно підтримує v — що критично для симуляцій тканини, де пружинна сила залежить від v (демпфування).

Velocity Verlet — стандарт для рушіїв типу Cannon-es, Bullet, Havok завдяки точному v і відмінній стабільності.

4. Leapfrog integrator

Leapfrog («жаб'ячий стрибок») — ще одна форма того ж симплектичного методу. Позиція та швидкість зсунуті на половину кроку відносно одна одної, як два жабяти що стрибають по черзі:

Алгоритм обчислює прискорення лише один раз за крок (на відміну від Velocity Verlet), тому він вдвічі дешевший за обчисленнями. Саме тому Leapfrog є стандартом для N-тіл та молекулярної динаміки.

Drift-Kick-Drift (DKD)

Leapfrog можна записати у так званій DKD-формі — «дрейф-удар-дрейф» — щоб стартувати з цілого часового кроку:

5. Runge-Kutta 4th order (RK4)

RK4 — найвідоміший загальний метод чисельного інтегрування ОДУ. Він оцінює похідну у чотирьох точках всередині кроку і береже зважене середнє:

RK4 має глобальну похибку O(dt⁴) — при зменшенні кроку вдвічі похибка зменшується у 16 разів. Він ідеальний для хаотичних систем (атрактор Лоренца, подвійний маятник), де потрібна висока точність траєкторії за відносно короткий час.

RK4 обчислює f() чотири рази за крок — це вчетверо дорожче за Leapfrog. Але він незамінний для жорстких систем (stiff ODEs) і там, де траєкторія важлива сама по собі (подвійний маятник, атрактор Лоренца).

6. Method comparison

7. When to use which

Real-time game physics (60 FPS)

Для ігрових рушіїв та інтерактивних симуляцій з великим dt (1/60 с = 16.7 мс) — Velocity Verlet або Symplectic Euler. Вони симплектичні, дешеві і достатньо точні. Cannon-es використовує Symplectic Euler зі split impulse для тіл.

Cloth and soft bodies

Velocity Verlet: демпфування пружин залежить від швидкості — потрібне явне v. Position-Based Dynamics (PBD) використовує скорочену форму Верле: оновлює позиції без пружних сил, потім виконує constraint projection.

N-body and molecular dynamics

Leapfrog / Störmer-Verlet. Мільярди кроків у MD-симуляціях — тому дешевизна (1 обчислення F/крок) та збереження енергії критично важливі.

Chaotic attractors, double pendulum

RK4. Траєкторія важлива сама по собі — потрібна висока точність на кожному кроці. Та й горизонт відносно короткий (секунди або хвилини).

8. Псевдокод

Velocity Verlet (JS)

// Velocity Verlet — тканина, маятник
function stepVelocityVerlet(particles, dt) {
  for (const p of particles) {
    // Крок 1: оновити позицію
    p.x += p.vx * dt + 0.5 * p.ax * dt * dt;
    p.y += p.vy * dt + 0.5 * p.ay * dt * dt;

    // Зберегти старе прискорення
    const ax0 = p.ax, ay0 = p.ay;

    // Крок 2: перерахувати сили з нової позиції
    computeForces(particles);

    // Крок 3: оновити швидкість через середнє прискорення
    p.vx += 0.5 * (ax0 + p.ax) * dt;
    p.vy += 0.5 * (ay0 + p.ay) * dt;
  }
}

Leapfrog DKD (JS)

// Leapfrog — N-тіла, планети
function stepLeapfrog(bodies, dt) {
  const hdt = dt / 2;

  // Drift (D): половина кроку позиція
  for (const b of bodies) {
    b.x += b.vx * hdt;
    b.y += b.vy * hdt;
  }

  // Kick (K): повний крок швидкість
  computeGravity(bodies);            // → b.ax, b.ay
  for (const b of bodies) {
    b.vx += b.ax * dt;
    b.vy += b.ay * dt;
  }

  // Drift (D): ще половина кроку позиція
  for (const b of bodies) {
    b.x += b.vx * hdt;
    b.y += b.vy * hdt;
  }
}

RK4 для ОДУ (JS)

// RK4 — атрактор Лоренца, подвійний маятник
function rk4(state, deriv, dt) {
  const k1 = deriv(state);
  const k2 = deriv(add(state, scale(k1, dt/2)));
  const k3 = deriv(add(state, scale(k2, dt/2)));
  const k4 = deriv(add(state, scale(k3, dt)));

  // Зважене середнє: 1/6 · (k1 + 2k2 + 2k3 + k4)
  return add(state, scale(
    add(k1, add(scale(k2, 2), add(scale(k3, 2), k4))),
    dt / 6
  ));
}

// Допоміжні функції для векторних операцій
const add   = (a, b) => a.map((v, i) => v + b[i]);
const scale = (a, s) => a.map(v => v * s);