Диссертация (1145359), страница 6
Текст из файла (страница 6)
В случае сглаженного потенциала таких проблем не возникает.Во-вторых — и это более важно — процедура сглаживания уменьшает “зернистость” распределения частиц и тем самым приближает потенциал модельнойсистемы к потенциалу системы с более гладким распределением плотности, т.е.системы, описываемой бесстолкновительным уравнением Больцмана.Очевидно, что величина не может быть слишком большой, поскольку, содной стороны, это приводит к значительному искажению потенциала, а с другой — ставит сильные ограничения на пространственное разрешение тех илииных структурных особенностей системы.
При этом важно иметь объективныйкритерий выбора параметра сглаживания в N -body экспериментах. В данномразделе такой критерий обосновывается, исходя из анализа временны́х изменений функций распределения для сферически-симметричных моделей при различных .Критерий выборапо МеритуМерит [56] предложил критерий выбора параметра сглаживания, основанный на минимизации средней иррегулярной силы, действующей на частицу.
Онввел величину M ISE (mean integrated square error), которая характеризует отличие силы, создаваемой дискретным набором N -тел, от силы, действующей всистеме с непрерывным распределением плотностиZM ISE() = Eρ(r)|F(r, ) − Ftrue (r)| dr ,2(1.25)где E — усреднение по многим реализациям системы, F(r, ) сила, действующаяна частицу единичной массы, находящуюся в точке r, со стороны N частиц сосглаженным потенциалом, Ftrue (r) — сила, действующая на эту же частицу в36системе с непрерывным распределением плотности.Рис. 1.7. Зависимость величинысM ISEот параметра сглаживания,для модели ПламмераN = 1000.Имеются два фактора, которые вносят вклад в M ISE :- флуктуации дискретного распределения; они очень велики при малых значениях , уменьшаются с ростом и для данного тем меньше, чем большеN (рис. 1.7, левая ветвь кривой);- ошибки в представлении потенциала (отличие сглаженного потенциала отпотенциала точки); эти ошибки, наоборот, очень велики при больших ,падают с уменьшением и не зависят от N (рис.
1.7, правая ветвь кривой).При некотором значении = m величина M ISE будет минимальной(см. рис. 1.7). Данное m Мерит предложил считать оптимальным выборомпараметра сглаживания в N -body экспериментах. Эта величина зависит отN и, например, для модели Пламмера аппроксимируется зависимостью m ≈0.58N −0.26 (если использовать вириальную систему единиц, в которой G = 1,полная масса модели Mpl = 1, полная энергия системы Epl = −1/4) в диапазоне N от 30 до 300 000 [57].37Постановка задачиКритерий Мерита основан на минимизации средней ошибки в представлении силы в равновесной системе в начальный момент времени.
Мы поставилизадачу вывести критерий выбора параметра сглаживания непосредственно издинамики системы. Если описывать равновесную устойчивую конфигурациючерез дискретный набор N тел, то из-за разнообразных ошибок параметры системы (например, эффективный радиус) будут со временем отклоняться от начальных значений. Чем меньше величина этих отклонений, тем более адекватномоделируется динамика системы. Можно при этом ожидать, что существует такое значение , которое минимизирует эти изменения.Нами исследовались две равновесные устойчивые модели: сфера Пламмеракак случай модели с распределением плотности, близким к однородному3Mplρ(r) =4πa2plr2 +a2pl5/2(1.26)и сфера Хернквиста [47] как пример модели с сильно неоднородным распределением плотностиahrMhr,(1.27)2π r(r + ahr )3где Mpl и Mhr — полные массы моделей Пламмера и Хернквиста, соответственρ(r) =но, а apl и ahr — характерные масштабы распределения плотности для моделиПламмера и Хернквиста, соответственно (для модели Пламмера внутри радиусаapl сосредоточено примерно 35% массы, для модели Хернквиста внутри радисуса ahr сосредоточено 25% массы).
Для обеих моделей предполагалось изотропное распределение по скоростям. Использовалась вириальная система единиц(G = 1, Mpl = Mhr = 1, Etot = −1/4). В этих единицах apl = 3π/16 ≈ 0.59,ahr = 1/3.Изпользовался пакет NEMO (http://bima.astro.umd/nemo/, [58]), основанный на методе структурирования данных в виде “иерархического дерева”(TREE) при вычислении силы взаимодействия между частицами [59]. При этом38параметр θ, отвечающий в алгоритме TREE за точность вычисления силы, дляприведенных ниже расчетов равнялся 0.75.
Тестовые эксперименты, проделанные с большей (θ = 0.5) и меньшей (θ = 1.0) точностью вычисления силыпоказали, что для данного интервала параметра θ наши выводы от θ не зависят. Как продемонстрировано в работе [60], ошибки в представлении силы вметоде TREE меньше ошибок, связанных со сглаживанием потенциала, поэтому данный метод может быть использован для решения поставленной задачи.Для интегрирования уравнений движения частиц использовалась стандартнаясхема с перешагиванием (leap frog), которая обеспечивает второй порядок точности по времени. Число частиц N = 10 000. Для задания начальных равновесныхконфигураций также использовался пакет NEMO.Мы применяли два вида сглаживания потенциала.
Первый метод основанна замене потенциала точки на потенциал сферы Пламмера, что приводит кмодификации силы взаимодействия между частицами согласно схеме (1.24). Вовтором — для сглаживания потенциала использовался кубический сплайн [32].Мы приводим результаты только для первого метода сглаживания. Однако всенаши выводы сохраняются и для сглаживания сплайном.Отслеживалось, как сохраняются различные глобальные параметры системы, в частности мы анализировали изменение функции распределения частицв пространстве и отдельно изменение функции распределения частиц по скоростям, т.е.
фактически следили за изменениями фазовой плотности.Отклонение распределения частиц в пространстве в момент времени t отначального мы характеризовали величиной ∆r , которая вычислялась следующим образом: модель разбивалась на сферические слои, в каждом сферическомслое вычислялась разность между количеством частиц в момент времени t и вначальный момент. Далее величина ∆r вычислялась как сумма модулей этихразностей и нормировалось на суммарное количество частиц в системе N . Толщина слоя составляла 0.1, максимальный радиус равнялся 2 (при этом дляобеих моделей в каждом из сферических слоев было статистически значимое39количество частиц — больше 50). Заметим (и это важно): если взять две случайные реализации какой-либо модели, то для них ∆r не будет равно нулю.
Уэтой величины существует заметный естественный шум, связанный с конечнымчислом частиц. Уровень шума оценивался как среднее значение величины ∆rдля большого количества пар случайных реализаций модели.Параметр ∆v вычислялся аналогичным образом, но для распределениячастиц в пространстве скоростей. Для приведенных результатов толщина слояравнялась 0.1, а максимальная скорость выбиралась равной 1.5 (такой выбормаксимальной скорости для обеих моделей обеспечивал в каждом из сферических слоев статистически значимое количество частиц — больше 30).Вычислялось также время парной релаксации для исследуемых моделей взависимости от .
Оно определялось как время существенного отклонения частицы от первоначальной радиальной орбиты (критерий отклонения и способоценки во многом схожи с теми, что применялись в работе [61]). Кратко методоценки можно описать следующим образом. Строилась случайная реализациясистемы (сферы Пламмера или сферы Хернквиста). В этой модели на расстоянии примерно в один эффективный радиус от центра выбиралась частица спочти радиальной орбитой со скоростью, направленной к центру.
Орбита этойчастицы подправлялась до строго радиальной. Прослеживалась эволюция всейсистемы до того момента, пока помеченная частица не пройдет расстояние в 1.5раза больше своего первоначального расстояния от центра. Используя время tp ,прошедшее от начала движения, и угол отклонения частицы от своей первоначально радиальной орбиты (обозначим этот угол через αp ), можно вычислитьвремя, которое понадобится частице, чтобы отклониться от первоначального направления движения на один радиан.
Усреднив полученное время по большомуколичеству подобных частиц, мы получаем оценку времени парной релаксации1trelax=αp√tp2,(1.28)здесь trelax — время парной релаксации, < . . . > — усреднение по большому40количеству пробных частиц. Формула (1.28) справедлива в диффузионном приближении, т.е. в предположении о малых углах единичного рассеяния.Результаты численного моделированияСфера ПламмераРисунки 1.8–1.10 иллюстрируют результаты расчетов для сферы Пламмера. Использовалось N = 10 000. Для этого числа частиц оптимальное значениепараметра сглаживания по Мериту m = 0.05 в вириальных единицах. Минимальное значение параметра сглаживания, до которого был подстроен шагинтегрирования dt = 0.006. Усреднение для каждого значения проводилосьпо нескольким моделям.
Шаг интегрирования dt = 0.01. Перечислим основныевыводы.Рис. 1.8. Зависимость∆rи∆vот для разных моментов времени для сферы Пламмера. Гори-зонтальная линия показывает уровень естественного шума величинлиниями отмечено положениеmиdt .∆rи∆v . Вертикальными41Рис. 1.9. Среднее расстояние до ближайшего соседа на различных радиусах для моделейПламмера и Хернквиста. Вертикальные линии соответствуют радиусам сфер, внутри которых находится 0.01, 0.1 и 0.5 массы всей системы. Горизонтальная линия показывает оптимальное значениеm .Рис.
1.10. Время релаксации для сфер Пламмера и Хернквиста как функция . Вертикальнаялиния отмечает оптимальное значениеm .Шаг интегрированияdt = 0.001.1. При > m отличия сглаженного потенциала от ньютоновского становятся значительными. Система подстраивается под новый потенциал, изменяя свою функцию распределения. Подстройка происходит практическисразу — на временах t 20 (единица времени примерно соответствует вре-42мени пересечения ядра apl данной модели). В дальнейшем в такой сильно“сглаженной” и далекой от начального состояния системе практически непроисходит никаких изменений (рис.