Сохранение инвариантов: симплектические методы и Верле

На коротком счёте RK4 непобедим по точности. Но запустите его на маятник или орбиту на тысячи периодов — и заметите беду: энергия медленно уплывает. Планета по чуть-чуть приближается к звезде или улетает прочь. Чтобы считать долго, нужны методы особого сорта — те, что берегут не точное значение, а структуру задачи.

Симплектический интегратор — численный метод для гамильтоновых (консервативных) систем, сохраняющий геометрическую структуру фазового пространства. Как следствие, энергия в нём не дрейфует со временем, а лишь колеблется около постоянного значения.

Проблема дрейфа энергии

Рассмотрим консервативную систему — такую, где полная энергия должна сохраняться вечно. Простейший пример — гармонический осциллятор (грузик на пружинке, или маятник при малых колебаниях):

x'' = -x        (это и есть закон Ньютона: ускорение = -сила)
Энергия:  E = 0.5 * (v^2 + x^2),   где v = x'
Точное решение: x = cos(t), и E постоянна всегда.

Запустим на этой системе обычный явный метод Эйлера и проследим за энергией на много периодов. Период колебаний равен 2*pi, то есть около 6.28; пройдём время до 100 — это примерно 16 периодов. Что произойдёт с энергией?

Беда в том, что явный Эйлер для осциллятора систематически добавляет энергию на каждом шаге. Грузик с каждым колебанием раскачивается всё сильнее — численная орбита раскручивается наружу спиралью. Это не случайная ошибка, которая то туда, то сюда: это направленный дрейф, и он накапливается неограниченно. На длинном счёте решение становится физически абсурдным.

Симплектический (полу-неявный) метод Эйлера

Удивительно, но спасает крошечное изменение. В обычном Эйлере мы обновляем координату и скорость по старым значениям:

явный Эйлер:    x_new = x + h*v        (по старой v)
                v_new = v - h*x        (по старой x)

А в полу-неявном (симплектическом) Эйлере мы сначала обновляем скорость, а потом координату — уже по новой скорости:

симплектический Эйлер:  v_new = v - h*x        (сначала скорость)
                        x_new = x + h*v_new    (потом координата, по новой v)

Перестановка двух строк и одна замена v на v_new — вот и вся разница в коде. Но эффект колоссальный: энергия перестаёт дрейфовать и начинает лишь слегка колебаться около константы. Орбита больше не раскручивается. Сравним оба метода на одном прогоне.

import math

# Гармонический осциллятор x'' = -x. Энергия E = 0.5*(v^2 + x^2).
# Точное решение: E постоянна. Старт: x=1, v=0 -> E0 = 0.5.

def energy(x, v):
    return 0.5 * (v * v + x * x)

h = 0.05
steps = 2000          # время от 0 до 100 -> около 16 периодов

# --- Явный метод Эйлера (обновляем по СТАРЫМ значениям) ---
x, v = 1.0, 0.0
E0 = energy(x, v)
for _ in range(steps):
    x_new = x + h * v       # по старой v
    v_new = v - h * x       # по старой x
    x, v = x_new, v_new
E_euler = energy(x, v)

# --- Симплектический (полу-неявный) Эйлер ---
# сначала скорость, потом координата по НОВОЙ скорости
x2, v2 = 1.0, 0.0
for _ in range(steps):
    v2 = v2 - h * x2        # сначала скорость
    x2 = x2 + h * v2        # потом координата, по новой v2
E_symp = energy(x2, v2)

print("E_начало             =", round(E0, 6))
print("E_конец (Эйлер)       =", round(E_euler, 6))
print("E_конец (симплектич.) =", round(E_symp, 6))

Вывод:

E_начало             = 0.5
E_конец (Эйлер)       = 73.745001
E_конец (симплектич.) = 0.510945

Цифры говорят сами за себя. Стартовая энергия — 0.5. У явного Эйлера за 16 периодов она раздулась до 73.7 — почти в 150 раз! Грузик «раскачался» из ничего, что физически невозможно для замкнутой системы. А у симплектического Эйлера энергия осталась 0.511 — отклонение от 0.5 всего на пару процентов, и это отклонение не растёт, а колеблется. Один и тот же шаг h, один и тот же осциллятор — разница только в порядке двух строк обновления.

