Задания

Изменяя $V_y$ в пульте, определите значение первой космической скорости.

При радиусе $x_0 = 2.5$ подбирайте скорость так, чтобы орбита стала круговой – расстояние $R$ перестало меняться. Показание «Расстояние R» вам в помощь: у круговой орбиты оно постоянно.

Сравните найденное значение с тем, что даёт формула круговой скорости из первой статьи.

Затем проверьте формулу ещё на трёх радиусах – 1, 4 и 9. Как падает скорость с ростом радиуса – во сколько раз она меняется, когда радиус растёт в девять раз? Для радиусов 4 и 9 предварительно расширьте границы графика в create_plot.

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

Изменяя $V_y$, определите вторую космическую скорость – минимальную, при которой планета улетает и не возвращается.

Сравните с формулой второй космической из первой статьи.

На симуляции граница определяется трудно: при скорости чуть меньше критической планета уходит очень далеко и возвращается через долгое время. Придётся расширить границы графика в create_plot – например до $\pm 30$ – и запастись терпением.

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

Постройте графики зависимости кинетической, потенциальной и полной энергий от времени.

Формулы для $K$, $U$ и $E$ – на единицу массы планеты – в первой статье.

Заведите второй график и рисуйте на нём три кривые с разными подписями. Понадобится переменная времени – её в этой программе нет, придётся добавить (float t = 0.f; и t += dt; в расчётной функции).

Что проверить:

1. в какой точке орбиты кинетическая энергия наибольшая, а потенциальная наименьшая? Сверьтесь с показанием «Расстояние R»: совпадают ли эти моменты;
2. полная энергия должна быть постоянной. Насколько она колеблется? Увеличьте шаг вдвое и посмотрите, как вырастет размах;
3. вычислите отдельно средние значения $\langle K\rangle$ и $\langle U\rangle$ за полный оборот. Для связанной кеплеровой орбиты теорема вириала даёт $2\langle K\rangle=-\langle U\rangle$. Обратите внимание: это отношение средних, а не среднее мгновенного отношения $K/|U|$.

Часть первая, на бумаге. Как выражается момент импульса $L$ через компоненты $x$, $y$, $v_x$, $v_y$? Вспомните определение $\vec{L} = \vec{r} \times \vec{p}$ и то, что для плоского движения от векторного произведения остаётся одна компонента.

Часть вторая, в программе. Вычислите $L$ на каждом шаге, выведите показанием и постройте график.

Проверьте:

1. сохраняется ли $L$? Насколько точно по сравнению с энергией?
2. поставьте круговую орбиту. Чему равен $L$ и совпадает ли с $r v$?
3. на вытянутой орбите сравните скорость в ближайшей и дальней точках. Выполняется ли $r_1 v_1 = r_2 v_2$? Это и есть второй закон Кеплера в самом простом виде.

Проверьте выполнение третьего закона Кеплера, $T_1^2 / T_2^2 = a_1^3 / a_2^3$.

Нужно измерить период и большую полуось для нескольких разных орбит.

Период проще всего ловить так же, как в маятнике, – по смене знака координаты $y$ с минуса на плюс (планета вернулась на ось $x$ с той же стороны):

if (y_prev < 0.f && y >= 0.f) {
    set_float_param("Период", t - t_prev);
    t_prev = t;
}
y_prev = y;

Большая полуось считается по $R_{max}$ и $R_{min}$ (определение в первой статье). Оба расстояния стоит копить в программе, а не ловить глазами.

Соберите таблицу из четырёх-пяти орбит и проверьте, что $T^2/a^3$ у всех одинаково. Затем посчитайте, какое значение этой величины предсказывает безразмерный третий закон Кеплера из первой статьи. Сошлось?

Задания на численный метод, а не на физику.

1. Поставьте очень вытянутую орбиту: $x_0 = 2.5$, $V_y = 0.15$. Планета будет проходить у самой звезды – там, где по третьей статье и накапливается почти вся ошибка шага. Что при этом происходит с формой орбиты за десяток витков? Измерьте эффект: на сколько примерно поворачивается эллипс за один виток?
2. Уменьшайте шаг (dt в начале файла) и повторяйте. Во сколько раз падает поворот при уменьшении шага вдвое? При каком шаге он перестаёт быть заметным?
3. Настоящие орбиты планет тоже поворачиваются – прецессия перигелия Меркурия составляет 43 угловые секунды за век и объясняется общей теорией относительности. Как отличить настоящую прецессию от численной?

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

Добавьте ещё одну планету, не взаимодействующую с первой.

Заведите второй набор из четырёх переменных состояния и повторите для них шаг метода. Пусть у второй планеты будут свои начальный радиус и скорость в пульте.

