Содержание к диссертации
Введение
2 Обзор 12
2.1 О фазовых моделях звездных систем 12
2.2 Оригинальный метод Шварцшильда 15
2.3 Методы, оперирующие частицами 18
2.4 Итерационный метод Родионова и Сотниковой 25
2.5 Моделирование гравитационного потенциала звездной системы . 26
3 Метод 28
3.1 Основные идеи метода Шварцшильда 28
3.2 Основные шаги построения модели 30
3.3 Разделение пространства модели на ячейки 32
3.4 Построение библиотеки орбит 35
3.5 Вычисление весов орбит 39
3.5.1 Симплекс—метод 40
3.5.2 Итерационный М2М алгоритм 42
3.6 Тестирование построенных моделей 45
3.6.1 Переход к дискретной модели N тел 45
3.6.2 Численное решение задачи N тел 47
3.7 Зависимость результатов от используемой библиотеки орбит . 48
4 Исследуемые потенциалы 51
4.1 Потенциал Галактики в модели Кутузова—Осипкова Т~~5Г~
4.2 Сферический потенциал Пламмера 52
4.3 Потенциал Галактики в модели Флинна и др 54
5 Результаты 56
5.1 Построение модели в потенциале Галактики Кутузова—Осипкова 56
5.2 Модель со сферическим потенциалом Пламмера 63
5.3 Построение модели в потенциале Флинна и др 70
5.3.1 Применение кинематических параметров 70
5.3.2 Обсуждение 76
5.4 Выводы 82
Заключение 85
Литература 91
- Оригинальный метод Шварцшильда
- Моделирование гравитационного потенциала звездной системы
- Построение библиотеки орбит
- Численное решение задачи N тел
Введение к работе
1.1 Актуальность темы
Проблемы строения, формирования и эволюции нашей и других галактик являются одними из самых актуальных в современной звездной астрономии. Одним из способов, позволяющим добиться понимания динамических процессов, происходящих в галактиках, является построение фазовых моделей звездных систем и численное моделирование системы частиц в рамках задачи N тел. Распределение пространственной плотности получаемых систем должно соответствовать наблюдательным данным или в общем случае наперед заданному закону, при этом обычно рассматриваются модели, близкие к равновесию. Одним из способов построения фазовых моделей является метод, предложенный в 1979 году Мартином Шварцшильдом. Метод относится к численным, и на ранних этапах использовался не очень широко из-за отсутствия достаточных вычислительных мощностей. Однако с развитием компьютеров эта проблема отошла на второй план, а возможность хорошо аппроксимировать заданное распределение плотности и некоторые другие параметры делают метод Шварцшильда привлекательным для решения задач по построению фазовых моделей звездных систем.
Методы численного моделирования и исследования системы частиц в рамках задачи N тел также важны. Часто используют термин "численный эксперимент", когда подразумевается численное моделирование и прослеживание эволюции системы гравитирующих N тел. В настоящей работе эти методы применяются для тестирования и исследования построенных модифицированным методом Шварцшильда моделей, так как при численном построении всегда имеются некоторые отклонения параметров в получающемся распределении плотности и соответственно от заданного гравитационного потенциала. Эти отклонения, в свою
очередь, могут привести к еще большему перераспределению плотности и уходу от заданного потенциала. Поэтому найденную модель необходимо проверить, насколько она получилась равновесной и устойчивой и сохраняет свои параметры.
В данной работе решаются вопросы реализации метода Шварцшиль-да на современном уровне и его модификации для расширения области применения. Разработаны и реализованы алгоритмы перехода от моделей, состоящих из набора фазовых траекторий, к дискретным моделям, представленным системой N точечных масс. Также дается описание построенных с помощью созданного программного комплекса фазовых моделей звездных систем.
1.2 Цели работы
В данной работе были поставлены следующие цели.
Модификация метода Шварцшильда для построения фазовых моделей звездных систем.
Разработка и создание программного комплекса для построения фазовых моделей звездных систем с применением алгоритмов модифицированного метода Шварцшильда.
Разработка метода перехода от моделей, состоящих из набора фазовых траекторий, к дискретным моделям, представленным системой N точечных масс, и получения начальных данных для моделирования задачи N тел.
Апробация разработанных алгоритмов и программ на модельной задаче построения фазовой модели со сферическим потенциалом Пламмера [13].
Построение фазовых моделей звездных систем на основе двухком-понентного потенциала Кутузова-Осипкова [9, 10] и трехкомпо-нентного потенциала Флинна и др [15] для Галактики.
Проверка полученных моделей на равновесность и устойчивость.
1.3 Научная новизна
В данной работе впервые применен метод Шварцшильда [16] с новым алгоритмом вычисления весов орбит, основанным на made-to-measure (М2М) алгоритме [1, 8]. Данная модификация позволяет учитывать кинематические характеристики при построении моделей. Добавлена возможность применения ячеек с переменными размерами и разной геометрической формы. С использованием модифицированного метода Шварцшильда были построены фазовые модели на основе потенциала Флинна и др. [15] для Галактики.
Разработан алгоритм дискретизации, позволяющий получать дискретные системы N тел (материальных частиц) из моделей, основанных на фазовых траекториях. Алгоритм используется для исследования построенных моделей на равновесность и устойчивость.
Впервые этим методом построены фазовые модели осесимметричных звездных систем на основе потенциалов Кутузова-Осипкова [6] и Флинна и др. [2] для Галактики. Исследована их равновесность и устойчивость. Показано, что разработанный метод пригоден для решения поставленных задач.
1.4 Научная и практическая значимость работы
В результате проделанной работы разработан и апробирован новый комплекс алгоритмов и программ, предназначенных для построения фазовых моделей звездных систем по заданному потенциалу. Данный комплекс основан на синтезе метода Шварцшильда [16] и М2М алгоритма [1, 8]. Новый модифицированный метод Шварцшильда позволяет строить фазовые модели звездных систем по заданному потенциалу и/или распределению плотности. Преимущества нового метода заключаются в возможности построения фазовых моделей на основе только распреде-
ления вещества без дополнительных предположений о наличии интегралов движения. С введением улучшенного метода расчета весов появилась возможность построения моделей с учетом заданных параметров скоростей. Результатом работы метода является модель в виде набора фазовых траекторий и в виде набора фазовых координат точечных масс, которые могут быть использованы как входные данные для численного решения задачи N тел.
1.5 Апробация работы
Результаты, изложенные в диссертации, докладывались на семинарах Кафедры небесной механики СПбГУ и лаборатории небесной механики и звездной динамики Математико—механического факультета СПбГУ, 2006-2009 гг.; на общегородском семинаре по звездной динамике им. К.Ф. Огородникова 2006-2009 гг.; на Всероссийской конференции "Звездные системы" к 100-летию П.П.Паренаго, Москва, 24-26 мая 2006 г.; на конференции "ВАК-2007", Казань, 17-22 сентября 2007 г.; на международной конференции "Dynamics of galaxies", Пулково, 2007 г.; на семинаре обсерватории Университета г. Турку, Финляндия, 2008 г.
1.6 Структура и объем диссертации
Оригинальный метод Шварцшильда
Мартин Шварцшильд, описав свой метод в статье 1979 года [49], стал пионером данного направления. Его метод основан на нахождении набора и определения весов звездных орбит, распределение которых воспроизводит заданный гравитационный потенциал, в котором они вычислялись. Использование этого метода не зависит от числа изолирующих интегралов движения, существующих для данной системы или заданного гравитационного потенциала. При этом предполагается, что система находится в стационарном или квазистационарном состоянии. В работе [49] строились трехосные модели эллиптических галактик. Метод позволял использовать как осесимметричные, так и трехосные потенциалы, так как наблюдательные данные и некоторые численные эксперименты указывали на то, что для части галактик могут лучше подходить модели с потенциалом на основе трехосного эллипсоида. Задача, стоявшая перед М. Шварцшильдом, состояла в разработке методики, позволяющей строить трехмерные самосогласованные модели с заданным распределением плотности. Описанный Шварцшильдом метод [49] не предполагает проверки на устойчивость, так как он был разработан только для получения таких моделей. Считается, что общий гравитационный потенциал стационарен, и все вычисления траекторий, требуемые в методе, проводятся в фиксированном потенциале. Здесь следует отметить, что это, на самом деле, является серьезным допущением, так как изначально неизвестно, возможно ли при взятом потенциале или распределении плотности устойчивое решение для фазовой модели.
Даже если распределение плотности основывается на реальных наблюдательных данных некоторой галактики, остается неизвестным, не имеется ли у этой галактики невидимая, скрытая масса, которая осталась не учтенной в распределении плотности. Зато метод, предложенный М. Шварцшильдом, имел серьезное преимущество перед аналитическими методами, позволяя строить модели для потенциалов сложной формы без рассмотрения вопроса о существовании дополнительного третьего интеграла. Отметим основные шаги при построении, рассмотренные в статье [49]. а) Выбор функции распределения плотности, в соответствии с кото рой предполагается построить модель. Принятое распределение играет роль входных данных метода, задается один раз в начале и в дальнейшем не меня ется. Данные предполагается получать, например, из распределения звездной плотности для галактик. Распределение представляет собой набор значений масс в ограниченных ячейках, которые могут быть как стационарны в про странстве, так и, например, вращаться как цельная система вокруг выбранной оси. Выбрав количество (Л се ) ячеек для построения и разбив пространство соответствующим образом, получим распределение плотности в виде вектора D(J), где J Є (l..NcdL). б) Получение гравитационного потенциала. Из данного распределе ния плотности потенциал вычисляется путем решения уравнения Пуассона. Здесь М. Шварцшильд [49] указывает на возможные трудности, связанные со сложностью получения алгебраического вида для потенциала. Если получение аналитического представления потенциала затруднительно или невозможно, то, как альтернативу, можно использовать табличные значения и интерполя цию для вычисления сил, действующих на частицы, однако это затрудняет программирование и ведет к росту времени счета. в)
Вычисление орбит. Для данного потенциала и с учетом особенностей распределения плотности находится набор траекторий. Каждой орбите ста вятся в соответствие начальные координаты и скорости, с которыми она вы числяется в рассматриваемом потенциале. Эти данные определяются один раз на этапе подготовки. Затем орбита "запускается" в заданном потенциале, и прослеживается на интервале порядка 100 времен пересечения системы для надежного представления орбиты как элемента галактики на длительном протяжении ее существования. В процессе вычисления каждой траектории запоминаются времена пребывания в ячейках, упомянутых в шаге а). В случае многократного пересечения орбитой одной и той же ячейки, времена суммируются. Итогом вычисления всей совокупности траекторий является набор (двумерная матрица) суммированных промежутков времени B(I, J), где / — номер орбиты, J — номер ячейки. Величины вдоль 1-й строки матрицы можно рассматривать как отражающие распределение плотности по ячейкам, соответствующее 7-й орбите. г) Контроль различия/сходства траекторий. Это требуется для уменьшения времени счета путем исключения орбит с "одинаковыми" функциями распределения плотности. Так, например, если у траектории нет дополнительного интеграла, то она заполнит все ячейки в области пространства, ограниченной интегралом энергии, точно также, как и другие орбиты, не имеющие дополнительного интеграла и с тем же значением интеграла энергии. Такие орбиты находятся из сравнения векторов B(I, J) и их можно объединить для уменьшения количества переменных на следующем шаге. е) Получение модели с заданным распределением плотности. После построения L орбит возникает вопрос, можно ли, запустив на каждую орбиту определенное количество "звезд", равное С(1), воспроизвести распределение плотности, заданное на шаге а)? Этот вопрос можно сформулировать через математические соотношения: имеет ли решение система с учетом ограничений на С(1): где D( J) - плотность в J-Pi ячейке, L — число траекторий. Условие неотрицательности С(1) связано с тем, что количество звезд или вес орбиты не должны быть отрицательными. Решение таких систем является типичной задачей линейного программирования.
Возникает вопрос о существовании решения. Если ответ "нет", то это может быть вызвано, например, неполным покрытием модели вычисленными L орбитами. В таком случае требуется добавление дополнительных орбит, которые увеличат область "занятого" объема пространства. Если же решение невозможно найти, несмотря на большое количество орбит, покрывающих предполагаемую область модели, значит, можно сделать предварительный вывод о том, что выбранное распределение плотности не может быть получено таким способом. Если решение удается получить, то выбранное на шаге а) распределение плотности может соответствовать конфигурации, находящейся в динамическом равновесии и может быть представлено в виде суперпозиции звезд на орбитах, в соответствии с результатами решения системы. Вопрос о единственности полученного решения оставлен автором [49] открытым. Описанный метод автор применяет к построению фазовой модели с плотностью, представляющей из себя суперпозицию сферической компоненты и трехосной эллиптической компоненты, которые моделировали хаббловский профиль плотности у эллиптических галактик, нужную сплюснутость и бар. Первые попытки решения системы уравнений не увенчались успехом, и поэтому был проведен анализ фазовых плотностей орбит. Выяснилось, что часть ячеек оказалась не заполнена, а некоторые орбиты имели похожие распределения. Набор орбит подвергся корректировке, часть орбит была исключена и добавлено некоторое количество новых орбит для заполнения пустых ячеек. Из них и была построена модель с необходимым распределением плотности. Подробного исследования свойств полученной модели не проводилось, однако сам факт возможности построения модели с необходимым распределением
Моделирование гравитационного потенциала звездной системы
Распределение масс в дисковых галактиках, в первую очередь, основано на анализе кривых вращения. Обзор определений масс звездных дисков (сплюснутых компонентов галактик) по кинематическим данным был сделан в работе [55]. Кривая вращения галактики может быть представлена графиком, на котором отображена зависимость скорости вращения центроидов звезд в галактике от расстояния до центра галактики. Данные наблюдений показывают, что в большинстве случаев галактики имеют плоские кривые вращения на периферии, то есть центроиды вращаются вокруг центра галактики с почти постоянной линейной скоростью в широком диапазоне расстояний от центра галак тики. Это означает, что звезды вращаются в периферийной области гораздо быстрее, чем ожидается при нахождении их под действием только гравитационных сил видимого распределения звезд и светящейся материи галактики. Проблема кривых вращения галактик — это несоответствие наблюдаемых скоростей вращения материи в дисках спиральных галактик и предсказаний динамики, учитывающих только видимую массу [27]. В настоящий момент считается, что это несоответствие выдает присутствие "скрытой массы", которая пронизывает галактику и простирается до галактического гало. Стандартный подход к моделированию кривой вращения галактики — суммирование ее основных компонент (одно- или многокомпонетного диска, балджа, гало). В общем случае результаты моделирования неоднозначны, и форма кривой вращения может быть удовлетворительно получена при различных соотношениях и распределениях протяженных дисковых и сфероидальных компонент. Изучение кинематики и пространственного распределения вещества в дисках других галактик показало, что поверхностная плотность и дисперсия скоростей экспоненциально уменьшаются с расстоянием от центра галактики, причем дисперсии остаточных скоростей пропорциональны поверхностной плотности [23].
Из общих соображений следует, что при построении модели по заданному потенциалу нужно стремиться включить в рассматриваемый объем модели область пространства, включающую 90-99% полной массы, соответствующей потенциалу. Тем самым удастся приблизить не только плотность в нужном районе пространства, но и общий гравитационный потенциал. Для построения фазовых моделей по заданному гравитационному потенциалу или распределению плотности был применен метод Шварцшильда [49]. В результате работы метода получаются фазовые модели в виде набора орбит, которые запоминаются в форме начальных данных. В этом дискретном наборе орбит каждая траектория, в общем случае не замкнутая, имеет определенную "массу" и рассматривается как бесконечная последовательность "звезд" (материальных точек), следующих друг за другом через равные промежутки времени в моделируемом потенциале. В пределе при уменьшении промежутков времени орбиту можно рассматривать как материальную нить с переменной плотностью на ее протяжении. Набор материальных "нитей" заполняет пространство модели, и подбирая их массы, требуется воспроизвести распределение вещества в модели, соответствующее моделируемому потенциалу. Для нахождения масс траекторий и контроля получаемого распределения плотности производится разбиение пространства модели на замкнутые ячейки и последующее вычисление масс внутри них, при этом сравниваются массы, вычисленные по заданному гравитационному потенциалу и получающиеся в процессе построения модели. Метод относится к так называемым "алгоритмам грубой силы", когда результат получается в процессе выбраковки элементов, не удовлетворяющих некоторой системе критериев. В идеале после выполнения всех операций остаются только необходимые элементы. В методе Шварцшильда такими элементами являются траектории в фазовом пространстве, а критерии — распределение пространственной плотности и, возможно, некоторые кинематические характеристики, если они включаются в набор условий. Такой подход к построению имеет как сильные, так и слабые стороны.
Преимуществом является то, что не нужно делать никаких предположений относительно существования дополнительных интегралов движения в заданном моделируемом потенциале и потенциал может быть практически любой. Минусом метода является необходимость использовать большое количество орбит и вычислять каждую траекторию до тех пор, пока не стабилизируется распределение интересующих нас характеристик в ячейках, пересекаемых рассматриваемой орбитой. Метод предполагает вычисление и обработку статистики для всех этих орбит, в том числе и для тех, которые впоследствии будут отброшены. Таким образом, какая-то часть вычислений оказывается ненужной в финальной модели, хотя отказаться от нее невозможно, потому что на этих этапах еще не известно, будет ли орбита использована в итоговом распределении. В случае недостаточного количества орбит в исходном распределении появляется риск получить недостаточно заполненные участки в фазовом пространстве, что отрицательно скажется на процессе построения модели или вообще не удастся удовлетворить критериям с требуемой точностью. Как правило, целью является построение равновесной и устойчивой модели с заданным потенциалом. При этом равновесной системой мы будем называть систему частиц—звезд, в которой макроскопические (средние по пространству) параметры с течением времени не меняются, в то время как сами частицы, составляющие систему, перемещаются, и их координаты и скорости постоянно меняются. Равновесная система может быть как устойчивой, так и неустойчивой. Устойчивой будем называть систему, которая при малом возмущении стремится вернуться в исходное состояние или имеет близкое состояние в малой окрестности от исходного. Если возмущение продолжает расти со временем, то такая система не устойчива.
В пользу получения равновесной модели, построенной методом Шварцшильда, имеются следующие предпосылки: набор орбит, являющихся базовой основой построения модели, вычисляется в заданном фиксированном потенциале, который должен быть воспроизведен в полученной модели. Распределение плотности вдоль траекторий считается обратно пропорциональным скорости материальной точки на ней. После отбраковки части орбит, суперпозиция оставшихся траекторий в соответствии с методом Шварцшильда должна аппроксимировать распределение массы в ячейках, которое вытекает из исходного заданного потенциала. В пределах погрешности потенциал суперпозиции "материальных" орбит будет близок к моделируемому потенциалу, при этом масса распределена вдоль орбит с учетом темпа движения. Поэтому "звезды", запущенные с начальными данными отобранных траекторий, будут двигаться в построенном потенциале и в моделируемом потенциале схожим образом, воссоздавая потенциал, который требуется при заданном распределении плотности. В итоге должна получиться самосогласованная система с потенциалом, близким к заданному гравитационному потенциалу. Однако поведение построенной модели может отличаться от описанного выше. Во-первых, потенциал суперпозиции материальных нитей—орбит не будет в точности соответствовать моделируемому потенциалу: расхождение будет тем больше, чем менее плотно расположены исходные орбиты в фазовом пространстве, и чем более грубая пространственная сетка ячеек используется для аппроксимации пространства модели и вычисления весов орбит. Во-вторых, в этом слегка искаженном гравитационном потенциале траектории будут несколько иначе заполнять пространство модели, что, в принципе, может привести к еще большему отклонению от исходного заданного потенциала и от самосогласованной модели. В общем случае процесс рассогласования может нарастать и параметры могут уйти от исходных значений очень далеко. Такая опасность есть у большинства методов построения фазовых моделей звездных систем, и метод Шварцшильда не гарантирует в этом смысле устойчивость. Поэтому на завершающем этапе построения необходимо дополнительное тестирование построенной модели.
Построение библиотеки орбит
Метод Шварцшильда подразумевает построение фазовой модели звездной системы на основе представленного набора орбит, из которых производится отбор таким образом, чтобы совокупность выбранных орбит воспроизводила моделируемый потенциал наилучшим образом, возможным для данного набо pa орбит. В общем случае, если нет дополнительных условий по кинематическим характеристикам, то приблизить заданный потенциал можно разными наборами орбит. Получившиеся фазовые модели будут отличаться кинематическими характеристиками, поскольку включают разные совокупности орбит, и в общем случае могут отличаться не только по кинематике, но и по устойчивости. "Материал" для нахождения того или иного набора орбит содержится в библиотеке орбит, вычисленных в заданном потенциале. Построение приемлемой фазовой модели будет возможно только в том случае, когда необходимые "компоненты" фазового пространства будут присутствовать в предоставляемом для метода Шварцшильда исходном наборе (библиотеке) орбит. Поэтому вычисление достаточно обширного набора траекторий (библиотеки орбит) с подходящими свойствами является важной и ответственной задачей. В стационарной осесимметричной системе в общем случае имеются два интеграла движения — интеграл энергии Е и интеграл площадей J, и естественно использовать эти интегралы движения как параметры для нахождения начальных координат орбит. Аналитический вид третьего интеграла движения найден только для нескольких частных случаев (например, третий квадратичный по скоростям интеграл Г.Г. Кузмина). Поэтому для перебора различных возможных вариантов орбит следует рассматривать различные начальные условия для траекторий при заданных значениях Е и J. Такой подход с использованием интегралов Е и J для создания библиотеки орбит с некоторыми видоизменениями был применен при построении всех моделей в данной работе. Основными действиями при формировании библиотеки орбит были следующие. На основе максимального заданного размера системы Rmod в экваториальной плоскости и по заданному потенциалу вычислялась верхняя граница диапазона интеграла площадей:
В диапазоне от 0 до Jmax выбиралось некоторое заданное количество значений интеграла площадей Jj и вычислялись соответствующие им значения Щ радиусов круговых орбит RCir, лежащих в экваториальной плоскости модели. Выли рассмотрены разбиения по интегралу площадей J и по радиусам Д„г как равномерные, так и с переменным шагом, в частности, с равномерным по логарифму, чтобы получить более плотное разбиение в центральной части модели и более редкое на периферии. Число значений Jj менялось в различных вариантах от 100 до 300; применялся как равномерный, так и переменный шаг. На следующем этапе вычислялся диапазон возможных значений интеграла энергии Е. Для каждого значения интеграла площадей Jj этот диапазон получается свой. Значение интеграла энергии ограничивалось снизу энергией, соответствующей движению по круговой орбите верхняя граница определялась условием нахождения траектории в пределах заданной границы модели Rm0d Конкретный набор значений интеграла энергии Ej выбирался из интервала {Emin,Emax). Разбиение интервала, как и в случае с интегралом площадей, использовалось как равномерное по Е, так и равномерное по логарифму Е. Результаты экспериментов показали, что целесообразнее использовать неравномерное дробление на 5-10 значений, при этом значение знаменателя прогрессии брать равным 2-3. Максимум фазовой плотности в полученных моделях при равномерном распределении по Е смещен к области со значениями энергии, соответствующими орбитам, близким к круговым; доля орбит с большими эксцентриситетами заметно меньше. Использование неравномерного распределения Е для построения библиотеки орбит существенно увеличивает разнообразие траекторий, к тому же, как показали численные эксперименты, использование такого несложного приема существенно упрощает работу модуля выбора траекторий и вычисления весов орбит. Задание двух интегралов движения Е{ и Jj не достаточно для определения траектории, так как оно определяет целую тороидальную область возможных движений. Дополнительными характеристиками, использующимися для получения и описания наборов орбит, в разных вариантах и в различных моделях были: а) расстояние от оси системы в экваториальной плоскости Rk, являвшееся начальным радиусом запуска траектории, который менялся от положения круговой орбиты до границы, определяемой поверхностью нулевой скорости для конкретной пары значений Jj, Ef, б) угол запуска орбиты по отношению к экваториальной плоскости XY (энергия вертикальных осцилляции). В варианте а) значения / являются аналогами дополнительного третьего интеграла, когда форма и ограниченность области, заполняемой витками траектории, указывает на его существование. Однозначная связь между траекторией и Rk определяется тем, что траектория при этом запускалась перпендикулярно экваториальной плоскости XY со значением интеграла энергии, равным Ei. Величина тангенциальной скорости определяется из интеграла площадей: Э = Jj/R . Орбита при таком способе запуска однозначно описывается тремя значениями /,, Ei, Rk, поэтому после вычисления весов при построении дискретной модели положения и скорости точечных масс соответствовали первоначальным орбитам. Однако у такого "вертикального" способа запуска имеется существенный недостаток: получается избыток траекторий с большими компонентами скорости по z и недостаток траекторий, прилегающих к экваториальной плоскости.
Это ведет к завышенному значению дисперсии вертикальной компоненты скорости в итоговом распределении, что в свою очередь ведет к завышенной полутолщине диска в построенной модели. Поэтому в варианте б) в качестве третьего параметра было решено использовать угол запуска а к экваториальной плоскости XY (5-7 диапазонов значений с неравномерным шагом). Сама начальная точка траектории выбиралась в плоскости XY (z = 0). Радиус запуска R в этом варианте не фиксировался, а выбирался равномерно случайно из интервала по R, определяемого внутренней и внешней границами контура нулевой скорости Е. + U{R о) _ А- = о. (3.4) Для того, чтобы плотнее заполнить область фазового пространства, соответствующую определенным значениям [Jj, Ei, а ), для каждой тройки запускалось по несколько орбит и статистика по ним суммировалась. В выполненных нами построениях моделей использовалось по 15-20 орбит. Недостатком при такой генерации орбит и наборе статистики является отсутствие запоминания параметров конкретных траекторий. Запоминаются только параметры случайной генерации траекторий и результаты статистики пересечений ячеек этими траекториями. В результате на этапе замены построенной модели на дискретную модель из точечных масс не получается абсолютно точное соответствие фазовому распределению использованных орбит из-за наличия в построении случайной составляющей. Однако использование большого числа орбит (15-20) и точечных масс (от 5 до 100 в зависимости от веса
Численное решение задачи N тел
Для исследования поведения системы N частиц под действием создаваемого ими гравитационного притяжения и, если требовалось в некоторых моделях, воздействия дополнительного внешнего гравитационного поля, использовался иерархический алгоритм (tree code) [56, 16]. Применение этого алго ритма для прослеживания эволюции системы N тел позволяет существенно сократить временные затраты на вычисления при больших значениях N. Соотношение числа операций в иерархическом алгоритме составляет N\og(N) против ./V2 при прямом вычислении всех сил между телами. Выбор параметров иерархического алгоритма был сделан на основе рекомендаций из работы Родионова и Сотниковой [15]. Исследование моделей Галактики проводилось на интервале времени примерно 20-40 оборотов звезды, находящейся на галактоцентрическом расстоянии Солнца (5-10 млрд. лет). В начальный момент интегрирования и затем через определенные промежутки времени вычислялись общие характеристики системы, профили плотности, скоростей и дисперсий скоростей. Параметры определялись с помощью программного модуля, применявшегося при построении библиотеки орбит, причем использовалась та же сетка ячеек для разбиения пространства модели. Таким образом, можно сравнивать одни и те же характеристики системы, вычисленные на различные моменты времени, и проследить эволюцию параметров модели. После проведения численных экспериментов с полным циклом построения фазовой модели для потенциала Кутузова—Осипкова при хорошей аппроксимации заданных потенциала и распределения плотности были выявлены некоторые кинематические особенности получившихся первых моделей: они имели большие дисперсии1 скоростей, величины которых были порядка 80-120 км/с на расстоянии 10 кпк от центра модели, а кривая вращения на плоском участке соответствовала линейной скорости порядка 150-180 км/с. Обе эти характеристики не согласуются с данными наблюдений для галактической окрестности Солнца. Причина появления таких кинематических характеристик была связана с начальной библиотекой орбит, из которой выбирался набор орбит для создания фазовой модели. В ней оказалось завышенное ко личество траекторий с большими значениями компоненты скорости по z, то есть, условно говоря, орбит с преувеличенной энергией Ez вертикальных колебаний.
Появление избыточного числа орбит с большими скоростями по Z и недостаточного числа прилегающих к экваториальной плоскости траекторий было связано с использованным алгоритмом формирования библиотеки орбит. Максимальное количество орбит и ячеек ограничивается сверху имеющимися вычислительными ресурсами, и в наших построениях использовались значения для них 104 с учетом приемлемых времен нахождения набора траекторий, аппроксимирующего заданные характеристики модели. С помощью такого количества орбит необходимо было заполнить по возможности большую часть фазового пространства модели. Для создания библиотеки орбит использовались распределения интеграла площадей J (радиусов круговых орбит Ясгг в экваториальной плоскости), интеграла энергии Е и набора значений радиальных расстояний /. В сопутствующей плоскости все траектории стартовали ортогонально оси R, то есть имели нулевую начальную компоненту скорости по Я и максимальную скорость по z при данных значениях J, Е и /4, что привело к большому числу траекторий с завышенными значениями энергии Ez "вертикальных" осцилляции. Это особенно важно для таких быстро вращающихся подсистем, как дисковая подсистема Галактики, где существенно распределение орбит по интегралу площадей. Оно определяет положение (расстояние от оси вращения) круговых траекторий в экваториальной плоскости. Некруговые траектории с одинаковым значением интеграла площадей находятся внутри тора, который окружает круговую орбиту с данным значением J и размеры которого зависят от интеграла энергии Е. Чем больше энергия Е при данном фиксированном интеграле площадей, тем сильнее траектория отличается от круговой.
Энергия Ez вертикальных колебаний определяет, будет ли орбита расположена вблизи плоскости диска или отходит от нее в область больших z. Использование в первоначальных вариантах одинакового числа делений по каждому из этих параметров J, Е и Es не дало удовлетворительного результата и привело к завышенной дисперсии скоростей в диске. Количество круговых орбит получалось порядка 20 при радиусе моделируемой области 40 кпк и круговые орбиты в окрестности Солнца в среднем отстояли друг от друга на расстоянии 2 кпк. Ячейки же были расположены гораздо более плотно, так как их размеры в радиальном направлении были 0.2-0.3 кпк. Чтобы заполнить ячейки, в которые не попадали круговые и близкие к круговым орбиты, приходилось привлекать орбиты с большой энергией (более "размазанные" в пространстве). Естественно, при этом дисперсия скоростей в полученной модели возрастала. Чтобы устранить описанный дефект модели, круговые орбиты, вокруг которых формируются области возможных движений, должны быть расположены чаще. Расстояние между ними можно оценить следующим образом: дисперсия радиальных компонент скоростей звезд в Галактике в окрестности Солнца оценивается величиной порядка 30-40 км/с, это соответствует колебаниям радиуса орбиты на величину порядка 1 кик, следовательно, и круговые орбиты должны быть расположены не реже. В последующих вариантах алгоритмы формирования библиотеки орбит были изменены, чтобы обеспечить достаточно широкий спектр траекторий.
