суббота, 28 декабря 2013 г.

Движение двух гравитирующих тел: численный и аналитический подходы


Этот пост навеян постом "Задача N тел".
Посвящается: школьника и студентам.

Задача ставится так: есть два точечных тела с массами m_{1} и m_{1}, которые притягиваются друг к другу с силой гравитационного притяжения \vec{F_{12}}=G\frac{m_{1}m_{2}}{r^{3}}\vec{r}_{12}.
(+ читать далее...)
В зависимости от начальных условий (т.е. от значений координат и скоростей этих тел в начальный момент времени), они будут двигаться по замкнутым траекториям (эллипсам), или же разлетятся и никогда больше не встретятся. Надо составить программу, в которую можно вводить какие угодно начальные условия, а она будет рисовать траектории движения этих двух тел. Причем, рисовать эти траектории программа должна двумя способами: по формуле - аналитическому решению задачи двух тел, и численно рассчитав эту траекторию.
Построив такую программу, нужно проанализировать выдаваемые ей результаты: при каких условиях численное решение будет давать большую ошибку, и траектории, рассчитанные аналитически и численно не будут совпадать.

План решения задачи:
1) Получить аналитическое решение задачи двух тел.
2) Написать в Mathcad полученные формулы. Рисовать по ним траектории движения для различных начальных условий.
3) Построить траектории движения численно.
4) Сравнить полученные результаты, сделать выводы.

1) Решение задачи двух тел есть во многих книгах. Например, в первом томе Ландау, Лифшиц, или в этой книге на с.367. Если коротко, то движение двух тел сводят к движению одного тела в определенном потенциале. Поняв как движется это тело, сразу становится понятно, как движутся два исходных тела.
Нажмите сюда чтобы свернуть или развернуть выкладки.

2) Полученные формулы записываем в систему Mathcad и строим графики. Тут можно посмотреть как это будет выглядеть, а тут - загрузить файл .mcdx.
В зависимости от начальных условий эксцентриситет может быть положительным или отрицательным. Его значение по модулю характеризует степень отклонения траектории от окружности:
  • е=0 - окружность
  • 0<е<1 - эллипс
  • е=1 - парабола
  • е>1 - гипербола
А знак эксцентриситета влияет лишь на расположение траектории:
Будем рассматривать следующие начальные условия:
x0 = 10;
y0 = 10;
Vx0 = -0.03;
Vy0= -0.01.
Эксцентриситет при этом равен: е = -0.997.

3) Для численного моделирования запишем в Mathcad систему, описывающую движение того одного тела, к которому можно свести движение двух тел:
\left\{ { \dot{x}=v_{x}\\\dot{v}_{x}=-\gamma\frac{x}{(x^{2}+y^{2})^{\frac{3}{2}}}\\\dot{y}=v_{y}\\ \dot{v}_{y}=-\gamma\frac{y}{(x^{2}+y^{2})^{\frac{3}{2}}} } \right.
Численно решаем эту систему методом Rkfixed для начальных условий, приведенных выше. Убедимся, что она (справа) совпадает с траекторией, построенной по формуле (слева):

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

4) Со временем траектория, рассчитанная численно, начнет давать все большую и большую погрешность:
На последней картинке тела вообще разлетелись в разные стороны.
Как хорошо видно ниже, на графике изменения х-ой составляющей скорости от времени, проблема тут в резких поворотах, просчитать которые довольно сложно, т.к. функция там меняется довольно быстро. Подробнее о проблемах интегрирования рассматриваемых систем читайте тут.

Систему можно решить и другими методами: Stiffb и Stiffr. Они более приспособлены для таких случаев.

Выводы: Задача о движении двух гравитирующих тел имеет точное решение. Решая эту задачу численно, можно "на практике" проверить состоятельность различных методов численного интегрирования.

Замечание: Результаты численного моделирования задачи трех тел (более сложный случай изложенной выше задачи) изложены тут.