Метод Верле (Verlet) и leapfrog

Симплектический Эйлер — первого порядка точности, грубоват. Для серьёзных расчётов (молекулярная динамика, небесная механика) используют его старшего брата — метод Верле, он же скоростной Верле или leapfrog («чехарда»). Он тоже симплектический, но имеет второй порядок точности и при этом всё так же дёшев. Для системы x''=a(x) (ускорение зависит только от координаты, как в задаче с потенциалом U: a=-grad U) скоростной Верле выглядит так:

a = ускорение в текущей точке x
x_new = x + h*v + 0.5*h^2 * a            (шаг координаты)
a_new = ускорение в новой точке x_new
v_new = v + 0.5*h*(a + a_new)            (шаг скорости по среднему)

Название «leapfrog» — оттого что в эквивалентной форме координаты и скорости считаются в шахматном порядке, «перепрыгивая» друг через друга на полушаге. Связь с физикой прямая: это просто аккуратная дискретизация второго закона Ньютона F=m*a для консервативной силы. Именно метод Верле десятилетиями держит на плаву симуляции движения молекул и планет — потому что бережёт энергию на миллионах шагов.

Как работает под капотом

Почему перестановка строк творит чудо? Дело в геометрии фазового пространства — плоскости (x, v). Точное решение осциллятора движется по окружности постоянного радиуса (радиус задаёт энергию). Численный метод на каждом шаге чуть искажает эту картину. Явный Эйлер на каждом шаге растягивает фазовый объём — площади раздуваются, и точка по спирали уходит наружу: энергия растёт. Неявный Эйлер, наоборот, сжимает — точка падает внутрь.

Симплектические методы устроены так, что точно сохраняют фазовый объём (это и значит слово «симплектический»). Точка не уходит ни наружу, ни внутрь — она остаётся на почти-окружности. Глубокая теорема говорит больше: симплектический интегратор точно сохраняет не исходную энергию, а энергию слегка возмущённой системы, очень близкой к настоящей. Поэтому численная энергия вечно колеблется в узком коридоре вокруг истинной и никогда не уплывает. Это бесценно для долгого счёта: можно интегрировать орбиту Солнечной системы на миллиарды лет, и она не развалится.

Важная оговорка: симплектические методы не делают каждый отдельный шаг точнее, чем RK4. На коротком интервале RK4 обычно даёт меньшую ошибку положения. Преимущество симплектических методов — качественное и долгосрочное: они берегут инвариант (энергию) и потому не порождают физически невозможного поведения за долгое время. Это пример общего принципа: иногда важнее сохранить структуру задачи, чем минимизировать локальную ошибку.

Частые ошибки

• Считать RK4 универсально лучшим. RK4 точнее на коротком счёте, но не симплектичен: на консервативной системе он медленно дрейфует по энергии. Для долгих симуляций это дисквалифицирующий недостаток.
• Перепутать порядок обновлений. В симплектическом Эйлере критичен порядок: сначала скорость по старой координате, потом координата по новой скорости. Если обновить по-старому — получится обычный явный Эйлер с дрейфом.
• Применять симплектические методы не к тем задачам. Их преимущество — на гамильтоновых (консервативных) системах без трения. На диссипативных задачах (с затуханием) сохранять энергию не нужно и даже вредно — там энергия и должна убывать.
• Ждать монотонного убывания ошибки энергии. Симплектическая энергия не сходится к константе — она колеблется около неё. Это нормально и правильно: важно, что амплитуда колебаний не растёт со временем.

Итоги

• На долгом счёте обычные методы (включая RK4) дрейфуют по энергии: численная орбита раскручивается или схлопывается.
• Симплектические интеграторы сохраняют структуру фазового пространства — энергия колеблется, но не дрейфует.
• Полу-неявный (симплектический) Эйлер отличается от явного порядком обновлений: сначала скорость, потом координата по новой скорости.
• Метод Верле (leapfrog) — симплектический метод второго порядка для x''=a(x), рабочая лошадка молекулярной и небесной динамики.
• Симплектические методы выигрывают не точностью одного шага, а сохранением инварианта на огромных интервалах времени.