Метод Верле и скоростной Верле

Золотой стандарт молекулярной динамики: метод Верле и его скоростная версия.

Метод Верле — симплектический интегратор второго порядка, вычисляющий новую координату по ускорению, а скорость — по среднему ускорению на концах шага; основа молекулярной динамики и небесной механики.

Зачем нужен Верле

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

Скоростной Верле

Самая удобная форма — скоростной Верле (velocity Verlet). За шаг он делает три действия: продвигает координату с учётом текущего ускорения, вычисляет новое ускорение, и обновляет скорость по среднему старого и нового ускорений:

a      = ускорение(x)
x_new  = x + v·dt + ½·a·dt²        (учёт кривизны через a·dt²)
a_new  = ускорение(x_new)          (пересчёт силы в новой точке)
v_new  = v + ½·(a + a_new)·dt      (скорость по среднему ускорению)

Ключевая идея — усреднение ускорения на концах шага. Это устраняет систематический перекос, из-за которого Эйлер копил ошибку. Сравним три метода на осцилляторе за одинаковое время и посмотрим на отклонение энергии от истинного значения 0.5:

import math
def run(method, steps=2000, dt=0.05):
    x, v = 1.0, 0.0
    for _ in range(steps):
        if method == "euler":
            xn = x + v*dt
            vn = v - x*dt
            x, v = xn, vn
        elif method == "symplectic":
            v = v - x*dt
            x = x + v*dt
        elif method == "verlet":
            a = -x
            x = x + v*dt + 0.5*a*dt*dt
            a_new = -x
            v = v + 0.5*(a + a_new)*dt
    return 0.5*v*v + 0.5*x*x
E0 = 0.5
for m in ("euler", "symplectic", "verlet"):
    E = run(m)
    print(f"{m:11s}: E = {E:.4f}  (отклонение {100*(E-E0)/E0:+.1f}%)")
print(f"Точное значение E0 = {E0:.4f}")

Вывод:

euler      : E = 73.7450  (отклонение +14649.0%)
symplectic : E = 0.5109  (отклонение +2.2%)
verlet     : E = 0.4999  (отклонение -0.0%)
Точное значение E0 = 0.5000

Результат говорит сам за себя. За 2000 шагов явный Эйлер раздул энергию в 147 раз — симуляция полностью разрушена. Полунеявный Эйлер удержал отклонение в 2%. Верле дал почти идеальные 0.4999 — практически точное сохранение. Один и тот же шаг, разные методы, несопоставимое качество.

Почему Верле точнее

Член ½·a·dt² в обновлении координаты учитывает кривизну траектории — то, как ускорение искривляет путь за шаг. Эйлер этот член игнорирует, поэтому идёт по прямой касательной. Усреднение ускорений ½(a+a_new) для скорости — это правило трапеций, точное до второго порядка. Вместе они дают симплектический метод второго порядка почти даром.

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

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

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

  • Забыть пересчитать ускорение в новой точке. Скоростной Верле требует двух значений ускорения за шаг — старого и нового.
  • Использовать одно и то же a для координаты и скорости. Скорость обновляется по среднему старого и нового ускорений, иначе теряется второй порядок.
  • Менять dt посреди симуляции. Классический Верле чувствителен к постоянству шага; для переменного шага есть отдельные схемы.

Итоги

  • Верле — симплектический интегратор второго порядка, основа МД и небесной механики.
  • Скоростной Верле усредняет ускорение на концах шага для скорости.
  • Член ½·a·dt² учитывает кривизну траектории — отсюда точность.
  • На одном шаге Верле сохраняет энергию несравнимо лучше Эйлера.
Проверьте себя
1. Какой порядок точности у метода Верле?
AНулевой
BПервый
CВторой
DДесятый
2. По какому ускорению скоростной Верле обновляет скорость?
AПо старому ускорению
BПо новому ускорению
CПо среднему старого и нового ускорений
DПо нулевому
3. Зачем в обновлении координаты член ½·a·dt²?
AДля красоты формулы
BОн учитывает кривизну траектории, которую Эйлер игнорирует
CЧтобы увеличить шаг
DЧтобы обнулить скорость