Численный метод: симплектический Эйлер

Производная – это скорость изменения. Если величина $f$ сейчас равна $f(t)$ и меняется со скоростью $f'(t)$, то через маленькое время $\Delta t$ она будет примерно равна

$$f(t + \Delta t) \approx f(t) + f'(t)\,\Delta t$$

Это и есть метод Эйлера: заменяем точное решение касательной и проходим по касательной короткий отрезок. Чем короче отрезок, тем меньше мы ошибёмся.

На C++ это записывается ещё проще, потому что там есть операция «прибавить к себе»:

alpha += v / l * dt;         // угол вырос на (скорость × время)
v     += -g * sin(alpha) * dt;   // скорость выросла на (ускорение × время)

Эти две строки задают один шаг используемой схемы.

В каком порядке их выполнять?

Порядок присваиваний определяет, какие значения – старые или уже обновлённые – используются в правой части. Рассмотрим два варианта.

Различие ровно в одном: во втором варианте, когда мы считаем новую скорость, переменная alpha уже обновлена. Мы берём силу не там, где маятник был, а там, где он оказался. В коде задачи стоит вариант Б.

У настоящего маятника без трения полная энергия постоянна:

$$E = \frac{v^2}{2} + (1 - \cos\alpha)$$

(в безразмерных величинах, где $m = l = g = 1$; первое слагаемое – кинетическая энергия, второе – потенциальная, отсчитанная от нижней точки).

Контроль энергии позволяет оценить накопление численной ошибки. Поведение двух схем различается:

Обе схемы имеют первый порядок глобальной точности, но по-разному ведут себя на больших временах. Явный Эйлер не сохраняет геометрическую структуру гамильтоновой системы и для этой задачи даёт заметный дрейф энергии. Симплектический Эйлер сохраняет симплектическую структуру; поэтому при достаточно малом постоянном шаге ошибка энергии обычно остаётся ограниченной.

Симплектическая схема не сохраняет точную энергию и не гарантирует точного совпадения траектории после периода. Её преимущество – отсутствие систематического ухода энергии на рассматриваемом интервале.

Строгое объяснение формулируется в терминах гамильтоновой механики. На фазовой плоскости точное движение сохраняет площадь. Симплектический Эйлер сохраняет соответствующее дискретное геометрическое свойство, а явный Эйлер – нет. Поэтому одинаковая локальная точность не означает одинакового поведения на длинном интервале времени.

Схемы, устроенные так, что они сохраняют структуру уравнений, называются симплектическими, и в задачах о движении под действием сил пользуются именно ими. Расчёт орбит на миллионы лет вперёд ведут схемами того же семейства, просто более точными.

В следующей задаче о движении в поле тяготения используется та же схема; там различие методов можно проверить по энергии и форме орбиты.

В программе используется шаг dt = 0.05.

Период малых колебаний в безразмерных единицах $T = 2\pi \approx 6.28$. Значит на один период приходится $6.28 / 0.05 \approx 126$ шагов на одно полное колебание. Этого хватает и для гладкой картинки, и для точности.

Для первоначальной оценки полезно иметь не менее 20–30 шагов на наименьший характерный период. Это не универсальная гарантия точности: результат следует проверять повторным расчётом с меньшим шагом.

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

В этой программе шаг зашит в код строкой float dt = 0.05f; и в пульт не выведен. Чтобы поэкспериментировать с ним, придётся править файл и пересобирать программу – либо вывести его в пульт самостоятельно.

Метод Эйлера относится к схемам первого порядка. Локальная ошибка за один шаг пропорциональна $\Delta t^2$, а за всё время расчёта – $\Delta t$; чтобы улучшить ответ вдесятеро, нужно вдесятеро больше шагов.

Для более точных расчётов применяют, например, метод Верле (симплектический метод второго порядка) и методы Рунге – Кутты. Классический RK4 имеет четвёртый порядок, но не является симплектическим.

В этой работе симплектический Эйлер выбран из-за простой записи: на нём можно отдельно исследовать влияние шага и порядка обновления переменных.