Leapfrog symplecticity: a proof, and why it bounds energy error forever

1. What "symplectic" means

Гамільтонова система еволюціонує у фазовому просторі (q, p) (координати, імпульси). Симплектична форма — це білінійна антисиметрична 2-форма ω = dq ∧ dp, яка вимірює орієнтовану площу в фазовому просторі. Потік φₜ точних рівнянь Гамільтона зберігає цю форму точно: φₜ* ω = ω для будь-якого t.

Тут J = ∂(q', p')/∂(q, p) — якобіан відображення одного кроку в фазовому просторі. Це узагальнення умови det J = 1 (збереження об'єму, теорема Ліувілля) — симплектичність строгіша: вона зберігає орієнтовану площу для кожної 2-площини, не лише повний об'єм.

2. Leapfrog as a composition of maps

Розкладемо Leapfrog (DKD-форма) на три елементарні відображення фазового простору:

Ключова ідея доведення: кожне з трьох відображень окремо симплектичне, а композиція симплектичних відображень завжди симплектична (це випливає з ланцюгового правила для якобіанів: J_composite = J₁·J₂·J₃, і якщо кожне Jᵢᵀ·Ω·Jᵢ = Ω, то й добуток задовольняє те саме).

3. Computing the Jacobian of Drift and Kick

Розглянемо 1D-систему (q, p) для простоти; узагальнення на багато вимірів — покомпонентне. Якобіан Drift-кроку:

Якобіан Kick-кроку (де сила F(q) залежить лише від позиції — типово для консервативних потенціалів):

Обидві матриці мають детермінант точно 1 — необхідна (але не достатня) умова симплектичності в 1D.

4. The proof: each map preserves ω

У 1D симплектична умова спрощується до det(J) = 1 (в одновимірному фазовому просторі 2×2 симплектична форма — це просто орієнтована площа). Перевіримо для Drift:

Аналогічно для Kick:

Оскільки детермінант добутку матриць дорівнює добутку детермінантів:

Це доводить, що кожен крок Leapfrog зберігає орієнтовану площу фазового простору точно — незалежно від dt, від нелінійності сили F(q), і без жодного наближення. На відміну від Euler чи RK4, тут немає похибки, яка накопичується — умова виконується рівно, аналітично.

5. Backward error analysis

Найглибший результат теорії симплектичних інтеграторів (Hairer, Lubich, Wanner) — backward error analysis: симплектичний метод точно розв'язує «збурений» гамільтоніан H̃ = H + h²·H₂ + h⁴·H₄ + … для деякої обчислюваної H₂, H₄, … Оскільки H̃ — це справжня консервативна система (сама симплектична), її енергія зберігається точно.

«Назавжди» — це трохи перебільшення: похибка обмежена лише на експоненційно довгому проміжку часу (за теоремою KAM/backward-error), але для практичних симуляцій (мільйони — мільярди кроків N-тіл) це фактично невідрізнено від «завжди».

6. Why RK4 is not symplectic

RK4 обчислює зважене середнє чотирьох оцінок похідної в різних точках — це відображення фазового простору не факторизується у послідовність канонічних (симплектичних) перетворень. Якобіан RK4-кроку не задовольняє Jᵀ·Ω·J = Ω за жодного вибору dt (окрім тривіального dt → 0). Це можна перевірити прямим обчисленням якобіана для простого гармонічного осцилятора — det(J_RK4) ≠ 1 при скінченному h.

7. Consequences compared

8. Псевдокод

Verifying det(J) numerically (JS)

// Чисельна перевірка det(J) = 1 для одного кроку Leapfrog
// на простому гармонічному осциляторі F(q) = −k·q
function leapfrogJacobianDeterminant(k, m, h) {
  // Drift: q' = q + p*(h/2)/m, p' = p
  const Jd = [[1, h / (2 * m)], [0, 1]];
  // Kick: p' = p + F(q)*h = p − k*q*h, q' = q  → F'(q) = −k
  const Jk = [[1, 0], [−k * h, 1]];

  matMul2 = (A, B) => [
    [A[0][0]*B[0][0] + A[0][1]*B[1][0], A[0][0]*B[0][1] + A[0][1]*B[1][1]],
    [A[1][0]*B[0][0] + A[1][1]*B[1][0], A[1][0]*B[0][1] + A[1][1]*B[1][1]]
  ];

  const J = matMul2(matMul2(Jd, Jk), Jd);
  const det = J[0][0]*J[1][1] − J[0][1]*J[1][0];
  return det; // always exactly 1, for any k, m, h
}