CC BY
20
информатика
Батюков Александр Михайлович
УДК 004.942
ОБ АЛГОРИТМАХ МОДЕЛИРОВАНИЯ ПРОЦЕССОВ АГРЕГАЦИИ, ОГРАНИЧЕННОЙ ДИФФУЗИЕЙ*
Аннотация
Описаны реализации классического алгоритма моделирования процессов агрегации, ограниченной диффузией, как на плоскости, так и на триангуляционной сетке. Разработаны и реализованы оптимизации этих алгоритмов, основанные на априорной оценке вероятности присоединения частицы к граничным точкам агрегата. Дана оценка сложности описываемых алгоритмов. Приведены результаты численных экспериментов.
Ключевые слова: процессы агрегации, ограниченной диффузией; алгоритмы на триангуляционной сетке; априорная оценка вероятности присоединения.
ВВЕДЕНИЕ
Многие процессы, происходящие в природе и обществе, обладают хаотической динамикой, которая характеризуется возникновением фрактальных структур. Подобные структуры, характеризуемые статистической (неполной) повторяемостью, часто называют хаотическими фракталами. Важный класс среди них составляют так называемые фрактальные кластеры — класс физических объектов, плотность которых уменьшается по мере их роста. Исследование фрактальных кластеров началось с появления теоретической модели агрегации, ограниченной диффузией (diffusion limited aggregation, DLA), которая описывает агрегацию частиц в условиях их случайного блуждания. Позже были предложены и другие модели, такие как cluster-cluster aggregation (CCA) и reaction-limited cluster aggregation (RLCA). Их подробный обзор можно найти в [3].
Модель diffusion-limited aggregation (DLA) описывает рост групп частиц (назы-
© Батюков А.М., 2014
ваемых агрегатами) при их случайных (броуновских) блужданиях. Эта модель была предложена Виттеном, Сандлером и Пинским [5, 7] и оказалась широко применимой для имитационного моделирования таких явлений, как осаждение металла при электролизе, распространение электрического разряда при пробое диэлектрика, протекание диффузионных процессов в жидкостях и газах.
При реализации алгоритмов моделирования DLA расчёты занимают значительное количество времени. Поэтому существуют различные варианты оптимизации алгоритма Э. Пинского, которые позволяют снизить время вычислений, требуемых для построения агрегата из большого количества частиц. Описания некоторых из них можно найти в [6]: одни основаны на том, что частицам позволяется за один шаг перемещаться на расстояния, большие 1, другие - на том, что при бросании частицы на плоскость случайным образом ей задается начальный вектор
* Работа выполнена при частичной поддержке гранта РФФИ 13-01-00782.
поступательного движения, который на каждом шаге блужданий добавляется к свободным перемещениям частицы.
В данной работе описана оптимизация классического алгоритма DLA на плоскости, основанная на использовании неполных данных при вычислениях с помощью вычисления априорной оценки вероятности присоединения бросаемой на плоскость частицы к каждой из граничных точек агрегата [1]. Как показывает численный эксперимент, она позволяет сократить время построения агрегата в среднем более чем в 10 раз.
Естественным обобщением алгоритма построения модели DLA на плоскости является его перенос на поверхность произвольной формы. Такая адаптация была выполнена в [2], где поверхность была приближена с помощью нерегулярной триангуляционной сетки, на которой, согласно модели, и протекают процессы диффузии. Воспроизведение алгоритма, описанного в этой статье, показывает, что для построения агрегата из достаточно большого количества частиц, как и в плоском случае, требуется значительное время вычислений.
В данной работе мы предлагаем адаптацию метода оптимизации процесса моделирования DLA на нерегулярной триангуляционной сетке с использованием априорной оценки вероятности. Метод учитывает описанные в статье [2] неравновероятные переходы между треугольниками триангуляции. Результаты эксперимента показывают, что при многократных построениях агрегата на той же поверхности оптимизированный ме-
Л ^
Рис. 1. Пример построенного агрегата
тод позволяет значительно сократить время вычислений.
1. ОПИСАНИЕ МОДЕЛИ DLA НА ПЛОСКОСТИ
Согласно классической модели DLA, блуждание частиц происходит по линиям дискретной прямоугольной сетки. Кроме того, предполагается, что частицы присоединяются к агрегату последовательно, одна за одной. Процесс построения агрегата начинается с единственной частицы, располагающейся в начале координат прямоугольной сетки. Случайным образом вычисляется первоначальное положение другой частицы. Эта частица начинает блуждание по линиям сетки, при этом на каждом шаге она может с вероятностью 1/4 переместиться либо на единицу вправо, либо на единицу влево, либо на единицу вниз, либо на единицу вверх относительно своего текущего положения. Блуждание продолжается до тех пор, пока частица не окажется соседней с какой-либо из частиц, входящих в агрегат. Остановившаяся частица соединяется отрезком с агрегатом и, начиная с этого момента, входит в него. После этого происходит бросание новой частицы. Описанный процесс повторяется многократно. Пример построенного агрегата показан на рис. 1.
При практической реализации вводят дополнительные ограничения: бросают частицу в заданную окрестность агрегата и ограничивают число шагов ее блуждания некоторым числом N. Если за N шагов частица не присоединилась к агрегату, ее отбрасывают.
2. ОПТИМИЗАЦИЯ МОДЕЛИ DLA НА ПЛОСКОСТИ
В работе [1] описана оптимизация, основанная на априорном вычислении оценки вероятности присоединения бросаемой на плоскость частицы к каждой из граничных точек агрегата.
Предположим, что при бросании частицы на координатную сетку ее координаты оказались равны (х + а, у + Ь). Сделаем оценку вероятности того, что данная частица присоединится к точке агрегата с координата-
ми (х,у). Описываемая ситуация показана на рис. 2. Если при блуждании частицы по плоскости она сделала п шагов в сторону от точки (х, у) по оси абсцисс, то, чтобы попасть в (х, у), частица должна сделать (п + a) шагов в сторону точки (х, у). Соответственно, если при блуждании частица сделала т шагов в сторону от точки (х, у) по оси ординат, то она должна сделать (т + Ь) шагов в сторону точки (х, у). На каждом шаге блуждания сторона следующего перехода определяется равновероятно между четырьмя вариантами. Это позволяет дать следующую оценку вероятности того, что через (2п + a + 2m + Ь) шагов координатами блуждающей точки окажутся (х, у):
1
^и 44+а 4т 4т+Ь •
Число шагов, через которые блуждающая точка попадет в (х, у), определяется значениями и и т. Классы траекторий, определяемых конкретными значениями и и т, представляют собой разбиение множества всех траекторий на непересекающиеся классы. Поэтому, чтобы получить оценку для всего множества траекторий, подсчитаем сумму ряда:
1
II
4П 4П+ а 4т4т+Ь
= 11-
т= 0 и= 0
I 2и+2т+а+Ь
1 а+Ь II
т=0 и=0
4а+ь ^^ ^^ 1 /гт+ и
т = 0 и = 0 16
= 4а+Ь 1 16т 1 16т = [ 15 ) 4а+Ь .
Полученное значение дает верхнюю оценку величины, определяющей среднюю скорость, с которой блуждающая точка с координатами (х + а, у + Ь) может попасть в точку (х, у).
При бросании частицы на плоскость вычисляют такие оценки между первоначальными ее координатами и всеми точками границы агрегата. Положим, 1
Р(а, Ь) = -
1 а+Ь
и полученный набор чисел
будем называть коэффициентами выбора
Рис. 2. Шаг процесса блужданий
присоединения. В соответствии с вычисленными коэффициентами строится функция распределения вероятности присоединения частицы к граничным точкам агрегата -«функция численного распределения» [4], из области значений которой случайным образом выбирается значение. Прообраз этого значения и определяет, к какой точке агрегата присоединится частица.
Нужно заметить, что, пользуясь априорными оценками, мы теряем часть информации, связанную с тем, каким именно путем частица из своего начального положения присоединилась к агрегату. От получения этой информации мы отказываемся ради ускорения вычислений.
При реализации этого алгоритма для ограниченной области из М точек на сетке предварительно вычисляются значения коэффициентов выбора присоединения для каждой пары точек сетки. Из них можно составить матрицу G размера М х М вероятностей переходов из одной точки в другую. Заметим, что для двух пар точек, таких что в каждой паре суммарное покоординатное расстояние между точками равно 5 = а + Ь, значение коэффициента выбора будет одинаковым и равным 1/45. Это означает, что построение матрицы G избыточно, и достаточно вычислить значение коэффициента выбора для всех возможных 5, что позволяет в каждый момент времени использовать не более, чем О (М) памяти и построить его за О (М) операций.
В условиях задачи при бросании Т частиц сложность построения агрегата по вычисленным априорным оценкам составляет О (ТМ).
т=0 и=0
1
1
3. АДАПТАЦИЯ МОДЕЛИ DLA
ДЛЯ НЕРЕГУЛЯРНОЙ ТРИАНГУЛЯЦИОННОЙ СЕТКИ
В работе [2] описана реализация алгоритма DLA на нерегулярной триангуляционной сетке. Точно так же, как и в плоском случае, процесс формирования агрегата начинается с одной частицы, занимающей один из треугольников триангуляции. Другие частицы последовательно бросаются на треугольники, составляющие триангуляцию, и начинают блуждать по ней. На каждом шаге частица может перейти в один из треугольников, соседних с треугольником, задающим ее текущее положение. Переходы в разных соседей не равновероятны - вероятность перехода из треугольника 1 в треугольник 2 определяется по формуле
Р(1 ^ 2) = --^--
1 а + 1/Ь + Vс ,
где величины а, Ь, c - длины сторон треугольника 1, отмеченные на рис. 3.
В качестве примера авторы [2] успешно моделируют процесс распостранения вещества на модели поверхности кости. Тем не менее, для расчетов агрегаций, составленных из большого количества частиц, как и в плоском случае, требуется значительное время работы вычислителя.
4. ОПТИМИЗАЦИЯ МОДЕЛИ DLA НА НЕРЕГУЛЯРНОЙ ТРИАНГУЛЯЦИОННОЙ СЕТКЕ
Построим на выделенной области поверхности триангуляционную сетку. Как и в [2], треугольники в ней могут иметь разные раз-
С
Рис. 3. Определение вероятности перехода между элементами триангуляционной сетки
меры. Пусть общее число треугольников равно М.
Каждому треугольнику ставим в соответствие вершину некоторого графа G. В этом графе между вершинами i и у есть дуга тогда и только тогда, когда соответствующие им треугольники М1 и Му имеют общую сторону. На дугах расставим веса, равные вероятности перехода из i в у согласно формуле (1).
Рассмотрим всевозможные пути на графе длины не более, чем некоторое N. Для любого пути между двумя произвольными вершинами i и у вероятностью этого пути назовем произведение вероятностей (весов) всех входящих в путь дуг. Вычислим некоторый коэффициент выбора для вершин i и у, равный сумме вероятностей всевозможных путей из i в у длины не более, чем N. Из всех коэффициентов выбора составим матрицу G размера М х М.
Вычисления организуем следующим образом. Изначально нам даны вероятности всех путей, состоящих ровно из одной дуги (это веса на соответствующих дугах). Будем последовательно вычислять некоторые матрицы Gk. В матрицу G1 на позицию (/',у) при наличии на графе дуги (/', у) запишем значение веса этой дуги. На каждом шаге будем из матрицы Gk составлять матрицу Gk +1 следующим образом. Пусть элемент (/, у) матрицы Gk равен у. Рассмотрим все дуги на графе, начало которых - у (таких дуг не более 3 -по числу возможных соседей). Пусть для некоторой из этих дуг ее конец - х. Тогда в матрицу Gk +1 на позицию (/', х) добавим число, равное у, умноженному на вероятности дуги (у, х). Матрицей G назовем сумму всех матриц Gk. Несложно заметить, что матрица G содержит значения всех коэффициентов выбора для путей длины не более, чем к.
Сложность указанного алгоритма в условиях нашей задачи О ^М2). В каждый момент времени в памяти необходимо держать матрицы G1, Gk и матрицу G (вследствие того, что наш процесс марковский, не нужно держать в памяти промежуточные G), поэтому размер используемой памяти пропорционален М 2.
Моделирование бросания частицы на триангуляционную сетку будем производить, как и в плоском случае, - в соответствии с коэффициентами выбора строить функцию неравномерного распределения вероятности присоединения частицы к границе агрегата, случайным образом выбирать точку из ее области значений и по прообразу определять, к какой точке агрегата присоединится частица. При бросании Т частиц сложность вычислений составит, как и ранее, О(ТМ).
Можно сказать, что построенная матрица G, как в плоском случае, так и в случае моделирования на поверхности, содержит в себе информацию о геометрии исходной поверхности.
Геометрия плоской поверхности проста, поэтому матрица G для нее имеет очень простой вид. В случае моделирования на нерегулярной триангуляционной сетке матрица G не избыточна и не может быть без потерь информации сохранена в структуре данных меньшей размерности.
Безусловно, подсчет матрицы G для больших М трудоемок, и в случае однократного вычисления агрегата на конкретной поверхности время вычислителя, может быть сравнимо с временем вычислений по классическому алгоритму. Но обычно задачу построения агрегатов решают для того, чтобы оценить динамику протекающих на
поверхности процессов. Это означает, что для конкретной поверхности для сбора статистики необходимо построить достаточное число (сотни) агрегатов. И в этом случае важно то, что для каждой поверхности матрицу G можно подсчитать 1 раз, после чего ее можно использовать без изменений для каждого нового вычисления. Это позволяет существенно экономить время вычислений при многократной операции построения агрегата на поверхности.
5. ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ
Реализованы два алгоритма построения агрегата на плоскости и два алгоритма построения агрегата на нерегулярной триангуляционной сетке. В каждой паре один из них соответствует классическому построению модели DLA, а второй - алгоритму с вычислением априорной вероятности.
Произведены две серии экспериментов. В первой были построены с помощью двух алгоритмов по 10 агрегатов на плоскости. Во второй с помощью двух алгоритмов были построены по 10 агрегатов на двух триангуляционных сетках для поверхностей
2 2 2 х , у 2 , х2 2 -
-+ ---г =1 и-+ у = 2 г .
4 4 4
Для вычислений были взяты параметры N = 300, М = 5000, Т = 1000. Время, затраченное на вычисления, показано в табл. 1.
Табл. 1. Время построения агрегатов в сериях вычислений
Классический алгоритм Оптимизиро ванный алгоритм
Время вычисления 1 агрегата на плоскости 37 мин. 43 сек 1 мин. 8 сек
Общее время вычислений 10 агрегатов на плоскости 6 ч. 15 мин. 22 сек. 1 мин. 54 сек.
Время вычисления 1 агрегата на поверхности 1 26 мин. 12 сек. 29 мин. 43 сек.
Общее время вычисления 10 агрегатов на поверхности 1 4 ч. 47 мин. 33 сек. 31 мин. 58 сек.
Время вычисления 1 агрегата на поверхности 2 24 мин. 59 сек. 28 мин. 11 сек.
Время вычисления 10 агрегатов на поверхности 2 4 ч. 11 мин. 37 сек. 32 мин. 4 сек.
6. ЗАКЛЮЧЕНИЕ
Разработаны и реализованы оптимизации алгоритма построения агрегата по модели DLA на плоскости и на нерегулярной триангуляционной сетке. Эти оптимизации особенно эффективны при необходимости
многократных построений агрегатов на одной исходной поверхности. Можно ожидать, что подобная оптимизация может быть использована для других моделей построения агрегата на плоскости и на нерегулярной триангуляционной сетке.
Литература
1. Батюков А.М., Ампилова Н.Б. О модификации модели агрегации, ограниченной диффузией // Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление, 2013. Вып. 3 (174). С. 233-238.
2. Евсеев А.А., Нечаева О.И. Клеточно-автоматное моделирование диффузионных процессов на триангуляционных сетках // Прикладная дискретная математика, 2009. Вып. 4(6). С. 72-83.
3. Смирнов Б.М. Физика фрактальных кластеров. М.: Наука, 1991.
4. Knuth D. The Art of Computer Programming. Vol.2. Seminumerical Algorithms, 3-ed edition. 2001.
5. Pinski A. Diffusion-limited aggregation // Physical Review. The American Physical Society. USA, 1981.
6. Vicsek T. Fractal Growth Phenomena // World Scientific Pub Co Inc; 2 Sub edition, 1992.
7. Witten T.A., Sander L.M. Diffusion-limited aggregation // Physical Review. The American Physical Society. USA, 1982.
ON MODELINGN OF DIFFUSION LIMITED AGGREGATION PROCESSES
Abstract
The implementations of diffusion-limited aggregation algorithms both on plane and a surface approximated by irregular triangulation net are described. The optimization of the method based on the a priori estimation of the sticking probability of a particle to the aggregation (joining coefficient) has been proposed and realized. The algorithm complexity estimation and results of numerical experiments are described.
Keywords: diffusion-limited aggregation processes, algorithms on triangulation net, a priori estimation of joining probability.
(q) Наши авторы, 2014. Our authors, 2014.
Батюков Александр Михайлович, аспирант кафедры информатики математико-механи ческого факультета СПбГУ, batsun@gmail. com.
