Поставив перед собой задачу компьютерного моделирования движения n тел с учетом их взаимного гравитационного воздействия для исследования динамики Солнечной системы, функционирования "черных дыр" и развития Вселенной, автор столкнулся множеством проблем.
В ходе работы пришлось последовательно шаг за шагом решить следующие задачи, каждая из которых по своему уникальна:
1. Разработана новая универсальная высокоточная переменная "BigDouble" повышенной разрядности и выделены особенности алгоритмов на ее основе;
2. Исследованы проблемы точности, впервые показано влияние разрядности на погрешность решения системы дифференциальных уравнений движения n тел;
3. Разработана модель гравитационного поля массивного кольца (тора) произвольных размеров и массы;
4. Разработана универсальная модель гравитационного поля Земли в любой точке с учетом формы эллипсоида вращения, с учетом заданного распределения плотности вещества планеты по объему и влияния других факторов в трехмерном приближении;
5. Разработана эскизная модель и анимация процесса циклического сжатия Вселенной до критического состояния и "Большого взрыва".
По сути это план работы и порядок публикаций в ДЗЕН автора по данной теме. Статьи не строго научные, поэтому в них максимально просто описаны результаты исследований по пунктам нашего плана. Постараемся не утомлять наших читателей сложными математическими формулами и выкладками, хотя математически все строго выверено и соответствует законам классической механики.
В предыдущей статье мы показали, что главное преимущество действительной переменной "BigDouble" состоит в высокой точности представления чисел в программе, а также в том, что разрядность мантиссы и предельную величину экспоненты можно задавать в широких пределах.
Проверку точности интегрирования дифференциальных уравнений движения системы n-тел мы провели на тест-задаче, в которой 12 тел с массой Солнца обращаются по кругу радиуса R = 6 астрономических единиц (примерно радиус орбиты Юпитера). В общем случае предусмотрено решение задач моделирования движения тел в трех измерениях, но в тест-задаче использованы только 2 координаты X и Y и, соответственно, скорости Vx и Vy.
Все тела расположены симметрично с угловым интервалом 30 градусов. Это позволило задать начальные координаты Xo,Yo и составляющие скорости тел Vxo и Vyo с требуемой точностью. Важное значение имеет то, что начальная скорость Vo каждого тела по касательной задана из условия равенства центробежной силы, действующей на каждое тело суммарной силе гравитационного притяжения по закону И. Ньютона.
Эта симметричная тест-задача для тел одинаковой массы весьма неустойчива и гораздо сложнее "задачи трех тел", но хороша тем, что, хотя интегрирование ведется по координатам X и Y, все тела движутся строго по окружности постоянного радиуса R=sqrt(X^2+y^2) и постоянного модуля скорости V=sqrt(Vx^2+Vy^2). Даже малейшее отклонение от точного расчета приводит к быстрому накоплению погрешности R и V. Относительную погрешность dV/V или dR/R можно легко определить на любом шаге интегрирования.
Надо заметить, что накопленная относительная погрешность модуля скорости V имеет тот же порядок величины, что и R, а в случае, когда шаг интегрирования по времени оптимальный, то значения dV/V и dR/R могут быть очень близкими по величине и противоположными по знаку. Это наблюдение возможно предстоит более тщательно исследовать с целью дальнейшего повышения точности расчетов.
На совмещенном графике рис.1 траекторий тел по кругу и скоростей показано изменение скорости Vx (синий цвет) и скорости Vy (красный цвет) одного выбранного тела по времени, модуль вектора скорости V показан зеленым цветом. Из графика видно, что составляющие вектора скорости Vx и Vy изменяются строго по гармоническому закону. Ничего другого для равномерного кругового движения и быть не может.
Полученный результат численного интегрирования уравнений движения тестовой симметричной периодической задачи 12 тел при разных значениях разрядности переменной показан на следующей картинке рис.2 .
Дальнейшее снижение погрешности ограничено применяемым методом интегрирования систем дифференциальных уравнений. В данном примере использован достаточно устойчивый и точный метод Рунге-Кутты. Однако, как известно, нет предела совершенству и над повышением устойчивости и точности метода интегрирования можно и нужно еще поработать.
Анализируя график рис. 2 можно сделать вывод, что применение обычных переменных Double с разрядностью 16 десятичных разрядов или даже Long Double с разрядностью 20 десятичных разрядов не обеспечивает точности интегрирования уравнений движения с шагом по времени от слова совсем.
Более того, исследования показывают, что эти переменные обладают отрицательным свойством быстро накапливать погрешности округления при использовании в циклах с большим количеством повторений. А в задачах интегрирования уравнений движения как раз могут применяться циклы с миллионами шагов по времени. У новых переменных BigDouble этот эффект полностью, даже в последних разрядах отсутствует, проверено до значений 100 млн. циклов.
Для интересующихся проблемой округлений в комментариях можно дополнительно дать небольшой фрагмент программы в кодах на языке С++, где эффект катастрофического накопления погрешностей округления обычных переменных хорошо виден.
В заключение хочу привести пример циклической задачи 12 тел, у которой тангенциальная скорость меньше равновесной. При такой скорости тела движутся по эллипсам и графики составляющих векторов скоростей Vx и Vy , изображенные на рис 3 уже отличаются от синусоиды, а также друг от друга.
.Чтобы не перегружать статью, автор не стал приводить примеры точности вычислений интегральных характеристик системы n-тел. При достаточно точном решении уравнений движения все законы сохранения выполняются автоматически, в том числе: закон сохранения энергии, закон сохранения импульса и момента импульса, а также теорема о движении центра масс системы.