CC BY
28
УДК 548.55-522:004.942-021
Г. Ю. Аверчук1, Е. С. Куркина2
СРАВНИТЕЛЬНЫЙ АНАЛИЗ МЕТОДОВ МОНТЕ-КАРЛО НА ПРИМЕРЕ РАСЧЕТА СЛОЖНОЙ ДИНАМИКИ РЕШЕТОЧНОЙ МОДЕЛИ ХИМИЧЕСКОЙ РЕАКЦИИ*
Методами Монте-Карло рассчитывается сложная колебательная динамика микроскопической модели химической реакции, происходящей на поверхности катализатора. Используется пять наиболее употребительных алгоритмов и оценивается сложность и производительность каждого из них. Показано, что время расчета может отличаться более чем на два порядка. Выявлен самый эффективный алгоритм Монте-Карло. С помощью него исследованы механизмы возникновения колебаний и волн на микроуровне в решеточной модели, содержащей миллионы узлов.
Ключевые слова: микроскопическая стохастическая модель, гетерогенно-каталитическая реакция, методы Монте-Карло, автоколебания.
1. Введение. Микроскопические имитационные модели являются наиболее полными моделями, описывающими динамику систем с вероятностной природой взаимодействий, таких, как химические реакции. С помощью этих моделей можно детально описать структуру катализатора и сложные взаимодействия веществ с катализатором и между собой в адсорбционном слое, исследовать механизмы возникновения макроскопических структур и волн на микроуровне и изучить влияние флуктуаций. В основе микроскопического описания лежит решеточная модель с определенным дискретным множеством микросостояний и заданным набором элементарных процессов, реализация которых приводит к изменению макросостояния системы. Эволюция системы рассматривается в приближении марковского случайного процесса с непрерывным временем и конечным фазовым множеством, элементы которого представляют собой микросостояния системы. Эволюция рассчитывается с помощью динамических стохастических алгоритмов, называемых методами Монте-Карло (МК); при этом получаются отдельные наиболее вероятные траектории марковского случайного процесса.
Основным препятствием для использования микроскопических стохастических моделей являются их огромная размерность, а значит, и запрашиваемые вычислительные ресурсы для реализации методом МК. Поэтому для расширения возможностей применения имитационных моделей актуально выявление наиболее эффективных алгоритмов. На примере расчетов сложной динамики микроскопической модели химической реакции, происходящей на поверхности катализатора, производится сравнение пяти наиболее употребительных алгоритмов МК. Оценивается сложность и производительность каждого из них. Выявлен самый быстрый алгоритм.
2. Микроскопическая стохастическая модель реакции. Рассмотрим микроскопическую стохастическую модель реакции N0 + СО, происходящей на поверхности платинового катализатора с правильной квадратной решеткой (грань Pt(100)), изображенной на рис. 1. Кинетическая схема реакции состоит из семи элементарных стадий, описывающих элементарные реакции с двумя адсорбционными частицами N0 и СО:
1) (Х< >),,,, + * ^ [NO] — адсорбция N0;
2) (СО)Гаэ + * ^ [СО] — адсорбция СО;
3) [СО] + [NO] % (С02)гаэ + (N)raa + 2 * — образование С02;
4) [NO] (ХО),,,, + * — десорбция N0;
1 Факультет ИКТ РХТУ, асп., e-mail: altermnQgmail.com
2 Факультет ВМК МГУ, вед. науч. сотр., д.ф.-м.н., e-mail: elena.kurkinaQcs.msu.su
* Работа выполнена при поддержке РФФИ, проекты № 11-08-00979, 11-01-00887 и в рамках государственного контракта № 11.519.11.4004.
5) [СО] % (С0)газ + * — десорбция СО;
6) [N0] + * * + [N0] — поверхностная диффузия N0;
7) [СО] + * —> * + [СО] — поверхностная диффузия СО.
п п п
1 2
□ ш и
1 2
□ □ □
п п п п
ш 1 1 ш
□ и и
ш 1 1 2
□ □ □
а б
Рис. 1. Модель квадратной решетки (первые и вторые соседи): а — для одноузель-ного процесса; б — для двухузельного процесса
Здесь * — свободный адсорбционный центр; [N0], [СО] — адсорбированные частицы; (ГТО)газ, (СО)газ, (С02)газ? (^)газ — молекулы в газовой фазе; • • • 5 — константы скоростей элементар-
ных стадий (с-1).
Адсорбционный слой поверхности катализатора представлен моделью двумерного решеточного газа размером N = N1 х N2 с периодическими граничными условиями. Совокупность всевозможных микросостояний адсорбционного слоя образует фазовое пространство системы X. Занумеруем микросостояния по порядку: Х\, Х2, . . . , Каждое состояние системы X^ складывается из состояний всех узлов решетки и характеризуется вероятностью Рхг Состояние г-го узла в момент времени £ описывается вектором заполнения = где д = 0,1,2 в зависимости от того, является ли узел свободным или занятым частицей [N0] или [СО], а состояние всей системы X^ характеризуется матрицей состояния:
= 5(Х,-) = д = 0,1,2, г = 1,...,ЛГ. (1)
Элементарным процессом, или событием, называется реализация одной из возможных элементарных реакций, описываемых кинетической схемой, приводящая к изменению состояния одного или нескольких соседних узлов. В результате поочередного осуществления элементарных событий происходит эволюция системы.
Изменение вероятностей всех возможных состояний решетки Рх{1)-> X Е X, в марковском приближении описывается основным кинетическим уравнением (ОКУ)
^Г = ^ - Р*у(х *')] (2)
аг х'
с заданным начальным распределением вероятности состояний Рх{1 = 0) = Р^- Здесь У{Х' —>► X) — интенсивность перехода системы из состояния X' в состояние X в момент Интенсивность перехода определяется суммой скоростей всевозможных элементарных процессов, выводящих систему из данного состояния, и зависит только от текущего состояния решетки.
В рассматриваемой реакции имеют место одноузельные и двухузельные процессы, происходящие на соседних узлах. Ранее было показано, что для адекватного описания сложной нелинейной динамики реакции необходимо учитывать влияние ближайших "соседей" на скорость протекания элементарного процесса в узле или паре узлов [1, 2]. Это влияние описывается с помощью энергетических параметров латеральных взаимодействий изменяющих энергию активации Е^ г-го элементарного процесса (р — сорт адсорбированной частицы, участвующей в реакции). Энергия активации в имитационной модели рассчитывается по формуле
= Ег + О + ^^г, СО
где пит — число частиц [NO] и [СО] соответственно, находящихся на расстоянии первого соседства с рассматриваемой ячейкой, состоящей из одного узла (или двух узлов в случае двухузельного процесса) (рис. 1).
Значения энергий активации и предэкспонент скоростей элементарных процессов, а также параметры eip приведены в таблице. Внешние параметры, отвечающие лабораторным экспериментам, лежат в диапазонах: /'хо ~ Ю-74"-5 мбар, Pqo ~ Ю-74"-5 мбар, /Ч.о > Pqo-, Т ~ 400 -г- 450 К.
Система ОКУ (2) имеет огромную размерность, поэтому для его решения используются приближенные методы. Отдельные траектории эволюции реакционной системы в пространстве состояний могут быть рассчитаны с помощью алгоритмов Монте-Карло, которые в западной литературе называют кинетическими ("Kinetic Monte Carlo") [3-8].
Значения энергий активации, предэкспонент скоростей элементарных процессов и параметры латеральных взаимодействий
№ Щ, с-1 Ei, ккал NO) ккал £»,СО, ккал
1 1.93 х 10ь 0 0 0
2 1.93 х 10ь 0 0 0
3 2.0 х 1015 24.5 -2.0 -0.8
4 2.0 х 101Ь 37.0 1.8 2.2
5 1.0 х 101Ь 37.5 2.2 1.0
3. Алгоритмы стохастического моделирования. Метод Монте-Карло заключается в построении марковской цепочки случайных событий [9], соответствующих решению ОКУ (2). Эта цепочка состоит из двух взаимосвязанных последовательностей — последовательности моментов времени переходов {¿о < ti < ... < tk < ■ ■ •} из одного состояния в другое и последовательности этих состояний {XQ < Xi < ... < Xk < ...}. Пуассоновский поток событий, происходящих со скоростью Vx^X'(t)-, порождает цепочку значений времени ожидания {5^}, определяющих интервалы между переходами. Пусть состояние реакционной системы в момент t характеризуется матрицей состояний S(t) =
= {«¿g} (!)• Вычислим скорости Vj всех возможных элементарных событий j I.....Л>. выводящих
систему из текущего состояния, составим матрицу скоростей W = {и/} и рассчитаем суммарную скорость
Nv
V = Y,VJ- (3)
з=i
Отметим, что при переходе к другому состоянию изменяется заполнение только одной ячейки или немногих соседних ячеек. Таким образом, изменения в матрице состояний S и матрице скоростей W носят локальный характер. Расчет эволюции системы методом МК состоит из нескольких этапов.
Первый этап состоит в задании начального состояния фрагмента решетки. Начальное состояние фрагмента S(0) либо получается в результате произведенного ранее расчета, либо готовится специальным образом (например, является результатом случайного разбрасывания частиц по узлам фрагмента). Далее осуществляется переход ко второму состоянию, а от второго к третьему и т.д. по одному и тому же алгоритму, включающему следующие основные стадии:
1) расчет времени пребывания At), системы в текущем состоянии X
2) выбор элементарного события, переводящего систему в другое состояние
3) обновление матриц состояний и скоростей.
3.1. Метод отказа (rejection). Все скорости элементарных событий, выводящие систему из текущего состояния, нумеруются и привязываются к конкретному месту на решетке, где они могут происходить [3, 10]. Далее последовательно производятся следующие действия.
1. Вычисление скоростей элементарных актов. Рассчитывается суммарная скорость V (3) и находится максимальная скорость v ^ Vj, j = 1,..., Nv.
2. Определение момента выхода системы из текущего состояния. Генерируется значение случайной величины г] € (0,1] (здесь и далее рассматриваются случайные числа [10], равномерно распре-
деленные на некотором интервале). Время пребывания системы в текущем состоянии X рассчитывается по формуле
Д< = (4)
3. Выбор события: а) выбирается случайное число £ € [О, Ж„); б) в качестве возможного элементарного события, осуществляющего выход из данного состояния, рассматривается событие под номером к: к = + 1; в) это событие осуществляется, если выполняется условие
(5)
V
Если же неравенство (5) не выполняется, то процедура выбора события повторяется начиная с пункта а, пока не будет успешной.
4- Изменение состояния фрагмента. Если случайно выбранное событие оказалось подходящим и условие (5) выполнилось, то производятся реализация этого события и переход к новому состоянию, локально изменяется матрица состояний и обновляется список скоростей, которые могут выводить систему из нового состояния. Далее происходит переход к п. 2 алгоритма.
Достоинством метода отказа является простота. Трудозатратность этого алгоритма определяется отношением общего числа попыток выбрать событие к числу успешных попыток. Обратная величина д определяет эффективность метода:
V 1 К
1' -г- V г-Х... V —
Она равна отношению средней суммы V всех скоростей элементарных событий, осуществляющих выход из текущего состояния, к верхней границе этой суммы V.
Эффективность описанного метода невелика для систем, в которых скорости элементарных событий сильно различаются (на порядки) и заметно преобладает число элементарных реакций с небольшой скоростью. Тогда алгоритм делает много безуспешных попыток, когда выбранное событие отвергается. Поэтому этот алгоритм называют методом отказа.
3.2. Метод частичных сумм (rejection-free). В базовом методе частичных сумм [4] пп. 1, 2
являются такими же, как в методе отказа, а п. 3 существенно отличается.
3. Выбор события также производится за несколько шагов: а) рассчитываются частичные сум-
то
мы скоростей Vm = ^ Vj, т = 1,..., Nv, и упорядочиваются в списке, тогда общая сумма (3) равна з=1
последней частичной сумме V = Vnv ; б) выбирается случайное число £ € [О, V); в) последовательно перебирается список частичных сумм до тех пор, пока случайное число не попадет в интервал Vm-i ^ £ < Vm; г) выбирается и осуществляется т-е элементарное событие.
4- Изменение состояния фрагмента. Далее производится обновление матрицы состояния, матрицы скоростей и переход к п. 2.
С точки зрения метода отказа в этом алгоритме каждая попытка выбрать элементарное событие является успешной. Однако метод довольно трудоемкий, поскольку просмотр списка частичных сумм и сравнение их со случайным числом требуют проведения большого числа операций с плавающей точкой, порядка 0(т) или 0(NV) как оценки сверху.
Затраты можно существенно уменьшить, если упорядочить список скоростей, разбив их на М групп и подсчитав сумму скоростей в этих группах: V\, V2, ■ ■ ■, V д /.
3.3. Динамический метод Монте-Карло. Объединим в группу одинаковые элементарные события. Число групп М обычно намного меньше, чем общее число скоростей в списке (М Пункты 1, 2 являются такими же, как в методе отказа.
3. Выбор группы событий: а) вычисляются суммы скоростей в каждой группе: V; = т? • VI, где гп1 — число элементарных событий в группе, осуществляющих 1-й элементарный процесс; далее
то
находятся частичные суммы Ут = ^ Т^, т = 1,..., Л/. и упорядочиваются в списке. Тогда суммарная
1=1
скорость равна V = Ум', б) выбирается случайное число £ € [О, V); в) последовательно перебирается список частичных сумм до тех пор, пока случайное число не попадет в интервал Ут-1 ^ £ < Ут; г) выбирается т-я группа, состоящая из тгц одинаковых элементарных событий.
4- Выбор события: а) выбирается случайное число С € [0, тг); б) в качестве элементарного события, осуществляющего выход из данного состояния, выбирается событие под номером к: к = +1.
5. Изменение состояния фрагмента. Производится обновление матрицы состояния и скоростей, осуществляется переход к п. 1.
На решетках большого размера динамический метод [11] позволяет быстрее локализовать элементарное событие, осуществляющее выход из текущего состояния, чем предыдущий алгоритм. Число действий, затрачиваемых на определение элементарного процесса, здесь намного меньше, чем в методе частичных сумм, и составляет О(М). И число операций, необходимых для выбора элементарного события в группе, также имеет порядок О(М), поскольку случайно выбирается одно из равновероятных событий.
3.4. Кинетический метод Монте-Карло. Этот метод часто используется для моделирования роста кристаллов, пленок, фуллеренов и других структур [6]. Время выхода системы из текущего состояния здесь рассчитывается по отличному от всех вышеописанных алгоритмов, но тоже соответствует течению реального физического времени [5]. Алгоритм состоит из следующих шагов.
1. Задание начального состояния системы.
2. Определение момента и места (элементарной ячейки) выхода системы из текущего состояния. Для этого: а) в каждой ;/'-й элементарной ячейке ] = 1,..., N рассчитывается сумма У^ скоростей элементарных событий, которые могут произойти и изменить состояние данной ячейки, а значит, и всей системы; б) для каждой ячейки рассчитывается время пребывания ее в текущем состоянии и момент выхода из данного состояния благодаря пуассоновскому потоку элементарных событий, связанных с данной ячейкой. Для этого генерируется значение случайной величины r]j € (0,1]. Время пребывания ячейки в текущем состоянии рассчитывается по формуле
= .У I.....(6)
уз
в) производится сравнение значений времени пребывания (6) ячеек в текущих состояниях и выбирается ячейка ]*, имеющая минимальное время задержки ^ А^-.
3. Выбор события: а) генерируется случайное число £ € [0, V}*); б) выбирается элементарное событие, осуществляющее выход ячейки ]* и всей системы из текущего состояния.
4- Изменение состояния фрагмента. Производится обновление состояний ячейки и соседних ячеек, которые затронул элементарный процесс. Только в этих ячейках заново считаются скорости возможных элементарных процессов и находятся их суммы. Далее осуществляется переход к п. 2, б.
Во многих задачах кинетический метод работает существенно быстрее динамического метода, поскольку не требует полного обновления матрицы скоростей событий; изменения имеют локальный характер и связаны только с одной ячейкой. Вычислительная сложность алгоритма обусловливается необходимостью генерировать N раз случайное число и N раз вычислять время задержки (6), что в принципе ускоряется распараллеливанием. Таким образом, чтобы найти ячейку и выбрать элементарное событие, нужно сделать О (Ж) действий. В свою очередь обновление матриц состояний и скоростей обойдется всего в 0{Ые х [Ып + 1)), где Ые — число элементарных стадий в кинетической схеме, а1„ — число первых соседей (рис. 1).
3.5. Многоуровневые методы. Рассмотрим К-уровневый метод, в котором решетка имеет ¡у = дК уЗЛ0В1 Разделим фрагмент на д частей, каждую часть также разделим на д частей [7]. Процесс деления будем продолжать до тех пор, пока на последнем К-м делении не получим отдельные ячейки. Структура списка скоростей в виде К-уровневого дерева показана на рис. 2. Максимально возможное число уровней дерева равно К = -/V, а дерево в таком случае является бинарным. Далее рассчитываются частичные и полные суммы скоростей на каждой ветви:
то
1) на нижнем К-м уровне вычисляются У^т = ^ Я!з-. т = ■ ■ ■ > г-, — частичные суммы,
з=1
состоящие из сумм скоростей элементарных процессов г-й ветки, % = 1,...,д, и общие суммы всех скоростей V '/4;
гтНп nVi nVi nVi п*п nVi nVi nVi nVi nVi nVi гтНп n*n nVi n*n nVi
Рис. 2. Структура данных "дерево" в многоуровневом методе с К = 3 и g = 4
2) на каждой ветви следующего уровня также вычисляются частичные суммы и полная сумма
га
скоростей = Y^ ^ = 1,..., g, г = 1,..., g, & = К,... , 1;
з=i
- _ q _
3) вычисляется общая сумма скоростей V = Îy_1 =
3 = 1
Алгоритм состоит из следующих шагов. Первый шаг — задание начального распределения, второй — определение времени выхода по формуле (4).
3. Выбор элементарного события производится за К этапов. На каждом уровне начиная с первого выбирается случайное число £ Е [О, V^), г = 1,..., g, к = 1,..., К, и методом частичных сумм выбирается ветвь. Спускаясь все ниже, доходим до ветви К-го уровня, которая совпадает с ячейкой. Здесь выбираем и осуществляем элементарное событие.
4- Изменение состояния фрагмента. Далее происходит обновление матрицы состояния и значений скоростей событий и сумм скоростей в дереве событий.
Очевидно, что подобная реализация метода Монте-Карло более требовательна к оперативной памяти компьютера [8], так как необходимо хранить промежуточные значения сумм скоростей на каждой
ветви, но такие затраты оправданны, поскольку сложность представленного многоуровневого метода
х
на один шаг по времени составляет всего О (К х NVK ) — на выбор события и О (К х Ne х (Nn + 1)) — на обновление состояний и скоростей.
4. Сравнение результатов производительности. Каждый из представленных методов Монте-Карло был реализован последовательным алгоритмом на языке С++ и запускался на одном ядре процессора Intel Xeon Х5570. Каждым методом рассчитывалось 50 с реального процесса в имитационной модели реакции со скоростью диффузии 50 с-1. Исследовалась зависимость времени расчета от размера фрагмента решетки с N элементарными ячейками. Сетка изменялась от 20 х 20 до 100 х 100 с шагом 20 по обеим координатам. Тестовые расчеты производились без использования операций ввода/вывода, поэтому результирующее время отражает производительность непосредственно алгоритма. Исходный код программы можно найти в [12].
Зависимость времени расчета Т каждым из методов от размера поверхности N представлена на рис. 3, а, б. На рисунке используются следующие сокращения: D — динамический метод; К — кинетический метод; R — метод отказа; RF — метод частичных сумм; FK — многоуровневый метод, где К — число уровней дерева; FB — многоуровневый метод с бинарным деревом. Мы видим, что самым производительным по сравнению с другими методами оказался малоизвестный многоуровневый метод с К = 4 (рис. 3, б). Он превосходит самый медленный (метод отказа) примерно в 600 раз, а самый быстрый среди остальных (кинетический) — на два порядка. Один из наиболее популярных методов, динамический, позволил сэкономить часть вычислительных ресурсов (рис. 4), но проиграл по производительности кинетическому методу более чем в три раза. Метод частичных сумм не дал прироста производительности при решении данной задачи по сравнению с методом отказа. Среди многоуровневых методов алгоритм, основанный на бинарном дереве, оказался самым неэффективным, в то время как дерево всего с двумя уровнями показало неплохой результат с учетом того, что двухуровневый алгоритм достаточно легко реализовать.
Сравнение требуемых объемов оперативной памяти m каждым из представленных методов в зависимости от размера решетки N представлено на рис. 4. Мы видим, что здесь результаты не сильно расходятся, и разница между самыми затратными многоуровневыми методами и самым экономичным динамическим методом составляет примерно два раза. Это не является существенным недостатком многоуровневых методов, если учитывать современные возможности расширения объемов оперативной памяти на вычислительных станциях.
т, Ч
40 30 20 10
о
т, Ч 0,08 "
0,06 -0,04 -0,02 -
О 2 4 6 8 Л^х 1000
Рис. 3. Зависимость значений времени расчета Т от размера решетки Л/": а — каждым из методов Монте-Карло; б — для многоуровневых методов с различными К
т, КБ
18 -14 -
ю L^-1-1-1-1-
0 2 4 6 8 Nx 1000 Рис. 4. Результаты сравнения использования оперативной памяти т
5. Исследование динамики реакции. Многоуровневый метод позволил исследовать динамику рассматриваемой реакции на больших фрагментах решетки. В расчетах использовались решетки, содержащие до 4 млн узлов. Время расчетов существенно зависело от скорости поверхностной диффузии. Чем больше скорость диффузии, тем больше сумма скоростей всех событий и, следовательно, исходя из (4) и (6) меньше продолжительность пребывания системы в одном из состояний, а значит, для достижения заданного физического времени требуется больше итераций. На решетке 1000 х 1000 при скорости диффузии 200 с-1 потребовалось около 150 ч. расчетного времени для исследования процесса реакции в течение 200 с реального времени с сохранением состояний всех ячеек каждую 0.1 с.
Исследование динамики имитационной модели проводилось в области автоколебаний согласованной с ней точечной модели (системе ОДУ) [1]. Были изучены механизмы возникновения кинетических колебаний на микроуровне. Анализ полученных анимаций эволюции системы привел к выводу, что причиной возникновения колебаний являются латеральные взаимодействия, которые играют роль отрицательных обратных связей в слое адсорбата, необходимых для реализации колебаний. Полученные
в результате расчетов анимационные фильмы показали, что в ходе реакции в разных частях поверхности то зарождаются, то разрушаются "вихри". При специально заданных начальных условиях в модели построены плоские бегущие и спиральные волны (рис. 5).
Рис. 5. Движение спиральной волны на поверхности размером 1000 х 1000 с интервалом 4 с при значении скорости диффузии 200 с;-1
СПИСОК ЛИТЕРАТУРЫ
1. Куркина Е. С., М а к е е в А. Г. Бифуркационный анализ четырех-компонентной математической модели реакции (NO + C0)/Pt(100) // Обратные задачи естествознания. М.: Изд. отдел ф-та ВМиК МГУ, 1997. С. 52 78.
2. Куркина Е. С., А вер чу к Г. Ю. Моделирование методом Монте-Карло колебательной динамики гетерогенной каталитической реакции с латеральными взаимодействиями // Прикладная математика и информатика. № 41. М.: МАКС Пресс, 2012. С. 20 41.
3. Schulze Т. P. Efficient kinetic Monte Carlo simulation //J. Сотр. Phys. 2008. 227. P. 2455 24G2.
4. Battaile С. C. The kinetic Monte Carlo method: Foundation, implementation, and application // Comput. Moth. Appl. Modi. Engrg. 2008. 197. P. 338G 3398.
5. Gillespie D.T. Exact stochastic simulation of coupled chemical reactions // J. Phys. Chem. 1977. 81. N 25. P. 2340 23G1.
6. Whitesides R., Fronklach M. Detailed kinetic Monte Carlo simulations of grapheno-edgo growth // J. Phys. Chem. A. 2010. 114. N 2. P. G89 703.
7. Maksym P. A. Fast Monte-Carlo simulation of MBE growth // Semiconductor Sci. Tech. 1988. 3. P. 594 59G.
8. Blue J.L., Beichl I., Sullivan F. Faster Monte-Carlo simulations // Phys. Rev. E. 1995. 51. P. 8G7 8G8.
9. Ross S.M. Simulation. Burlington: Academic Press, 200G.
10. Gentle J.E. Random Number Generation and Monte Carlo Methods. Fairfax: Springer, 2005.
11. Куркина E. С. Автоколебания, структуры и волны в химических системах. Методы математического моделирования. М.: Изд. центр РХТУ им. Д. И. Менделеева, 2012.
12. https://github.com/newmon/monto-carlo_techs 1G. 12.2012.
Поступила в редакцию 19.12.12
COMPARATIVE ANALYSIS OF MONTE CARLO METHODS FOR SIMULATIONS OF COMPLEX DYNAMICS IN THE LATTICE MODEL OF THE CHEMICAL REACTION
Averchuk G. Yu., Kurkina E. S.
The complex oscillatory behaviour of the microscopic model of the chemical reaction on a catalyst surface is calculated by means of various Monte Carlo techniques. Five most popular algorithms are used to estimate the order of growth and real performance each of them. It is shown that the computation time can differ more than two orders of magnitude. The most efficient Monte Carlo algorithm was identified. Mechanisms of oscillations and waves were investigated on the microscopic level in the lattice model by means of this algorithm.
Keywords: microscopic stochastic model, heterogeneous catalytic reaction, Monte Carlo methods, self-sustained oscillations.
