CC BY
54
АСТРОНОМИЯ
УДК 524.3/4-32 А. А. Башаков
ТЕСТИРОВАНИЕ ФАЗОВЫХ МОДЕЛЕЙ ЗВЕЗДНЫХ СИСТЕМ, ПОСТРОЕННЫХ МЕТОДОМ ШВАРЦШИЛЬДА*
1. Введение
На сегодняшний день имеется большое количество динамических моделей Галактики, в разной степени соответствующих реальности. Для их построения используется несколько основных методик: численные методы, теоретические алгоритмы, а также их комбинации. Каждый из указанных подходов имеет как сильные, так и слабые стороны, и универсального подхода на сегодняшний день не выработано. В данной работе описано использование несколько модифицированного метода Шварцшильда. К сожалению, найденная методом Шварцшильда фазовая модель, как правило, не является единственной, так как в основе метода лежит аппроксимация только пространственного распределения плотности, при этом кинематические характеристики зависят от сделанных допущений и могут сильно отличаться в различных моделях, построенных для одного и того же гравитационного потенциала. К тому же метод Шварцшильда не гарантирует, что найденное решение для фазовой модели устойчиво.
Надежда получить устойчивую модель основывается на том, что в методе Шварц-шильда используется стационарный потенциал, который задается или который однозначно определяется по заданному распределению плотности. В этом фиксированном гравитационном потенциале вычисляются траектории, и создается библиотека орбит, каждая из которых определяется (просчитывается) на достаточно большом интервале времени. Но в самом методе Шварцшильда имеется существенное огрубление результатов, связанное с разбиением пространства модели на конечное число ячеек. Такая дискретная модель может обеспечить только некоторое приближение исходного потенциала. Даже если бы в построенной модели удалось получить средние плотности в ячейках равными заданной плотности, то из-за дискретности представления пространства мы будем иметь только некоторую аппроксимацию потенциала, и, в общем случае, останутся определенные различия с моделью с непрерывным пространством. В реальности
* Работа выполнена при финансовой поддержке РФФИ (грант №08-02-00361а) и Совета по грантам Президента РФ для государственной поддержки молодых российских ученых и ведущих научных школ (грант №НШ-1323.2008.2).
© А.А.Башаков, 2008
и сами средние значения плотности в ячейках удается приблизить только с некоторой точностью, поскольку требуется аппроксимировать плотности в ячейках одновременно во всем объеме модели. К тому же довольно сложно точно учесть распределение плотности в разреженной периферийной части модели и в центральных областях, где при наличии большого градиента плотности также возникают трудности из-за огрубления на ячейках. Поэтому при построении фазовой модели пространственную плотность и гравитационный потенциал удается только приблизить к первоначально заданным.
Метод Шварцшильда не позволяет исследовать, насколько устойчива построенная модель к отклонениям от исходного гравитационного поля, к его флуктуациям, поскольку в методе построение модели проводится для фиксированного стационарного потенциала (распределения плотности). Но если имеется отклонение от заданного потенциала, то это влечет за собой изменение траекторий в построенной фазовой модели: меняется темп прохождения различных участков орбиты и меняется область пространства, заметаемая каждой конкретной звездой при движении. Это, в свою очередь, может привести к еще большему отклонению распределения плотности от исходного, а значит и изменению общего потенциала системы. Насколько чувствительна модель к отклонениям в распределении плотности или потенциала, является открытым вопросом и, следовательно, остается открытым вопрос об устойчивости построенной модели. Поэтому полученную модель системы нужно обязательно проверять на устойчивость, и эта проверка должна быть достаточно аккуратной и корректной.
В статье представлено краткое описание метода Шварцшильда, особенностей его программной реализации и процесса построения и исследования модели на устойчивость. Для проверки устойчивости ставится и решается следующая задача: по известным параметрам начальных данных и весов орбит получить набор материальных точек, которые при усреднении воспроизводили бы плотность и кинематические характеристики найденной методом Шварцшильда фазовой модели, и затем использовать методы задачи N тел, например, иерархический алгоритм, реализующий численное моделирование дискретной системы. Весь цикл алгоритмов испытан на сферическом потенциале Пламмера. Для него известно аналитическое представление фазовой модели, которая устойчива, и кинематические характеристики которой описываются сравнительно простыми уравнениями. Также в статье дается описание построения и исследования на устойчивость фазовой модели, построенной на основе потенциала Галактики Кутузо-ва—Осипкова. Полученные результаты показывают возможность использования данного метода для построения фазовых моделей, но, вместе с тем, требуется некоторая осторожность или наложение определенных ограничений в выборе начальных данных при реализации метода.
2. Предпосылки
В работе [1] методом Шварцшильда была построена фазовая модель двухкомпонентного потенциала Галактики, предложенного Кутузовым и Осипковым [2, 3]. Пространственная плотность воспроизводилась хорошо, а кинематические характеристики полученной модели существенно отличались от соответствующих наблюдательных параметров Галактики в окрестности Солнца [1]. Несоответствие реальных и полученных в модели кинематических параметров могло быть связанно с различными причинами: с одной стороны, кинематика Галактики хорошо исследована только в окрестности Солнца, а, с другой стороны, модель Кутузова—Осипкова гарантирует только соответствие хода наблюдательной плотности (круговой скорости), остальные кинематические параметры авторами модели не исследовались. Также существуют определенные огра-
ничения в методе Шварцшильда, связанные, в первую очередь, с неоднозначностью нахождения фазовой модели через набор траекторий.
Существует несколько аналитических моделей с известными фазовыми плотностями. Одна из простейших таких моделей — сферическая модель Пламмера [4]. Полит-ропная модель Пламмера была использована для проверки всех алгоритмов и программ, используемых для получения и исследования на устойчивость фазовой модели. На первом шаге методом Шварцшильда была построена фазовая модель, состоящая из совокупности непрерывных орбит. Затем она была по определенным правилам заменена набором материальных точек. На заключительном этапе для проверки ее устойчивости и сохранения исследуемых параметров была прослежена эволюция полученной дискретной системы как задачи N тел. В этом случае были получены удовлетворительные результаты как по построению собственно фазовой модели, так и по исследованию ее эволюции.
3. Метод Ш^варцшильда и его программная реализация
В основе метода построения фазовых моделей методом Шварцшильда лежит предположение о пропорциональности вклада звезд с рассматриваемыми типами орбит в общий гравитационный потенциал стационарной системы суммарному времени, в течение которого звезда на заданной траектории находится в данном участке пространства [5]. Метод построения фазовой модели можно представить в виде нескольких основных шагов.
• Разбиение пространства на заданное число ячеек конечного размера.
• Вычисление массы модели в этих ячейках в соответствии с исходным потенциалом с использованием уравнения Пуассона.
• Вычисление траекторий звезд и создание библиотеки орбит в заданном потенциале. При этом для каждой траектории и каждой ячейки запоминается суммарное время пребывания частицы и усредненные компоненты скорости, с которыми частица пересекает ячейку.
• Нахождение весовых коэффициентов орбит с учетом того, что масса, создаваемая звездой в ячейке, пропорциональна времени, в течение которого она находится в данной ячейке. При этом ставится задача получить наилучшее соответствие масс в ячейках, найденных из уравнения Пуассона, и масс, полученных суммированием орбит с найденными весовыми коэффициентами. Для нахождения весовых коэффициентов используется симплекс-метод.
• Вычисление кинематических характеристик модели на основе ранее сохраненных данных о скоростях звезд в ячейках с учетом найденных весовых коэффициентов.
В результате получаем модель в виде набора орбит, а также профили пространственной плотности и кинематические характеристики модели. Метод Шварцшильда не гарантирует, что найденная фазовая модель для заданного распределения плотности единственная. Найденные симплекс-методом наборы траекторий могут различаться между собой и давать отличающиеся кинематические характеристики модели при одном и том же пространственном распределении. Это дает возможность через варьирование исходных библиотек орбит и наложение дополнительных условий в симплекс-методе для получения нужной кинематики найти наиболее подходящие наборы орбит.
Однако остается открытым об устойчивости фазовой модели, заданной выбранным набором траекторий. Для проверки полученной модели предлагается заменить ее системой тяготеющих точек, при этом свойства дискретной модели должны с максимально возможной точностью воспроизводить свойства исходной непрерывной модели: ее пространственную плотность и кинематические характеристики.
Таким образом, тестирование полученной фазовой модели предполагает представление ее в виде набора «звезд» и, далее, рассмотрение эволюции полученной системы N тел. Выделим следующие два важных шага.
• Дискретизация —представление найденной фазовой модели из набора непрерывных траекторий с найденными «весами» в виде набора материальных точек.
• Исследование эволюции построенной дискретной модели методами задачи N тел.
3.1. Разбиение пространства модели на ячейки. Деление пространства при построении модели потенциала Кутузова—Осипкова и модели Пламмера было реализовано по-разному.
При построении фазовой модели Кутузова—Осипкова использовались ячейки в форме прямоугольных параллелепипедов. Размеры области пространства, которая разбивалась на ячейки, а, соответственно, и модели, брались порядка 20-30 кпк. В итоге получалось, что при количестве ячеек порядка 104 (20 х 30 х 20), размер ребра ячейки был порядка 1 кпк. Такая точность была недостаточна для моделирования центральных областей модели, и в то же время некоторая доля массы вообще не могла быть учтена, так как находилась за границами области, в которой проводился набор статистики. А вследствие того, что модель Кутузова—Осипкова представляет собой осесимметричную систему с существенным сжатием вдоль полярной оси, часть ячеек (до 40%) оказывалась вовсе незаполненной (траектории в них не заходили). Ограничение сверху на количество ячеек накладывается симплекс-методом, с помощью которого вычисляются веса орбит. Так, при числе ячеек порядка 10000 и количестве орбит порядка 3000, получалась матрица размером 10000 х 3000, и для стандартных вариантов решения симплекс-методом при вычислении весов накапливались существенные ошибки счета, и время вычисления было очень большим.
При исследовании модели Пламмера с учетом ее сферической симметрии метод разбиения пространства был модифицирован. Был произведен переход от прямоугольных к криволинейным ячейкам: концентрическим сферическим слоям при дроблении по радиусу. Также была предусмотрена возможность разбиения исследуемого пространства по азимутальному и зенитному углу. Существенной модификацией стала возможность не равномерного, а экспоненциального дробления на ячейки по радиусу, позволяющая постепенно уменьшать размеры ячеек по мере приближения к центру модели.
В рассмотренном варианте построения и исследования модели Пламмера разбиение по азимутальному и зенитному углу не требовалось вследствие изотропности распределения в слоях по координатам и скоростям. По радиусу было выбрано 1000 ячеек, с экспоненциальным коэффициентом равным 1.0 (постоянный шаг). При таких параметрах и размере исследуемой области в 10 единиц, радиус центральной ячейки был порядка 0.01, что обеспечивало достаточную точность для моделирования данной области.
3.2. Построение траекторий и начальные данные. На следующем шаге в заданном потенциале просчитывались базовые орбиты. Для получения траекторий звезд
использовался метод численного интегрирования уравнений движения, предложенный Эверхартом [6]. Вычисления проводились с двойной точностью.
Одной из проблем на этом шаге является выбор стартовых параметров для базовых орбит. Сложность состоит в том, чтобы максимально полно охватить все многообразие траекторий в модели, не нарушив при этом, в частности, изотропности распределения для модели Пламмера.
В случае потенциала Кутузова—Осипкова для вычисления начальных данных для запуска траекторий было использовано дробление по трем параметрам: по интегралу площадей J (определяет радиус круговой орбиты), интегралу энергии Е (ширина области допустимых движений) и радиусу К (в области, где движение возможно, с учетом первых двух параметров). При этом, когда производился переход к дискретной модели, траектории, заменяемые набором точек, воспроизводились однозначно. Минусом при таком подходе является то, что кинематические параметры получаемой модели существенно зависят от способа задания начальных параметров библиотеки орбит. Так, одним из недостатков конкретных моделей, полученных при таком задании начальных данных в потенциале Кутузова—Осипкова, была слишком большая дисперсия скоростей в направлении, перпендикулярном плоскости Галактики. Получалась «горячая» модель, которая, возможно, возникала из-за запуска траекторий для библиотеки орбит перпендикулярно к плоскости симметрии, то есть «вверх» в сопутствующей плоскости. Так на расстоянии орбиты Солнца дисперсия скоростей по оси Z в построенной модели получалась порядка 80 км/с и более при том, что наблюдаемые значения в Галактике оцениваются в 20-30 км/с.
При исследовании модели Пламмера алгоритм вычисления стартовых параметров был изменен. Дробление проходило по двум параметрам: радиусу К (от почти нулевого радиуса до границ модели) с неравномерным шагом и по интегралу энергии Е (верхняя граница определялась условием не покидания частицами заданных границ модели). Каждая пара этих параметров давала радиус запуска и модуль скорости частицы. Чтобы обеспечить изотропность модели, для каждой пары значений радиуса и модуля скорости запускалось несколько десятков траекторий (до 50), при этом координаты точки запуска выбирались случайным образом на сфере заданного радиуса, направление вектора скорости также выбиралось изотропно. При этом статистика (время пребывания в ячейке, скорости и дисперсии скоростей) вычислялась по всем орбитам, соответствующим исходной паре стартовых параметров. Таким образом, модель получалась более изотропна. Минусом такого подхода является то, что при переходе к дискретной модели невозможно в точности воспроизвести орбиты по паре стартовых параметров. Поскольку статистические характеристики фазовой модели в виде набора орбит, найденного симплекс-методом, должны совпадать с характеристиками модели, представленной в виде набора материальных точек, процедура дискретизации модели Пламмера оказалась несколько сложнее, чем для случая потенциала Кутузова—Осип-кова.
3.3. Набор статистики. Определение статистических характеристик происходило следующим образом: по фазовым координатам, полученным до и после вызова подпрограммы-интегратора, строился полином от времени, который использовался для нахождения моментов пересечения границ ячеек и, соответственно, интервалов времени, в течение которых частица находилась внутри ячеек. В случае потенциала Ку-тузова—Осипкова, где использовались ячейки сравнительно простой формы, строился полином четвертой степени на основе координат и скоростей до и после вызова
подпрограммы-интегратора (жп, «„, хп+1, «„+1). Для получения приемлемой точности необходимо было не более 10 точек на обороте траектории. При исследованиях модели Пламмера форма ячеек изменилась, в связи с чем схема построения полинома также была изменена. Для его построения использовались только пространственные координаты, и полином был второй степени. Для сохранения приемлемой точности определения времен пересечения границ ячеек шаг интегратора был уменьшен (использовалось ^50 построений полинома на один оборот траектории).
При построении модели Пламмера с учетом ее сферической симметрии были сделаны некоторые упрощения в программном модуле набора статистики. Так, было решено отказаться от разбиения ячеек по «долготе» и «широте», и для набора статистики сохранялись только две компоненты скорости. Радиальная скорость (вдоль вектора положения тела) может иметь как положительные, так и отрицательные значения. В случае если ее усредненное значение близко к нулю, можно говорить о стационарности модели (отсутствие систематического движения: расширения или сжатия модели). Вторая исследуемая компонента скорости — тангенциальная — вычислялась по полной и радиальной скорости, при таком способе вычисления она имеет только модуль и не имеет направления.
Для получения достоверных значений пространственной плотности и кинематических характеристик необходимо просчитать траекторию каждой «звезды» в течение продолжительного времени. Для контроля был применен следующий способ: при каждом посещении ячейки сравнивалось время, в течение которого частица находилась в ячейке, и уже накопленное за предыдущие посещения время. В случае если для большинства ячеек их отношение становилось меньше некоторого заданного числа (порядка нескольких сотых), вычисления прекращались.
3.4. Нахождение набора орбит симплекс-методом. Для нахождения выборки траекторий и их весовых коэффициентов, с помощью которых наилучшим образом аппроксимируется плотность, соответствующая заданному гравитационному потенциалу, использовался симплекс-метод, в котором решение системы из Ь линейных уравнений с К неизвестными ищется с учетом ряда ограничивающих условий [7, 8]. В нашем случае Ь равно числу пространственных ячеек, К — числу траекторий в библиотеке орбит. Налагались следующие условия на веса т* траекторий:
А — массив, содержащий суммарные времена а^, в течение которых г-я траектория находилась в ^’-й пространственной ячейке, Ш — столбец с искомыми весовыми коэффициентами орбит, Р — столбец с величинами, пропорциональными массам в ячейках, полученным из уравнения Пуассона △ и = — 4пОр.
3.5. Переход к дискретной модели (дискретизация). В результате работы комплекса программ, реализующих метод Шварцшильда, получается фазовая модель,
к
^ тах, АШ < Р, т > 0, г € {1, К},
где А
представленная в виде набора орбит с соответствующими весами. Для проверки равновесности и устойчивости производится переход к дискретной модели в виде набора «звезд», и проводится исследование эволюции полученной системы N тел.
Алгоритм вычисления фазовых координат тел в дискретной модели был следующим. Сначала, основываясь на ранее найденных весах орбит wi, вычислялось количество звезд на каждой орбите (количество пропорционально весу), при этом из-за необходимости округления количества звезд на орбите до целого общее число тел могло немного измениться (так при исследовании сферической модели Пламмера N изменилось со 100 000 до 99 743). Массы всех тел полагались равными. Далее, в случае потенциала Кутузова—Осипкова по трем характеризующим орбиту параметрам (J, E, R) вычислялись стартовые координаты и скорости. Затем время, в течение которого звезда на данной орбите совершает примерно от 50 до 100 и более оборотов, делилось на количество тел на данной орбите. В процессе интегрирования траектории делались «остановки» через полученные интервалы, при этом координаты и скорости, которые имело тело в тот момент, запоминались, и в дальнейшем использовались как начальные данные для задачи N тел. Подобная операция проводились для всех орбит с ненулевым весом. Сравнить параметры исходной и дискретной моделей можно по графикам на рис. 1. Исследуемые нами характеристики в дискретной модели получались близкими к характеристикам, подсчитанным в модели, состоящей из непрерывных орбит.
Для модели Пламмера алгоритм вычисления начальных данных для системы N тел был изменен, поскольку в данном случае нельзя для каждой пары (R, E) указать однозначно определенную орбиту. Процесс вычисления был следующим: для каждой пары параметров орбиты (R, E) вычислялись радиус запуска и модуль стартовой скорости — они определялись однозначно. Количество звезд на орбите полагалось пропорциональным весу орбиты. Для каждой звезды положение выбиралось случайным образом на поверхности сферы заданного радиуса, направление вектора скорости также выбиралось случайно. Далее, траектория просчитывалась в течение нескольких оборотов в заданном потенциале, после чего брались текущие координаты и скорости, которые и использовались впоследствии для задания начальных условий в задаче N тел. В данном алгоритме используется следующее предположение: так как для каждой пары стартовых параметров (R, E) запускается большое количество «звезд», можно считать, что они достаточно равномерно заполнят область пространства, которую занимали траектории, использовавшиеся в методе Шварцшильда, и воспроизведут их кинематические параметры. Просчет траектории в течение некоторого времени после запуска был необходим для «ухода» тел с поверхностей сфер. Результат дискретизации представлен на рис. 3.
В данном случае соответствие характеристик оказалось несколько хуже, чем в модели Кутузова—Осипкова, хотя, в общем, все особенности непрерывной модели присутствуют и в модели, представленной набором N тел (рис. 3).
3.6. Исследование эволюции системы. Для исследования эволюции дискретной модели использовался иерархический алгоритм численного решения задачи N тел. Применение этого алгоритма оправдано при N более 10000, алгоритм позволяет существенно сократить время счета, так как при прямом решении временные затраты O(N2), а в случае иерархического алгоритма O(N log N). В модуле исследования эволюции системы была использована готовая программная реализация алгоритма [9]. Все вычисления проводились с двойной точностью. Параметры алгоритма для потенциала Кутузова—Осипкова были взяты £=0.25 кпк, 0=1.0, для модели Пламмера £=0.025,
0=1.0. При исследовании потенциала Кутузова—Осипкова были рассмотрены два варианта: N =100 000 и N=200 000. Результаты эволюции оказались сходными, поэтому было принято решение при исследовании модели Пламмера использовать N =100 000 тел.
4. Рассмотренные потенциалы
Описанным методом были исследованы модели двухкомпонентного потенциала Галактики Кутузова—Осипкова [2, 3] и сферического потенциала Пламмера [4]. Потенциал Кутузова—Осипкова имеет вид
и(Д, г) = ио[а^>а(£) + Ьу>ь(г)], а + Ь =1
где ио = и(0,0) —потенциал в центре, а и Ь — веса компонентов, <^>а и — их безразмерные потенциалы. Компонента а, сплюснутая и более протяженная, моделирует распределение в диске и гало:
<Ра(0 = а « - 1 + Vх! + ;
компонента Ь с быстрым падением плотности представляет ядро и балдж:
у>ь(г) = (1 + Аг)-1.
Здесь
С2 = г2 + 2(1 — е) л/е2 + С2 — е , г = л/р2 + С2з
р = Д/го, С = г/го,
Д и г —цилиндрические координаты, го —единица длины.
Модель содержит два масштабных (го, ио) и пять структурных (а, к, е, А, а) параметров, определяющих ход круговой скорости, сплюснутость и веса компонент. При построении фазовой модели использовались следующие значения параметров [3], оценка которых была выполнена по наблюдательным данным Галактики: ио = 3, 81 • 105 км2/с2, го = 2, 54 кпк, а = 0, 37, а = 4,17, к = 0, 69, е = 0,11, А = 6, 60.
Потенциал модели Пламмера имеет вид
тт < \ СШр1 иР1(г)
В процессе построения фазовой модели были использованы следующие значения параметров: ОЫр1 = 1,Дрг = 0.6. Размеры модели были ограничены сферой радиусом 10 единиц.
5. Результаты построения и исследования эволюции фазовых моделей
5.1. Модель потенциала Кутузова—Осипкова. При построении фазовой модели на основе потенциала Кутузова—Осипкова была использована библиотека из 2400 орбит. Исследуемая область пространства с размерами 35 х 35 х 20 кпк была разбита на 12000 ячеек (30 х 20 х 20). Симплекс-методом было отобрано 1008 орбит (42%), при том, что модель заняла примерно 7200 ячеек (60%). Масса модели, полученная суммированием орбит с найденными весами, отличалась от теоретической массы менее чем на
5%. О соответствии распределения массы в пределах модели можно судить по графику на рис. 1. Для кинематических параметров такое сравнение произвести нельзя, так как авторами модели ее кинематические свойства не исследовались.
Дискретизация найденной фазовой модели проводилась для N =100 000. Сравнение графиков на рис. 1 показывает хорошее согласие характеристик в дискретной и непрерывной моделях.
10000
1000
100
10
1
0 5 10 15 20 25 30 35 40
К (крс)
Рис. 1. Графики хода плотности в потенциале Кутузова—Осипкова.
Проверка на устойчивость показала, что полученная дискретная модель в результате исследования ее как системы N тел довольно быстро эволюционирует. Для оценок в качестве характерного интервала времени использовалось время оборота точки на круговой орбите, соответствующей Солнцу (галактический год). После пяти галактических лет распределение пространственной плотности в целом сохраняется, за исключением некоторого «размазывания» границ модели и сближения между собой значений дисперсий скоростей по всем трем компонентам. До пятнадцати галактических лет эволюция носит тот же характер, доля ушедших частиц около 1%, и ее можно объяснить тесными кратными взаимодействиями при решении задачи N тел. Далее эволюция замедляется, модель стабилизируется, и графики для 20 галактических лет уже практически повторяют соответствующие графики, построенные для 15 галактических лет.
5.2. Модель сферы Пламмера. Для построения фазовой модели Пламмера были использованы следующие значения параметров: размеры модели были ограничены 10 единицами, при наборе библиотеки орбит для вычисления начальных данных траекторий использовалось 100 значений радиусов сфер, с которых происходил запуск, для каждой сферы вычислялись 10 значений энергии Е, и для каждой пары Д, Е запус-
К (крс)
Рис. 2. Азимутальная скорость в потенциале Кутузова—Осипкова.
калось по 40 траекторий. Исследуемая область пространства была разбита на 1000 сферических ячеек. Количество тел, использующееся для представления дискретной модели, было равно 100 000.
Такие начальные данные давали 1000 орбит на 1000 ячеек, симплекс-методом было отобрано 502 орбиты (50%), на основе которых и была построена фазовая модель. Для сравнения на рис. 3 приводятся профили плотности, скоростей и дисперсий скоростей на различных этапах построения модели.
Так как кинематические параметры модели, сохраняемые в модуле сбора статистики, не были полными и предполагали потерю части информации, прямое сравнение с аналитической моделью сферы Пламмера было затруднено. Поэтому более удобным в качестве «стандарта» было бы использование некоторого дискретного, протестированного представления аналитической модели Пламмера, построенного способом, отличным от метода Шварцшильда.
В качестве «стандартного» варианта была взята дискретная модель сферы Пламмера, генерируемая в примененном иерархическом алгоритме. В блоке, реализующем иерархический алгоритм, его авторами был встроен механизм тестирования, основанный на построении сферы Пламмера и расчете ее эволюции [10]. Для проверки модуля набора статистики данная тестовая модель была просчитана в блоке, исследующем эволюцию системы, при этом вычислялись все исследуемые нами характеристики (профили плотности, скорости и дисперсии скоростей).
Из анализа результатов, полученных при построении и исследовании эволюции модели, построенной нами на основе потенциала сферы Пламмера, можно сделать следующие выводы.
01 23456789 10
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 ■0.1
01 23456789 10
vR — V0 —
- \ ое —
\
--— _
■ ;■
! \
01 23456789 10
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 ■0.1
0 1 2 3 4 5 6 7 8 9 10
vR v0
Ц
.. \ ое —
чу
і ‘"•'-Г' !
■ ■
100
80
60
40
20
0
obtained
\ thee reti за!
\
\
4
0 1 2 3 4 5 6 7 8 9 10
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
vR
V0
S0
і ---ч ■\
чу
ч>
-;х.
0 1 2 3 4 5 6 7 8 9 10
Рис. 3. Модель Пламмера: профили плотности (слева) и скоростей и дисперсий скоростей (справа). Сверху вниз: «стандартный вариант»; модель, полученная в результате работы симплекс—метода (суперпозиция орбит с найденными весами); дискретная модель (начальные условия для задачи N тел); после 20 средних времен пересечения.
Метод Шварцшильда дает хорошее приближение распределения пространственной плотности, но кинематические характеристики очень сильно зависят от библиотеки орбит, использовавшейся при построении модели. С учетом некоторых априорных знаний о распределении скоростей удалось сформировать подходящую библиотеку орбит и по ней построить фазовую модель. Кинематические характеристики этой модели в целом близки к «стандартному варианту» внутри сферы радиусом 5 единиц (Кр\ = 0.6), на
большем удалении характеристики уже существенно различаются. Причина этого в недостаточном количестве орбит, лежащих в этой области и использующихся при построении модели. Масса в модели Пламмера, находящаяся за пределами сферы радиусом 5 единиц, мала (~2%), и поэтому вполне допустимо примириться с несоответствием кинематики в этой части модели, по сравнению со «стандартным» вариантом.
В результате применения предложенного способа дискретизации была получена модель в виде набора материальных точек. Графики на рис. 3 дают возможность оценить степень ее соответствия исходной модели. В данном случае имеются отклонения как в профиле распределения плотности, так и в профилях кинематических характеристик, хотя в целом можно сказать, что полученная дискретная модель сохраняет основные свойства исходной непрерывной модели. Причина, по которой имеются некоторые несоответствия характеристик, в том, что при использованном нами способе дискретизации орбиты воспроизводились случайным образом с сохранением только интеграла площадей и энергии, при этом сами орбиты на этапе сбора статистики для симплекс-метода могли отличаться от соответствующих орбит на этапе дискретизации. Для уменьшения расхождений параметров дискретной и непрерывной моделей следует увеличивать количество запускаемых траекторий на каждую пару параметров (Д, Е). В данной реализации запускалось по 40 траекторий. Возможно в дальнейшем этот параметр будет увеличен, но на данном этапе полученная точность была приемлема.
По полученным результатам эволюции системы N тел можно заключить, что построенная дискретная модель устойчива, по крайней мере, на протяжении прослеженных 50 средних времен пересечения.
6. Выводы и заключение
Из анализа полученных результатов можно заключить, что в целом комплекс программ на основе метода Шварцшильда пригоден для построения фазовых моделей, хотя и требует проверки полученной модели на равновесность и устойчивость, так как алгоритм чувствителен к используемой библиотеке орбит. Для тестирования устойчивости предлагается алгоритм дискретизации фазовой модели и последующего исследования дискретной модели как задачи N тел.
Метод был применен для построения фазовых моделей на основе потенциала Галактики Кутузова—Осипкова и сферического потенциала Пламмера. Для каждого потенциала была построена фазовая модель (в виде суперпозиции орбит) и, затем, на ее основе дискретная модель, которая проверялась на устойчивость (задача N тел). Полученные дискретные системы эволюционировали, но в целом сохраняли свои характеристики.
Для того чтобы избавиться от сильной зависимости кинематических характеристик полученной фазовой модели от библиотеки орбит, в дальнейшем планируется усовершенствовать модуль вычисления весов орбит с целью аппроксимации не только пространственной плотности, но и кинематических характеристик. Это позволит упростить выбор начальных параметров траекторий для библиотеки орбит и в будущем сделать проще применение метода на новых исследуемых потенциалах.
В заключении выражаю благодарность Николаю Петровичу Питьеву и Виктору Владимировичу Орлову за полезные консультации и помощь в работе.
1. Башаков А. А., Питьев Н. П. Построение самосогласованных моделей звездных систем методом Шварцшильда // Вестник СПбГУ, 2007. Сер. 1. Вып. 3. С. 151-159.
2. Кутузов С. А., Осипков Л. П. Двухкомпонентная модель гравитационного поля Галактики // Астрон. журн., 1989. Т. 66. Вып. 5. С. 965-973.
3. Кутузов С. А., Осипков Л. П. Оценка параметров двухкомпонентной модели Галактики интервальным методом // Вопросы небесной механики и звездной динамики. Алма-Ата, 1990. С. 110-116.
4. Binney J., Tremaine S. Galactic Dynamics // Princeton University Press, 1987. P. 223.
5. Schwarzschild M. A numerical model for a triaxial stellar system in dynamical equilibrium // Astrophys. J., 1979. Vol. 232. P. 236-247.
6. Everhart E. Implicit single-sequence methods for integrating orbits // Celest. Mech., 1974. Vol. 10. P. 35-55.
7. Ашманов С. А. Линейное программирование. М., 1981.
8. Васильев Ф. П., Иваницкий А. Ю. Линейное программирование. М.: Факториал Пресс, 2003.
9. http://ifa.hawaii.edu/ barnes/treecode/treeguide.html
10. Aarseth S. J., Henon M., Wielen, R. A comparison of numerical methods for the study of star cluster dynamics // Astronomy and Astrophysics, 1974. Vol. 37. P. 183-187.
Статья поступила в редакцию 18 мая 2007 г.