Затем – главное:

1. как выглядит движение второй планеты при наблюдении её с первой? Постройте график разности координат $(x_2 - x_1, y_2 - y_1)$ – это и есть вид с первой планеты;
2. чем отличаются картины, наблюдаемые с внутренней и внешней планет?

Петли на относительной траектории соответствуют попятному движению планеты при наблюдении с другой движущейся планеты. Такой график связывает расчёт в гелиоцентрической системе с наблюдаемым движением на небесной сфере.

Введите взаимодействие между планетами.

Теперь каждая планета притягивается и к звезде, и к другой планете. Ускорение первой:

Здесь $\alpha_2$ – параметр тяготения второй планеты, а $R_{12}$ – расстояние между планетами. Для второй планеты формула симметрична.

Возьмите массу планеты в десять раз меньше массы звезды, а массу спутника – в десять раз меньше массы планеты. Тогда $\alpha_{planet} = 0.1$, $\alpha_{satellite} = 0.01$.

Это не отношения масс Солнца, Земли и Луны: настоящая Земля легче Солнца в 330 тысяч раз, и на экране такой спутник был бы не виден. Числа выбраны так, чтобы взаимодействие планеты со спутником было заметно за десяток витков. Плата за наглядность: планета в $0.1$ массы звезды – это уже слишком много, чтобы считать звезду неподвижной, поэтому модель остаётся качественной, а тела здесь называются планетой и спутником, а не Землёй и Луной.

Начальная оценка условий для последующего подбора:

x_planet = 2.5,  Vy_planet = sqrt(alpha / x_planet)
x_satellite = 2.7
Vy_satellite ≈ sqrt(alpha_planet / (x_satellite - x_planet)) + Vy_planet

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

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

Подберите параметры так, чтобы система прожила хотя бы десяток оборотов и не распалась. Что происходит, если сделать спутник тяжелее?

Рассмотрите движение частицы в поле вида $F \sim R^n$. При каких $n$ получаются замкнутые орбиты?

В программе меняется одна строка:

a = alpha * pow(R, n);   // было alpha / (R * R), то есть n = -2

Заведите ручку для $n$ и попробуйте значения от $-3$ до $2$. Для каждого значения используйте несколько невырожденных связанных начальных условий: круговая орбита возможна и в других полях и сама по себе не проверяет теорему.

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

Как выглядит замкнутая орбита в каждом из найденных полей и где в ней сидит центр притяжения – в фокусе эллипса или в его середине?

Осторожно: при $n$ около $-3$ и близком подлёте расчёт разваливается быстро. Это уже знакомая проблема шага, а не физика.

Учтите влияние на спутник торможения в атмосфере. Простейший способ – ввести силу сопротивления $\vec{F} = -k v \vec{v}$, пропорциональную квадрату скорости и направленную против движения.

float v = sqrt(v_x * v_x + v_y * v_y);
a_x += -k * v * v_x;
a_y += -k * v * v_y;

В приведённом коде слагаемое добавляется к ускорению, поэтому $k$ понимается как коэффициент сопротивления, уже делённый на массу. Начните со значения порядка 0.001 и проверьте зависимость результата от шага по времени.

Запустите почти круговую орбиту и дайте ей медленно снижаться. Постройте графики $K$, $U$ и $E$ и проследите за каждой величиной отдельно.

Сила сопротивления совершает отрицательную работу, значит полная механическая энергия обязана падать. А что происходит со скоростью спутника? Ответ, скорее всего, вас удивит – объясните его.

Для заметно вытянутой орбиты мгновенная кинетическая энергия ведёт себя сложнее, поэтому вывод делайте по усреднённым за оборот значениям.

Переделайте шаг в явный: обе величины считаются по старым значениям.

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

// потом обе величины
x += v_x * dt;   y += v_y * dt;
v_x += a_x * dt; v_y += a_y * dt;

Запустите круговую орбиту (скорость возьмите из задания 1) и подождите десяток витков. Что происходит с радиусом – измерьте, на сколько процентов он уходит. Затем верните симплектическую схему и повторите измерение на той же орбите.

Вопросы:

1. постройте график полной энергии для обеих схем. Чем принципиально отличаются кривые – не величиной ошибки, а её поведением?
2. уменьшите шаг вдесятеро. Спираль исчезла или просто стала медленнее? Что это говорит о том, можно ли «победить» явную схему мелким шагом?
3. симплектическая схема тоже не сохраняет энергию точно. Остаётся ли численная орбита строго замкнутой? Измерьте смещение перицентра после нескольких оборотов.

Не забудьте вернуть правильную схему обратно.