CC BY
139
2015 Дискретные модели реальных процессов № 1(27)
ДИСКРЕТНЫЕ МОДЕЛИ РЕАЛЬНЫХ ПРОЦЕССОВ
УДК 621.391.1:004.7
РЕЖИМЫ ФУНКЦИОНИРОВАНИЯ АСИНХРОННЫХ
КЛЕТОЧНЫХ АВТОМАТОВ, МОДЕЛИРУЮЩИХ НЕЛИНЕЙНУЮ ПРОСТРАНСТВЕННУЮ ДИНАМИКУ1
О. Л. Бандман
Институт вычислительной математики и математической геофизики СО РАН,
г. Новосибирск, Россия
Сдвиг научного интереса от физических явлений, подчиняющихся законам термодинамики, к нелинейным диссипативным процессам, содержащим химические и биологические превращения, привёл к аналогичному повороту в математическом моделировании: от решения дифференциальных уравнений к прямому дискретному стохастическому моделированию. Математическим фундаментом дискретного моделирования является асинхронный клеточный автомат — стохастический аналог клеточного автомата фон Неймана. Систематической методологии построения асинхронного клеточного автомата, моделирующего процессы, состоящие из многих действий, совместно преобразующих общее дискретное пространство, пока не существует. Нет ответа на вопрос, насколько и чем различаются результаты моделирования при разных способах организации (режимах) взаимодействий локальных операторов, составляющих сложный процесс. В работе делается попытка ответить на этот вопрос путём проведения серии вычислительных экспериментов по моделированию трёх типовых реакционно-диффузионных процессов при разных асинхронных режимах и сравнительного анализа их эволюций. Результат состоит в том, что качественный характер процессов не зависит от способа композиции, а количественные различия могут быть скорректированы.
Ключевые слова: дискретное математическое моделирование, асинхронный клеточный автомат, режимы функционирования, реакционно-диффузионные процессы, пространственная самоорганизация.
FUNCTIONING MODES OF ASYNCHRONOUS CELLULAR AUTOMATA SIMULATING NONLINEAR SPATIAL DYNAMICS
O. L. Bandman
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk,
Russia
E-mail: bandman@ssd.sscc.ru
1 Работа поддержана грантами Президиума РАН (Проект 15.9, 2014) и СО РАН (Интеграционный проект ИП-47, 2011).
The shift of scientific interest from physical phenomena obeying laws of thermodynamics towards nonlinear dissipative processes containing chemical and biological transformations stimulates a similar turn in mathematical modeling: from differential equation solution to direct and stochastic simulation. A foundation for discrete simulation is the asynchronous cellular automaton — a stochastic analogue of von-Neumann's cellular automaton. For the time being, there is no systematic methodology for constructing asynchronous cellular automata simulating processes composed of many actions transforming a common discrete space. It is not known, how different are simulation results obtained by different ways of composing simple operations for organizing a complex computational process. In the paper, an attempt is made to answer this question by means of performing a series of simulation of three typical reaction-diffusion processes with different asynchronous modes of functioning, and comparative analysis of their evolutions and invariants. The obtained result shows that qualitative character of the process under simulation does not depend on the composition mode, and quantitative differences may be corrected.
Keywords: discrete mathematical modeling, asynchronous cellular automaton, modes of functioning, reaction-diffusion processes, spacial self organization.
Введение
Наряду с теоретическими и экспериментальными видами научной деятельности развивается третья её составляющая — компьютерное моделирование исследуемых процессов. В частности, с помощью имеющихся мощных компьютеров стало возможным изучение механизмов химических и биологических процессов на микро- и на-ноуровнях, которые трудно или невозможно моделировать, основываясь на традиционных моделях математической физики, так как исследуемые явления нелинейны, разрывны и диссипативны. Поиски новых нетрадиционных моделей привели к переосмысливанию идеи клеточного автомата (КА) [1] и появлению множества его модификаций. Среди них значительное место занимают КА-модели явлений, которые могут быть представлены композицией простейших действий типа «движение» и «превращение» элементарных частиц и составляют класс реакционно-диффузионных процессов. Движения обычно подчиняются законам диффузии или конвекции, а превращения имитируют химические реакции. Все действия представляются как случайные процессы [2], поскольку в естественных условиях нет для них специальной синхронизации. Известно несколько методов моделирования случайных процессов, происходящих в пространстве. В разных предметных областях они имеют свои названия и отличаются способами объединения элементарных действий в единый пространственно-распределённый процесс. В материаловедении и микроэлектронике — это «кинетический метод Монте-Карло» [3, 4]. В каталитической химии — это просто «метод Монте-Карло» [2, 5]. На самом деле они все являются вариантами асинхронного КА, отличающимися способами построения правил перехода.
Математические модели реакционно-диффузионных явлений, выраженные в терминах асинхронных КА [6, 7], а также исследование методов композиции КА [8] позволили обобщить накопленный опыт построения сложных КА-моделей и выделить три главных режима взаимодействия элементарных операторов, называемых подстановками :
— стохастический режим, когда каждая случайно выбранная подстановка применяется к случайно выбранной клетке, и так для всех подстановок и всех клеток;
— локальная суперпозиция, когда все подстановки в случайном порядке применяются к одной и той же случайно выбранной клетке, и так для всех клеток;
— глобальная суперпозиция, когда каждая случайно выбранная подстановка применяется ко всем случайно выбираемым клеткам, и так для всех подстановок.
Для правильного выбора способа композиции при синтезе КА-модели полезно знать, как этот выбор влияет на результат моделирования, т. е. насколько различны эволюции асинхронных КА, отличающихся только способом композиции подстановок. Есть предположение, что влияние это незаметно или очень мало. В каких случаях это предположение подтверждается — на этот вопрос далее делается попытка ответить путём проведения вычислительных экспериментов над несколькими типовыми КА-моделями.
Работа организована следующим образом. В п. 2 дана формальная постановка задачи. В п. 3, 4 и 5 приводятся результаты экспериментального моделирования для процессов движения фронта, агрегации вещества и самоорганизации соответственно. Заключение содержит сравнительный анализ полученных результатов.
1. Формальное представление асинхронных КА-моделей
Формально асинхронный КА задаётся тройкой понятий К = (А, X, ©(X)), где А — алфавит состояний клеток, т. е. множество символов или чисел, обозначающих вещества, участвующие в процессе, или условия, определяющие их свойства; X = {хи : к = 1,..., N} — множество имён клеток, которое обычно задаётся множеством координат дискретного пространства конечных размеров; в(Х) — глобальный оператор, определяющий функционирование асинхронного КА. Пара (и,х), где и Е А, х Е X, называется клеткой. Множество клеток П = {(ии, хи) : ии € А,хи € X, У к, I (к = I ^ хк = х\)} образует клеточный массив, а перечень состояний всех клеток массива ПА = (и1; и2,..., и|Х|) — глобальное состояние.
В пространстве X определены подмножества имён клеток, называемые шаблонами соседства:
ЗД = {х,х + аь ... ,х + а„_1},
где ау —вектор смещения; п = |Тп(х)|. Шаблон Тп(хи) соседства выделяет в клеточном массиве П подмножество клеток-соседей для конкретной клетки , состояния которых составляют локальную конфигурацию
£(хи) = {ио,и1,... ,и„_1},
где и0, и1;... , ип-1 — состояния клеток хи, хи + а1;... , хи + ап-1 соответственно.
Элементарная операция, изменяющая состояние клеток-соседей х Е П, имеет вид подстановки
9(х) : £(х) £'(х),
где р — вероятность применения подстановки, вычисляемая исходя из известных скоростей моделируемых элементарных действий.
Применение подстановки 9 к клетке хи сводится к замене всех или некоторых состояний ц Е £ (хи) на новые значения и Е £'(хи):
и = ¡у(ио,... ,и„_1), п = |£(хи)|, (1)
где (и0,... ,ип-1) — функция перехода.
Применение подстановки к клетке Хк € X считается успешным, если все клетки с именами из Тп(хк) содержатся в П:
|(Мо,Хк), . . . , (Ип-1,Хк + Оп-1)} С П. (2)
При этом клетка (и,х) с переменным состоянием и считается принадлежащей П, если и определено в А.
Простым далее называется такой КА, у которого глобальный оператор в(Х) = = Ф(в) состоит из одной подстановки. Асинхронный режим простого КА предполагает следующий алгоритм применения подстановки:
1) с вероятностью р = 1 /|Х| выбирается клетка (и,Хк) С П;
2) к выбранной клетке применяется подстановка в(х);
3) если условие (2) для в(х) выполнено, то состояния и3- € Б (хк) немедленно меняются на значения функций переходов (1);
4) условно принимается, что |Х| повторений п. 1 и 2 составляет одну итерацию. Сложная асинхронная КА-модель обычно задаётся не одной, а набором подстановок в = ,..., вт} и способом их композиции Ф(в). Набор подстановок соответствует набору элементарных действий в моделируемом процессе, причём так как скорости К1,... , Кт выполнения этих действий различны, то каждая подстановка вг снабжается вероятностью своего применения
_ Кг
р = К1 + К2 + ... + Кт • (3)
т
Поскольку У)рг = 1, значения вероятностей в интервале (0,1) можно расположить так, 1
(г-1 г
^Рз , Р3 3=1 3=1
подстановки вг € в далее называется такой её выбор, что случайное число п € пг.
Применение глобального оператора в(Х) к массиву П(£) изменяет его глобальное состояние П(£) ^ П(£ + 1), составляя одну итерацию. Последовательность
П(0),..., П(*),..., П(Т)
называется эволюцией КА.
Способ композиции подстановок Ф^ (в) в глобальном операторе определяет алгоритм их применения к клеткам клеточного массива П(£) в течение одной итерации. Далее изучаются три основных способа композиции: стохастический (г = Б), локальный (г = Ь) и глобальный (г = С), на которых могут строиться ещё более сложные композиционные конструкции.
Стохастическая композиция Ф^ предусматривает следующий алгоритм применения подстановок:
1) с вероятностью р = 1/|Х| выбирается клетка (и,хк) € П;
2) случайно выбирается подстановка вг, которая сразу применяется к хк;
3) условно принимается, что |Х| • т повторений п. 1 и 2 составляют одну итерацию. Из приведённого алгоритма следует, что при стохастической композиции каждая случайно выбранная подстановка вг € в в течение одной итерации выполняется рг • |Х| раз, всякий раз для случайно выбранной клетки хк € Х.
Локальная суперпозиция Фь предусматривает следующий алгоритм применения подстановок:
Случайным выбором
1) с вероятностью р = 1 /IX| выбирается клетка (и,хи) Е П;
2) случайно выбирается подстановка 9г Е в и применяется к клетке х&, затем выбирается следующая подстановка 9у Е в, которая применяется к той же самой клетке, и так т раз;
3) условно принимается, что IX| повторений п. 1 и 2 составляют одну итерацию.
Из приведённого алгоритма следует, что при локальной суперпозиции каждая случайно выбранная подстановка 9г Е в в течение одной итерации выполняется рг•IX| раз, как и при стохастической композиции. Однако отличие состоит в том, что здесь все подстановки применяются к одной и той же случайно выбранной клетке.
Глобальная суперпозиция Фс предусматривает следующий порядок применения подстановок:
1) случайно выбирается подстановка 9г Е в;
2) 9г применяется последовательно к случайно выбираемым с вероятностью р = = 1 /IX| клеткам (и,хи) Е П;
3) условно принимается, что т повторений п. 1 и 2 составляют одну итерацию.
Из приведённого алгоритма следует, что при глобальной суперпозиции каждая случайно выбранная подстановка 9г Е в в течение одной итерации выполняется рг•IX| раз, так же, как в обоих предыдущих случаях. Однако отличие состоит в том, что здесь одна и та же случайно выбранная подстановка применяется ко всем случайно выбираемым клеткам.
Три приведённых способа композиции подстановок в асинхронных КА-моделях различаются степенью вводимого в алгоритмы порядка, который можно оценивать вероятностью ^ каждого применения подстановки в течение одной итерации. Эта вероятность равна обратной величине числа возможностей. Легко подсчитать, что для каждого способа композиции она равна
_ 2(0 - 2)! т!(0 - т)!
= о ' ^ = = о '
где 0 = IXI • т.
Сравнение эволюций КА-моделей при трёх приведённых способах композиции локальных операторов выполнялось для трёх типовых реакционно-диффузионных процессов: распространение фронта, агрегация вещества и самоорганизация в системе «хищник — жертва». Вычислительные эксперименты были поставлены следующим образом:
1) для каждого процесса были выбраны одна или несколько основных величин, характеризующих его эволюции;
2) каждый процесс моделировался трижды; КА-модели отличались только способом композиции подстановок. На каждой итерации выводились значения выбранных характеристик;
3) исследовались различия полученных значений для трёх способов композиции.
2. Распространение фронта двумерной волны
Первая математическая модель процесса распространения фронта волны исследована в [10, 11]. В [10] процесс называется «диффузией, соединенной с возрастанием количества вещества». Результаты оказались необычайно плодотворными не только в части математического анализа этого нелинейного явления, но и в части применения их для моделирования подобных процессов, например распространения огня,
распространения сорняков и эпидемий, химических превращений веществ в активных средах. КА-модели процесса распространения фронта также предлагались и изучались для разных применений [12]. Для современного моделирования процесс «распространение фронта» может рассматриваться как типовая составляющая более сложных явлений [13] и, следовательно, должен быть изучен подробно.
Кинетика процесса представлена как множество частиц, распределённых по пространству с неравномерной плотностью. Частицы диффундируют из области с большей плотностью на свободные места. Однако плотность в каждой точке изменяется не только в результате диффузионного перемещения, но также из-за реакции среды на эти изменения. Реакция описывается логистической функцией, имеющей вид полинома п-го порядка. В простейшем случае это функция следующего вида:
^(и,х) = а(и(х))(1 - (и(х))), (4)
где (и(ж)) —плотность вещества в клетке ж, равная осреднённому значению состояния по окрестности осреднения:
Е и(Х)
/ ( \\ ХеАу(х) ГкЛ
(и(х)) =-тт~\-• (5)
|Ау|
Здесь Ау(х) —окрестность осреднения, состоящая их множества клеток, отдалённых от ж не более чем на заданное /.
Асинхронный КА, моделирующий распространение фронта волны в двумерном пространстве Н = (А,Х, в(Х)), в своем классическом виде работает с булевым алфавитом А = {0,1}, множеством имён клеток X = {(г,3) : г,3 = 0,... , N} и глобальным оператором 0(Х) = Фх(6^(х), 6г(х)), где Ф^ —один из способов композиции Ф^, Ф^ или Ф^; —локальный оператор диффузии, 6Г — локальный оператор реакции, которые могут быть выражены одноимёнными элементарными подстановками.
Подстановка 6^(х) определена на пятиточечном шаблоне
ЗД,3) = {(г,3), (г,3 + 1), (г + 1,3), (г,3 - 1), (г - 1,3)} (6)
и имеет вид
: {u0, Ui, U2, U3, U4} —^ {u0, uiu2, u3, u4} (7)
Применение вd к случайно выбранной клетке (и0,ж^) £ П состоит в замене её состояния u0 на состояние « одной из соседних клеток, выбранных равновероятно, т. е.
(«0 = ul) & (ul = u0), если (/ — 1)/4 < rand1 < //4 & rand2 < pd, / = 1, 2, 3, 4,
где rand1, rand2 — случайные числа в интервале (0,1).
Начальным глобальным состоянием выбран клеточный массив размера 701 х 701. В центре массива находится круг радиуса R = 30, в котором все клетки имеют состояния v = 1, в остальных клетках массива состояния v = 0. При функционировании КА круг расширяется, постепенно заполняя единицами всю область. Вычислительный эксперимент состоял из запуска параллельной программы моделирования эволюций трёх КА, отличающихся способами композиции глобального оператора. Каждому способу композиции Ф z, z £ {Ф^, Ф ь, Фс}, отводился один OpenMP-поток на четырёхядерном компьютере Intel Core i7 920. В эксперименте выполнялось по 100 прогонов с разными начальными значениями генератора случайных чисел и последующим вычислением математического ожидания U(x) результирующих состояний. Для каждого способа
композиции на ¿-й итерации выводилось по два числа, характеризующих моделируемый процесс:
Г (¿) = Е и(х), V (¿) = (-Щ) - -Г ( - 1}) ,
хех 4 7
где Гг — суммарное состояние; Уг — скорость распространения фронта. Вычислительные эксперименты различались соотношением вероятностей диффузии и реакции (ра;Рг), которые равнялись (0,1; 0,9), (0,2; 0,8), (0,5; 0,5), (0,9; 0,1), (0,95; 0,05). Полученные значения Гг (¿) и Уг (¿) (рис. 1) показывают, что все три способа композиции глобального оператора на качественном уровне неразличимы, так как разность между суммарными состояниями Г?(¿) — Гь(£), Г?(¿) — Гс(£) ни в одном эксперименте ни при одном £ не превышает 3%. На размер разницы не влияет также соотношение вероятностей (ра; рг), от которого зависит скорость распространения фронта (рис. 2).
Рис. 1. Результаты вычислительного эксперимента при ра = рг = 0,5: а — суммарные состояния в зависимости от времени; б — скорости распространения фронта волны от времени
Рис. 2. Зависимость скорости распространения фронта волны V от значения вероятности применения подстановки диффузии ра. Значения V? ^ь, VG для всех случаев совпадают
3. Агрегация, ограниченная диффузией
Явление, которое называется «агрегация, ограниченная диффузией» (diffusion limited aggregation, DLM), далее —просто агрегация, не может быть представлено в виде дифференциального уравнения. Его первая модель была описана в виде взаимодействий и движений частиц в дискретном пространстве и сразу реализована на компьютере [14]. В [14] модель не была причислена ни к какому направлению математического моделирования, хотя по сути является асинхронным КА. Модель агрегации является абстракцией таких процессов, как электрогальванизация [15], образование кристаллических структур [16], рост городов [17] и др. Наиболее известны асинхронные КА-модели, в которых диффузионная составляющая описывает случайное блуждание частиц в пространстве, а реакционная составляющая преобразует блуждающую частицу в неподвижную, если она оказывается достаточно близко к другой неподвижной (явление «прилипания»).
Процессу соответствует эволюция асинхронного КА N = (A,X, в(Х)), где A = = {0,1, b}, X = {(i,j)}. Локальный оператор 0(i, j) = (9d(i,j)А(i,j)), где 9d(i,j) — подстановка наивной диффузии (7), 9r(i,j) —подстановка прилипания. Обе подстановки вероятностные, значения вероятностей pd, pr определяются по (3), исходя из известных значений коэффициентов диффузии и прилипания. Подстановка прилипания построена на шаблоне T5 (6) и имеет тот же вид, что и 9d(i,j), отличаясь только функцией перехода, которая здесь равна
, \ b, если (u0 = 1) & (ul = b) & (rand < pr), , _
Uo — \ l — 1, 2, 3, 4.
0 I Uo в остальных случаях,
Исходное глобальное состояние П(0) —равномерное с плотностью 0,5 распределение «единиц» по всему массиву, кроме пятна из клеток-зародышей в состоянии v(i,j) = b размера 2 х 4, расположенного в середине нижней границы массива. В процессе эволюции вокруг зародышей образуется растущая древовидная структура (рис. 3), похожая на коралл или мох.
Рис. 3. Древовидные структуры, построенные асинхронным КА агрегации при стохастической композиции подстановок 0d(i, j) и 9r(i, j): а — при (Pd; Pr) = (0,5; 0,5) и t = 520; б — при (pd; pr) = (0,9; 0,1) и t = 2500
Вычислительные эксперименты проводились для клеточного массива размера 500 х 200. Промоделированы три случая, отличающиеся соотношениями вероятностей (pd;pr) = (0,5; 0,5), (0,8; 0,2) и (0,95; 0,05). Во всех случаях каждому способу композиции Фг, z £ {S,L,G}, отводился один OpenMP-поток на четырёхъядерном компьютере Intel Core i7 920. В эксперименте выполнялось по 100 прогонов с разными начальными значениями генератора случайных чисел и последующим вычислением
математического ожидания состояний и(х). Для каждого способа композиции выводились зависимости Дz (¿) для г Е {£, Ь, С} (рис. 4) и вычислялись значения 8г (¿) при двух значениях радиуса Я =30 и 60:
Дг(¿) = £ и(х), (¿) = log(Nz(Я))/log(пЯ2),
хех
где Дг (¿) — суммарная масса на ¿-й итерации; 8г — фрактальная размерность построенной структуры; Nz (Я) — количество клеток в состоянии Ь, находящихся на расстоянии не более чем Я от зародыша (таблица).
а бе
Рис. 4. Зависимости суммарных масс Д^(£), Да(Ь) и Д^(¿) для трёх случаев соотношения интенсивностей диффузии и прилипания: а — (Ра;Рг) = (0,5; 0,5); б — (Ра; Рг) = (0,8; 0,2); в — (рА;Рг) = (0,95; 0,05)
Фрактальные размерности древовидных структур, построенных в ходе моделирования процесса агрегации
К 30 60
Рг) (0,5; 0,5) (0,95; 0,05) (0,5; 0,5) (0,95; 0,05)
Фь 1,65621 1,64746 1,65243 1,69751
1,65953 1,64608 1,6517 1,70208
1,64146 1,6474 1,65126 1,69346
Полученные результаты показывают, что кривые суммарных масс для локальной Ф^, глобальной Ф^ и стохастической Ф^ композиций для каждого случая (ра;рг) имеют одинаковый характер, но различаются при разных соотношениях интенсивности диффузии и прилипании. Так, при (ра; рг) = (0,5; 0,5) и (0,8; 0,2) скорость роста суммарной массы при локальной композиции больше, чем при стохастической, примерно на 5%, тогда как при (ра; рг) = (0,95; 0,05) разницы в скоростях роста суммарной массы незаметно.
Различия в значениях фрактальных размерностей 8Х (¿) от способа композиции не заметны (< 1%). То же самое можно сказать и по поводу пространственного образа древовидной структуры (см. рис. 3).
Общий вывод из проведённого вычислительного эксперимента состоит в том, что свободный выбор способа композиции локальных операторов допустим, если нет сильных ограничений по точности. Кроме того, всегда возможна корректировка масштабов при переходе от модельных к физическим шагам по времени и пространству. Важно, что инвариант КА одинаков при разных способах композиции.
4. Самоорганизация в системе «хищник — жертва»
Исследование взаимоотношений хищников с жертвами [18], обитающих совместно в замкнутом пространстве, является типовой задачей самоорганизации [19], которая привлекает внимание не только экологов, но и математиков, физиков, химиков и даже социологов. В традиционной математике системы «хищник — жертва» описываются дифференциальными уравнениями в частных производных с полиномиальными нелинейными функциями. Поведенческие свойства изучаются при помощи качественной теории нелинейных систем. Компьютерная имитация выполняется путём решения систем нелинейных дифференциальных уравнений [18]. Возможности этих методов ограничены, во-первых, размерностью фазового пространства, во-вторых, трудностями решения нелинейных уравнений на суперкомпьютерах. Клеточно-автоматные модели систем «хищник — жертва» изучались также в синхронном [20] и асинхронном вариантах [21].
Взаимодействия между хищниками и жертвами происходят следующим образом. Хищники (щуки, лисы) питаются жертвами (планктоном, зайцами). Если в каком-то месте пищи достаточно (плотность жертв больше, чем плотность хищников), то плотность хищников увеличивается: хищники размножаются с вероятностью, пропорциональной их плотности. Если пищи не хватает (плотность жертв меньше, чем плотность хищников), то плотность хищников уменьшается: они умирают от голода с вероятностью, пропорциональной недостатку пищи. Жертвы всегда пытаются размножаться с вероятностью, пропорциональной логистической функции [18]. И хищники, и жертвы перемещаются в пространстве по закону диффузии, причём хищники более подвижны, чем жертвы.
Асинхронный КА, моделирующий систему из одного вида хищников и одного вида жертв, описывается параллельной композицией из двух диффузионно-реакционных КА = (А,Хь, 6ь(X)) и Ни = (А,Хи, 6и(Х)), где А — булев алфавит; Хь = {(г,3)ь} и Хи = {(г,3)и} —два изоморфных друг другу множества, определяющих клеточные массивы и Пи, на которых отображаются эволюции хищника и жертвы соответственно. Пара клеток, связанных отношением изоморфизма, обозначается (г,3) (без индексов) и называется парой клеток-близнецов.
Глобальный оператор 6(Х) является композицией следующих подстановок:
— ^¿ь((г,3)ь) —диффузия хищника по массиву ;
— $аи((г,3)и) —диффузия жертв по массиву ;
— ((г,3)ь) —реакция (изменение плотности) хищников;
— 0ги((г,3)и) —реакция жертв.
Подстановки наивной диффузии (7) ((г, 3)ь) и #аи(г,3)и) применяются каждая на своём клеточном массиве, различаясь только коэффициентами диффузии ^ Ди).
Подстановки реакций построены на шаблоне Т2 = {(г,3)ь, (г,3)и}, содержащем пары близнецов:
((г,3)ь): ((V), («))
М(г,3)и) : ((«), (V))
Рг , ( / а —> (V',«'),
рг ,11 /\
где (V) и (и) — осреднённые по окрестности Ау = {(г + а, 3 + Ь) : а, Ь = —г,..., 0,... , г} величины, имеющие смысл плотности (5) организмов в точках, соответствующих координатам клеток (г,3)ь и (г,3)и.
Новые состояния в (8) зависят от плотностей обоих организмов следующим образом:
, \ 0, если (и) > (v) & rand < pv^0, vn = ^ 1 / \ . / \ n l. (У)
1, если (v) > (и) & rand
где вероятности pv^0 и pv^ определяются на основе следующих соображений. Если (v) > (и), то избыток хищников ((v) — (и)) на площади Av((i, j)v) не имеет пищи и умирает. Следовательно, ((v) — (и)) клеток из T((i,j)v) должны сменить состояния v = 1 на v = 0. Если (и) > (v), то хищнику достаточно пищи, он размножается в соответствии с логистической функцией Fv ((v)) вида (4). Отсюда вероятности в (9) равны
pv^0 = (v) — (и), если (v) > (и), Pv^i = Fv ((v)), если (и) > (v).
Функция перехода 9ru выглядит так же, как (8), отличаясь только вероятностями pu^0 и pu^, значения которых определяются на основе следующих соображений. Если (и) < (v) (жертва поедается хищником), её плотность уменьшается пропорционально плотности хищника. Если же (и) > (v), то жертва размножается с вероятностью, пропорциональной ((и) — (v)):
Pu^o = (v), если (v) > (и),
Pu^i = F((и) — (v)), если (и) > (v).
Вычислительные эксперименты проводились для клеточного массива П(0) = = Qv(0) U Пи(0) размера |Xv| = |Xu| = 400 х 800 (рис.5). В начальном состоянии жертвы распределены по пространству равномерно: У(и, (i,j)u) Е Пи(0) ((и) = 0,4). Хищник же собран в плотную полосу в середине области:
,, Г 0,6 для всех (i,j )v : i = 0,... , 399, j = 369,..., 439, [ 0,2 в остальных (i, j)v Е Hu(0).
Константы скоростей для диффузии и реакций приняты равными Kdv = 10, Kdu = 0,1, Krv = 0,2, Kru = 0,2, что соответствует следующим вероятностям (3) применения локальных операторов 6dv, 0du, , $ru:
Pdv = 0,952; pdu = 0,00952; prv = 0,01904; pra = 0,01904.
Рис. 5. Начальное состояние 0(0) КА, моделирующего самоорганизацию в системе «хищник — жертва» при появлении «косяка» хищника в виде полосы шириной 60 клеток и плотностью (г) = 0,6 при исходной плотности хищника (г) = 0,2 (а) и исходной плотности жертвы (и) = 0,4 (б)
Моделирование проводилось для трёх исследуемых способов композиции подстановок, которые для случая двух организмов и двух клеточных массивов выполняются по следующим алгоритмам.
1. Стохастическая композиция Ф^. Случайно выбирается клетка (г, 3) € X. В соответствии с распределением вероятностей из четырёх подстановок 9аь, 9аи, 9гь, 9ги выбирается одна, которая применяется к соответствующей клетке-близнецу ((г,3)ь или (г,3)и). Итерация содержит 4|Х| таких шагов.
2. Локальная композиция Ф^. Случайно выбирается клетка (г, 3) € X. Все четыре подстановки применяются к этой клетке последовательно в случайно выбранном порядке, каждая — к соответствующему близнецу. Итерация содержит Х таких шагов.
3. Глобальная композиция Ф^. Случайно выбирается одна из подстановок и применяется ко всем случайно выбираемым клеткам соответствующего ей подмассива. Итерация содержит 4 таких шага.
В процессе эволюции происходит постепенное перераспределение хищников и жертв по пространству, в результате которого достигается устойчивое равновесие в виде скоплений хищников и жертв, имеющих форму круглых или овальных пятен (рис. 6).
Рис. 6. Глобальные устойчивые состояния хищника (a) и жертвы (б) при стохастической композиции, полученные в результате моделирования
Каждому способу композиции Фг, z Е {S,L,G}, отводился один OpenMP-поток на четырёхъядерном компьютере Intel Core i7 920. В эксперименте выполнялось по 100 прогонов с разными начальными значениями генератора случайных чисел и последующим вычислением математического ожидания u(x) и v(x). На каждой итерации в файл выводились суммарные массы хищников и жертв:
(t)= Е (v, (i,j)v), лг(t) = E (u,(i,j)v).
(i,j)v ex„ (i,j)u
Сравнение эволюций для трёх способов композиции проводилось по суммарным массам Sz (t) хищников и Лх (t) жертв (см. рис. 7) , а также по количеству элементов связности (плотных скоплений) oz и Az в устойчивом состоянии. Сравнение по суммарным массам показывает максимальную разницу не более 3 %. Количество скоплений в устойчивых состояниях везде равно 23.
Заключение
Проведено исследование способов построения глобального оператора асинхронного КА, моделирующего реакционно-диффузионный процесс. Цель исследования — ответить на вопрос, насколько отличаются эволюции асинхронного КА при трёх способах композиции подстановок при построении глобального оператора: стохастическом, локальном и глобальном. Исследование проведено методом вычислительного
Рис. 7. Зависимости суммарных масс хищников (а) и жертв (б) от времени, полученные в результате экспериментального моделирования
моделирования для трёх типовых реакционно-диффузионных процессов: распространение фронта, агрегация вещества из раствора, самоорганизация в системе «хищник — жертва». Возможность выбора способа композиции полезна при синтезе сложных КА-моделей, так как позволяет разработчику согласовать вычислительный процесс с архитектурой суперкомпьютера, а также рационально встраивать его в систему более сложных моделей.
Главный вывод из проведённого исследования состоит в том, что результаты моделирования при выбранных способах построения глобального оператора качественно совпадают. Более того, количественные различия заметны (> 5 %) только для процесса агрегации. Следует отметить, что в случае моделирования естественных процессов качественные совпадения часто вполне достаточны, а расхождения количественных значений получаемых характеристик до 10 % считаются допустимыми. В тех особых случаях, когда полученные количественные расхождения недопустимы, можно внести коррекцию в масштабы шагов по времени или пространству. Качественная идентичность эволюций позволяет это сделать, хотя потребует дополнительных исследований.
ЛИТЕРАТУРА
1. Фон Нейман Дж. Теория самовоспроизводящихся автоматов. М.: Мир, 1971. 384 c.
2. Metropolis N. and Ulam S. The Monte Carlo method // Amer. Statist. Assoc. 1949. V. 44. No. 247. P. 335-341.
3. Chatterjee A. and Vlaches D. G. An overview of spatial microscopic and accelerated kinetic Monte-Carlo methods // J. Computer-Aided Mater. Des. 2007. V. 14. P. 253-308.
4. Nurminen L., Kuonen A., and Kaski K. Kinetic Monte-Carlo simulation on patterned substrates // Phys. Rev. B63, 035407. 29 December 2000.
5. Matveev A. V., Latkin E. I., Elokhin V. I, and Gorodetskii V. V. Turbulent and stripes wave patterns caused by limited C0ads diffusion during CO oxidation over Pd(110) surface: kinetic Monte Carlo studies // Chem. Eng. J. 2005. V. 107. P. 181-189.
6. Бандман О. Л. Дискретное моделирование физико-химических процессов // Прикладная дискретная математика. 2009. №3. C. 33-49.
7. Kireeva A. Parallel implementation of totalistic cellular automata model of stable patterns formation // 12th Int. Conf. "Parallel Computing Technologies", St.-Petersbourg, 2013. LNCS. 2013. V. 7979. P. 347-360.
8. Бандман O. Л. Методы композиции клеточных автоматов для моделирования пространственной динамики // Вестник Томского государственного университета. 2004. Приложение №9(1). С. 183-193.
9. Achasova S., Bandman O, Markova V., and Piskunov S. Parallel Substitution Algorithm. Theory and Application. Singapore: World Scientific, 1994. 180 p.
10. Колмогоров А. Н., Петровский И. Г., Пискунов И. С. Исследование уравнения диффузии, соединенной с возрастанием количества вещества, и его применение к одной биологической проблеме // Бюллетень МГУ. Секция А. 1937. Вып. 6. С. 1-25.
11. Fisher R. A. The genetical theory of natural selection. Oxford: Univ. Press, 1930. 58 p.
12. Szakaly T., Lagzi I., Izsak F., et al. Stochastic cellular automata modelling excitable systems // Central Eur. J. Phys. 2007. No. 5(4). P. 471-486.
13. Van Saarloos W. Front propagation into unstable states // Phys. Rep. 2003. V. 386. P. 29-222.
14. Witten T. A., Jr. and Sander L. M. Diffusion-limited aggregation, a kinetic critical phenomenon // Phys. Rev. Lett. 1981. V. 47. No. 19. P. 1400-1403.
15. Ackland G. J. and Tweedie E. S. Microscopic model of diffusion limited aggregation and electrodeposition in the presence of leveling molecules // Phys. Rev. E73, 011606. 26 January 2006.
16. Bogoyavlenskiy V. A. and Chernova N. A. Diffusion-limited aggregation: a relationship between surface thermodynamics and crystal morphology // Phys. Rev. E. 2000. No. 61(2). P. 1629-1633.
17. Batty M. and Longley P. Urban growth and form: scaling, fractal geometry, and diffusion-limited aggregation // Environment and Planning. 1989. V. 21(11). P. 1447-1472.
18. Свирежев Ю. М. Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987. 368 с.
19. Nicolis G. and Prigogine I. Self-organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations. N.Y.: Wiley, 1977.
20. Kutson J. D. A survey of the use of cellular automata and cellular automata like models for simulating a population of biological cells. Iowa State University, 2011. Graduate Thesis and Dissertations. Paper 10133. http://lib.dr.iastate.edu/efd
21. Chen Q., Mao J., and Li W. Stability analysis of harvesting strategies in a cellular automata based predator - prey model // Cellular Automata. LNCS. 2006. V.4173. P. 268-376.
REFERENCES
1. Von Neumann J. Teoriya samovosproizvodjashhihsja avtomatov. Moscow, Mir Publ., 1971. 384 p. (in Russian)
2. Metropolis N. and Ulam S. The Monte Carlo method. Amer. Statistical Assoc., 1949, vol.44, no.247, pp. 335-341.
3. Chatterjee A. and Vlaches D. G. An overview of spatial microscopic and accelerated kinetic Monte-Carlo methods. J. Computer-Aided Mater. Des., 2007, vol.14, pp. 253-308.
4. Nurminen L., Kuonen A., and Kaski K. Kinetic Monte-Carlo simulation on patterned substrates. Phys. Rev. B63, 035407. 29 December 2000.
5. Matveev A. V., Latkin E. I., Elokhin V. I, and Gorodetskii V. V. Turbulent and stripes wave patterns caused by limited C0ads diffusion during CO oxidation over Pd(110) surface: kinetic Monte Carlo studies. Chem. Eng. J., 2005, vol. 107, pp. 181-189.
6. Bandman O. L. Diskretnoe modelirovanie fiziko-himicheskih processov. Prikladnaya Diskretnaya Matematika, 2009, no. 3, pp. 33-49. (in Russian)
7. Kireeva A. Parallel implementation of totalistic cellular automata model of stable patterns formation. 12th Int. Conf. "Parallel Computing Technologies", St.-Petersbourg, 2013. LNCS, 2013, vol. 7979, pp. 347-360.
8. Bandman O. L. Metody kompozicii kletochnyh avtomatov dlja modelirovanija prostranstvennoj dinamiki. Vestnik Tomskogo gosudarstvennogo universiteta, 2004, Prilozhenie no. 9(1), pp. 183-193. (in Russian)
9. Achasova S., Bandman O, Markova V., and Piskunov S. Parallel Substitution Algorithm. Theory and Application. Singapore, World Scientific, 1994. 180 p.
10. Kolmogorov A. N., Petrovskij I. G., Piskunov I. S. Issledovanie uravnenija diffuzii, soedinennoj s vozrastaniem kolichestva veshhestva, i ego primenenie k odnoj biologicheskoj probleme. MSU Bull., Sec. A, 1937, no. 6, pp. 1-25. (in Russian)
11. Fisher R. A. The genetical theory of natural selection. Oxford, Univ. Press, 1930. 58 p.
12. Szakaly T., Lagzi I., Izsak F., et al. Stochastic cellular automata modelling excitable systems. Central Eur. J. Phys., 2007, no. 5(4), pp. 471-486.
13. Van Saarloos W. Front propagation into unstable states. Phys. Rep., 2003, vol.386, pp.29-222.
14. Witten T. A., Jr. and Sander L. M. Diffusion-limited aggregation, a kinetic critical phenomenon. Phys. Rev. Lett., 1981, vol.47, no. 19, pp. 1400-1403.
15. Ackland G. J. and Tweedie E. S. Microscopic model of diffusion limited aggregation and electrodeposition in the presence of leveling molecules. Phys. Rev. E73, 011606. 26 January 2006.
16. Bogoyavlenskiy V. A. and Chernova N. A. Diffusion-limited aggregation: a relationship between surface thermodynamics and crystal morphology. Phys. Rev. E. 2000, no. 61(2), pp. 1629-1633.
17. Batty M. and Longley P. Urban growth and form: scaling, fractal geometry, and diffusion-limited aggregation. Environment and Planning, 1989, vol. 21(11), pp. 1447-1472.
18. Svirezhev Ju. M. Nelinejnye volny, dissipativnye struktury i katastrofy v jekologii. Moscow, Nauka Publ., 1987. 368 p. (in Russian)
19. Nicolis G. and Prigogine I. Self-organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations. N.Y., Wiley, 1977.
20. Kutson J. D. A survey of the use of cellular automata and cellular automata like models for simulating a population of biological cells. Iowa State University, 2011. Graduate Thesis and Dissertations. Paper 10133. http://lib.dr.iastate.edu/efd
21. Chen Q., Mao J., and Li W. Stability analysis of harvesting strategies in a cellular automata based predator — prey model // Cellular Automata, LNCS, 2006. vol.4173, pp.268-376.
