Сегодня мы еще раз рассмотрим решение задачи 3-х тел. С точки зрения физики три тела в пространстве, разнесенные на некоторые расстояния и обладающие начальными скоростями, будут двигаться под действием Ньютоновских сил притяжения и с Ньютоновскими ускорениями. Третий закон Ньютона также участвует тем, что действие первого тела на второе сопровождается противоположно направленным действием второго на первое. Главными уравнениями являются:
Fij = G * Mi * Mj / R^2 (1)
где Fij – сила взаимодействия между телами i и j;
G – гравитационная постоянная;
Mi и Mj – массы тел i и j;
R – расстояние между телами i и j.
Следующее уравнение об ускорении W под действием силы F на тело массой M
Wi = -Fij / Mi; Ускорение «тормозит» тело i, поэтому минус.
Wj = Fij / Mj;
или
Wi = -G * Mj / R^2; (2)
Wj = G * Mi / R^2;
Сразу заметим, что силы, ускорения, скорости и координаты тел – вектора в трехмерном пространстве и должны вычисляться по правилам векторов.
С точки зрения математика тела движутся по своим траекториям так, что в каждой точке движения по траектории действует зависимость координат K от времени t
Kt = (W0*t^2)/2 + V0*t + K0 (3)
где K0 – координаты тела в момент времени t0;
V0 – вектор скорости тела в момент времени t0;
W0 – вектор ускорения тела в момент времени t0;
Здесь стоит отметить, что при сближении тел, каждое из них становится источником центростремительной силы, действующей на партнера так, чтобы он поворачивал свою траекторию вокруг некоторого мгновенного центра. Радиус кривизны и координаты центра мы постараемся определить из аналитической геометрии, математического анализа и других наук, наработанных человечеством за века.
Законы Ньютона утверждают, что силы взаимодействия можно определить для любого расположения тел в пространстве, а вот для координат тел и скоростей в любой момент времени законов нет. Поэтому ориентироваться будем на ускорения.
Рассмотрим тело, имеющее координаты K1 (вектор), движущееся со скоростью V (вектор) и с ускорением от действия совокупной силы притяжения других тел W (вектор)
. Рассматривать будем только проекцию на плоскость XY, так как рассуждения для всех плоскостей и осей одинаково. Векторная алгебра позволяет переносить начала координат векторов, поэтому совместим начала координат векторов ускорения и скорости в начало координат нашей диаграммы, изображенной на рис.1.
Здесь показано взаимодействие ускорения W с текущей скоростью тела. Чтобы понять, что же происходит с телом, разложим ускорение на две составляющие – Wd параллельно вектору скорости V и Wr перпендикулярно V. Дело в том, что вектор Wd просто ускоряет или тормозит тело, а вот вектор Wr является центростремительным ускорением, изменяющим траекторию движения тела. Именно из-за него формула (3) работает в этой программе очень плохо и требует предельного укорачивания отрезков интегрирования.
Из подобия треугольников
Wdx / Wdy = Vx / Vy (4)
(Wy – Wdy) / (Wdx-Wx) = Vx / Vy (6)
Отсюда
Wdx = Wdy * Vx / Vy
После подстановки
Wy – Wdy = (Wdy * Vx / Vy - Wx) * Vx / Vy
Wdy(1+Vx^2/Vy^2) = Wy + Wx * Vx / Vy
И наконец
Wdy = Vy * (Wy * Vy + Wx * Vx) / (Vy^2+ Vx^2)
Wdx = Vx * (Wy * Vy + Wx * Vx) / (Vy^2+ Vx^2)
Вектор Wr = ((Wx – Wdx); (Wy – Wdy)).
Модуль | Wr | = sqrt(Wrx^2 + Wry^2)
Уравнение центростремительного ускорения
Wc = V^2 / R
Где Wc – центростремительное ускорение;
V – тангенциальная скорость (по касательной к окружности) ; R – радиус окружности.
Скорость V нам известна, и будем считать, что она ничтожно изменится за время шага, ускорение Wc известно, можно определить радиус окружности, по которой полетит тело в предстоящий шаг интеграции.
R = (Vy^2+ Vx^2) / | Wr |;
Координаты центра вращения (Rx, Ry) можно определить исходя из координат тела, направления вектора Wr и длины радиуса, хотя использовать эти координаты нет необходимости.
Rx = Kx – R / | Wr | * Wrx;
Ry = Ky – R / | Wr | * Wry;
Тело полетит по дуге окружности радиуса R с центром (Rx, Ry), но очень быстро «слетит с рельсов». Ведь лететь оно будет с ускорением Wd, значит скорость возрастет/уменьшится и центр вращения сместится. Да и другие тела, создающие центростремительное ускорение, на месте не стоят. Поэтому ограничим дугу углом Pi/65536 То есть длина полета за один шаг будет
L = R / 65536;
А теперь определим время шага
L = W * T^2 / 2 + V * T = R / 65536;
T = (-V sqrt(V^2 – 4 * W * R / 65536)) / (2 * W);
Выполнение такого расчета наукообразно, но не дает никаких результатов. Длина ведь выбрана достаточно условно. Для упрощения расчетов можно определить время только по V, поэтому
T = L / V = R / 65536 / V;
Для прямой линии отрезок T будет бесконечным, для малых радиусов, где вдобавок огромные скорости, это будет очень малое время. Поэтому из трех периодов выберем наименьший и по нему будем выполнять шаг.
Посчитаем проекции на оси координат того пути, который пройдет тело с учетом начального и тангенциального ускорений и скоростей и обозначим их Dx и Dy.
Dx = (Wx + Wdx) * (T^2) / 4 + Vx * T
Dy = (Wy + Wdy) * (T^2) / 4 + Vy * T
Выберем большую из них, чтобы обеспечить максимальную точность и рассчитаем вторую координату на дуге по формуле
Dx = R – sqrt(R^2 – Dy^2), если Dy >=Dx
Dy = R – sqrt(R^2 – Dx^2), если Dx >=Dy
На рис.2 поясняется, откуда взялись эти формулы.
Естественно, что для трехмерного пространства нужно добавить в формулы определения векторов третье измерение, но в основном логика останется вышеизложенной.
Теперь по поводу того, как это лучше растолковать процессору, что он не только делал, что надо, но делал это без особого напряжения.
Мантисса числа с плавающей запятой имеет ограниченный размер, в большинстве случаев это 53 бита для формата DOUBLE и 64 бита для формата TBYTE. Это 16 или 19 десятичных знаков и не думаю, что координаты тел надо искать с большей точностью.
Чтобы минимизировать искажения при расчетах, связанных с шагом времени, а таковых большинство, следует интегрируемый интервал приравнивать к степени числа (1/2), тогда умножение любого числа на такое значение заменяется сложением порядка множимого со степенью интервала без изменения мантиссы. Это позволит достичь максимальной точности вычислений.
Описанный в статье алгоритм дает очень хорошие результаты для первого приближения искомой точки. Далее, необходимо сделать проверку на то, что рассчитанные координаты действительно соответствуют тем ускорениям, которые должны быть в следующей точке. При неточном попадании следует выполнить подгонку методом Эйлера или любым другим методом.
Для меня самым трудным был рисунок II.Ci.c.300 с периодом прохождения 198.6920179177 сек. Если другие версии программы проходили его за 40-60 минут, этот алгоритм сократил время до 2 минут. Особо придирчивым сообщу, что выход на начальные координаты чуть-чуть сдвинулся за счет «длинного» шага на последней итерации.
Все, желающие получить файл с программой для WINDOWS-64, могут скачать игрушку AccDbl36.exe здесь.
Спасибо, что дочитали до конца.