Численный метод: движение по плоскости – Спецкурс

Состояние – четыре числа

У маятника состояние задавалось двумя числами: угол и скорость. Здесь тело движется по плоскости, и состояние – четыре числа:

float x, y;                    // где планета
float v_x = 0.0f, v_y = 0.5f;  // куда и как быстро летит

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

Всё остальное – расстояние, ускорение, энергия – из этих четырёх чисел вычисляется и хранить его не нужно.

Шаг численного метода

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

// Шаг первый: координаты по старым скоростям
x += v_x * dt;
y += v_y * dt;

// Сила в НОВОЙ точке
R = sqrt(x * x + y * y);
a = alpha / (R * R);
a_x = -x / R * a;
a_y = -y / R * a;

// Шаг второй: скорости по этой силе
v_x += a_x * dt;
v_y += a_y * dt;

Порядок здесь – не оформление, а сам метод:

координаты сдвигаются по старым скоростям;
по новым координатам считается расстояние и ускорение;
скорости обновляются по этому ускорению.

Если сначала обновить скорость по силе в старой точке, а затем координату по новой скорости, получится второй вариант симплектического Эйлера. Явный Эйлер требует отдельно вычислить и новую координату, и новую скорость из одного старого состояния, а затем одновременно заменить переменные.

Сравнение с явным Эйлером

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

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

Куда именно уедет расчёт, зависит от схемы, шага и начальных данных: универсальной монотонности у явного Эйлера нет, и предсказать знак дрейфа заранее нельзя – его меряют.

Локальная ошибка обеих схем имеет порядок $\Delta t^2$, а глобальная – порядок $\Delta t$. Различается долговременное поведение: симплектический метод сохраняет геометрическую структуру гамильтоновой системы, но не точную энергию и не точную орбиту.

Проверка расчёта законами сохранения

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

В идеальной задаче двух тел сохраняются две величины – полная энергия и момент импульса (статья 1):

$$E = \frac{v^2}{2} - \frac{\alpha}{R}$$

Численная схема о них не знает – она знает только уравнения движения. Значит, посчитав $E$ и $L$ на каждом шаге и посмотрев, постоянны ли они, мы измеряем качество расчёта напрямую.

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

Именно поэтому графиков энергии и момента импульса в программе нет, хотя добавить их несложно.

Ограничения постоянного шага

Шаг dt = 0.02 подобран под орбиту по умолчанию: период около 15.4, то есть примерно 770 шагов на оборот. Запас большой.

Но у этой задачи есть слабое место, которого не было у маятника: вблизи звезды всё происходит быстро. Ускорение растёт как $1/R^2$, и на вытянутой орбите планета проходит перицентр за небольшое число шагов. Ошибка сосредоточена главным образом в этой части орбиты.

Значит, качество расчёта на вытянутой орбите решается в тех нескольких десятках шагов, что планета проводит у звезды, – а форму орбиты мы смотрим на всём витке. На что это похоже на экране, видно сразу, стоит поставить $V_y$ заметно меньше круговой скорости.

Один из способов уточнить прохождение перицентра – уменьшить шаг. Адаптивный шаг делает его малым в области быстрого движения и крупным вдали от центра. Однако произвольное изменение шага нарушает симплектические свойства обычной схемы Эйлера; для таких расчётов нужен специально построенный адаптивный метод либо расчёт с несколькими постоянными шагами.

Смягчение поля вблизи центра

При $R \to 0$ ускорение $\alpha/R^2$ обращается в бесконечность. Это не численная проблема, а свойство самой формулы: точечная масса создаёт бесконечное поле в своей точке.

Физически звезда имеет конечный размер, и внутри неё формула $1/R^2$ для точечной массы неприменима. Если расчётная точка пройдёт слишком близко к центру, дискретная схема даст очень большое изменение скорости за один шаг – планету выбросит в произвольную сторону с произвольной скоростью.

В программе стоит грубая защита: расстояние снизу ограничено числом 1e-4 (как это пишется на C++ – в статье о приёмах языка). Она предотвращает деление на слишком малое число, но меняет модель: если планета до ограничителя добралась, расчёт в этом месте всё равно неверен, и лечить надо шаг по времени или саму модель.

Правильные решения два: специально построенный адаптивный метод или смягчение потенциала – бесконечность убирают из самой формулы, добавляя в знаменатель постоянную $\varepsilon$:

$$a_x = -\frac{\alpha x}{(R^2 + \varepsilon^2)^{3/2}}$$

и так же для $a_y$. При $R \gg \varepsilon$ это прежнее тяготение, а в центре ускорение конечно. Приём стандартный для расчётов, где частицы могут сблизиться вплотную, – в задаче о молекулярной динамике он встретится снова.