В зависимости от начальных условий (т.е. от значений координат и скоростей этих тел в начальный момент времени), они будут двигаться по замкнутым траекториям (эллипсам), или же разлетятся и никогда больше не встретятся. Надо составить программу, в которую можно вводить какие угодно начальные условия, а она будет рисовать траектории движения этих двух тел. Причем, рисовать эти траектории программа должна двумя способами: по формуле - аналитическому решению задачи двух тел, и численно рассчитав эту траекторию.
Построив такую программу, нужно проанализировать выдаваемые ей результаты: при каких условиях численное решение будет давать большую ошибку, и траектории, рассчитанные аналитически и численно не будут совпадать.
План решения задачи:
1) Получить аналитическое решение задачи двух тел.
2) Написать в
Mathcad полученные формулы. Рисовать по ним траектории движения для различных начальных условий.
3) Построить траектории движения численно.
4) Сравнить полученные результаты, сделать выводы.
1) Решение задачи двух тел есть во многих книгах. Например, в первом томе Ландау, Лифшиц, или в
этой книге на с.367. Если коротко, то движение двух тел сводят к движению одного тела в определенном потенциале. Поняв как движется это тело, сразу становится понятно, как движутся два исходных тела.
Нажмите сюда чтобы свернуть или развернуть выкладки.
- Перейти в систему координат с центром в их
барицентре. Благодаря такому переходу формулы упрощаются, и задача сводится к движению одного тела в центрально-симметричном поле. Его радиус-вектор
так связан с радиус-векторами рассматриваемых тел
и
:
- Получается, что чтобы найти траекторию движения двух гравитирующих тел, достаточно найти траекторию всего одного тела в центрально-симметричном поле. В декартовых координатах можно записать:
Для решения задачи двух тел удобно перейти в полярные координаты 
и
:
В этих координатах вектор ускорения равен:
- Приравняем перпендикулярные к радиус-вектору компоненты в декартовых и полярных координатах, получим:
где h - постоянная.
- Приравняем радиальные компоненты в декартовых и полярных координатах, получим для
уравнение:
чьим решением будет:
откуда получаем, что:
где 
.
- Выразим постоянные
и e через начальные условия. Для этого рассмотрим систему:
и выразим из нее
и e:
Осталось выразить h,
,
. Для этого обозначим за a угол между радиус-вектором 
и
,
а за b - угол между скоростью
и
в начальный момент времени. Тогда:
2) Полученные формулы записываем в систему Mathcad и строим графики.
Тут можно посмотреть как это будет выглядеть, а
тут - загрузить файл .mcdx.
В зависимости от начальных условий эксцентриситет может быть положительным или отрицательным. Его значение по модулю характеризует степень отклонения траектории от окружности:
- е=0 - окружность
- 0<е<1 - эллипс
- е=1 - парабола
- е>1 - гипербола
А знак эксцентриситета влияет лишь на расположение траектории:
Будем рассматривать следующие начальные условия:
x0 = 10;
y0 = 10;
Vx0 = -0.03;
Vy0= -0.01.
Эксцентриситет при этом равен: е = -0.997.
3) Для численного моделирования запишем в Mathcad систему, описывающую движение того одного тела, к которому можно свести движение двух тел:
Численно решаем эту систему методом Rkfixed для начальных условий, приведенных выше. Убедимся, что она (справа) совпадает с траекторией, построенной по формуле (слева):
Посмотреть как выглядит программа можно
тут, скачать ее можно
тут или
тут, где она объединена с аналитическим решением.
4) Со временем траектория, рассчитанная численно, начнет давать все большую и большую погрешность:
На последней картинке тела вообще разлетелись в разные стороны.
Как хорошо видно ниже, на графике изменения х-ой составляющей скорости от времени, проблема тут в резких поворотах, просчитать которые довольно сложно, т.к. функция там меняется довольно быстро. Подробнее о проблемах интегрирования рассматриваемых систем читайте
тут.
Систему можно решить и другими методами:
Stiffb и Stiffr. Они более приспособлены для таких случаев.
Выводы: Задача о движении двух гравитирующих тел имеет точное решение. Решая эту задачу численно, можно "на практике" проверить состоятельность различных методов численного интегрирования.
Замечание: Результаты численного моделирования задачи трех тел (более сложный случай изложенной выше задачи) изложены
тут.
, то с очень большой вероятностью они разлетятся, и никакой стабильной системы не образуют.
, 
) можно свести эту систему к системе обыкновенных дифференциальных уравнений первого порядка:
и
, которые притягиваются друг к другу с силой гравитационного притяжения 
.