Как говорится, новое – это хорошо забытое старое. Причем чем старее, тем лучше. Интегрирование траектории трех притягивающихся тел в идеальном изолированном пространстве – задача в аналитическом виде не решаемая. Зато в численной форме решений много. Предлагаю алгоритм очень быстрого интегрирования без применения высокоточных переменных и каких-то расширенных форматов. Кроме того, вместо многоминутных расчетов, этот алгоритм затрачивает на все не более 15 секунд.
Суть алгоритма заключается в использовании подстановки кривой третьего порядка в качестве траектории на интегрируемом участке. Естественно, кривая строится в трех координатах по шагам. На каждом шаге идет контроль за изменением векторов скоростей и, в случае значительных отклонений, интервал шага делится на два, после чего шаг пересчитывается. После удачного шага размер интервала восстанавливается.
Итак, подробно. Будем считать, что координаты следующей точки отвечают выражениям:
X1 = Ax * t^3 + Bx * t^2 + Cx * t + Dx
Y1 = Ay * t^3 + By * t^2 + Cy * t + Dy (1)
Z1 = Az * t^3 + Bz * t^2 + Cz * t + Dz
Чтобы в дальнейшем не загромождать текст, буду писать все рассуждения только для Х, для остальных координат они такие же. Скорость тела – производная координаты по времени
V = dX/dt = 3*Ax*t^2 + 2*Bx*t +Cx (2)
Ну и ускорение обозначим через W
W = dV/dt = 6*Ax*t +2*Bx (3)
Если t=0, то Dx=X0; Cx = V0; Bx = W0/2.
Координаты X0 и V0 заданы в начальных условиях, W0 можно просчитать для момента T=0. Для расчета таинственного Ax (не забудьте Ay, Az) надо перейти в следующую точку. Для этого нам задан шаг интегрирования, например 0,00001 сек и примем в начальный момент, что в начале итерации Ax = 0. Подставим данные в формулу координат, получим первое приближение следующей точки. Для момента времени T1 формула вектора ускорений будет:
W1 = 3*Ax*T1 + W0 (4)
Рассчитаем фактический вектор ускорений в этой точке, исходя из сил притяжения. Полученное значение позволяет определить неизвестный коэффициент:
Ax = (W1 – W0)/(3*T1) (5)
Зная Ax (а также Ay и Az), можно уточнить координаты новой точки для T1. Проделав это уточнение пару раз, получим истинное положение новой точки. Определить составляющие скорости в точке 1 также несложно:
V1 = 3*Ax*T1^2 + W0*T1 + V0 (6)
Перенесем координаты, скорости и ускорения новой точки вместо исходной и повторим расчет столько раз, сколько требуется по заданию. Все было бы хорошо, но дьявол кроется в мелочах. Если интервал велик, можно пропустить момент «близкой встречи» тел или момент остановки тела и разворота движения обратно. Тогда тела просто разлетятся и вряд ли вернутся в расчетную область. Чтобы этого не произошло, необходимо контролировать изменение вектора скорости каждого тела. Неважно, поворачивает ли оно, тормозится или ускоряется, такие моменты надо обсчитывать более тонко, то есть с меньшим шагом T, например делением его пополам.
Для пояснения математики вернемся к трем координатам. Модуль разности скоростей |V1-V0| равен
Rv = sqrt((V1x-V0x)^2+(V1y-V0y)^2+(V1z-V0z)^2) (7)
Модуль скорости |V0|= sqrt(V0x^2+V0y^2+V0z^2) (8)
Если Rv / |V0| > sigma, то необходимо делить интервал. Для уменьшения затрат процессора, я использовал соотношение
Rv^2 / (|V0|)^2 < 1/32768
Кроме этот в программе имеется счетчик числа делений интервала, чтобы не произошло зацикливание при случайных залетах. Считается, если лимит делений исчерпан, значит получена достаточная точность.
Еще раз подчеркну, программа написана на самых обычных переменных формата DOUBLE и может быть воспроизведена на любом современном языке программирования. Я писал на любимом мною Ассемблере.
Большое спасибо за прочтение статьи!