CC BY
81
2012 Дискретные модели реальных процессов №3(17)
ДИСКРЕТНЫЕ МОДЕЛИ РЕАЛЬНЫХ ПРОЦЕССОВ
УДК 621.391.1:004.7
ИНВАРИАНТЫ КЛЕТОЧНО-АВТОМАТНЫХ МОДЕЛЕЙ РЕАКЦИОННО-ДИФФУЗИОННЫХ ПРОЦЕССОВ1
О. Л. Бандман
Институт вычислительной математики и математической геофизики СО РАН,
г. Новосибирск, Россия
E-mail: bandman@ssd.sscc.ru, bandman@academ.org
Вводится понятие инварианта клеточно-автоматной модели (КА-модели). Инвариант равен некоторой безразмерной характеристике моделируемого процесса, которая не зависит от его математического описания и может быть выражена как в модельных, так и в физических величинах. Знание инварианта важно для практического моделирования, поскольку оно позволяет вычислить масштабирующие коэффициенты для перехода от модельных величин к реальным и обратно. Приводятся примеры инвариантов нескольких известных КА-моделей реакционно-диффузионных процессов и предлагается основанный на инвариантах общий подход к решению задачи масштабирования КА-моделей.
Ключевые слова: клеточный автомат, клеточно-автоматное моделирование, нелинейная пространственная динамика, реакционно-диффузионные процессы, масштабирующие инварианты.
Введение
Современные компьютерные системы позволяют моделировать сложные физикохимические явления, имитируя превращения и движения реальных или абстрактных частиц в дискретных времени и пространстве. Широкий класс таких явлений составляют реакционно-диффузионные процессы. В них перемещения веществ представляются подчиняющимися законам диффузии, а их превращения (химические, фазовые, биологические) — нелинейными или пороговыми функциями, а также простыми подстановками вида «было — стало». Нелинейность и разрывность реакционной составляющей в такого рода явлениях делают затруднительным или совсем невозможным использование традиционных математических моделей, основанных на дифференциальных уравнениях. Поиски альтернативы привели к развитию дискретного моделирования, в частности к использованию идеи клеточного автомата (КА) фон Неймана [1].
Клеточно-автоматные модели реакционно-диффузионных процессов вызывают большой интерес [2, 3], поскольку позволяют представить сложный нелинейный процесс с помощью простых функций перехода одного конечного автомата. В теории динамических систем такие модели относятся к так называемым сложным системам (complex systems). Термин «сложные системы» введён в [4]. Его содержательный смысл состоит в том, что много простых вычислителей, функционируя совместно, проявляют себя как единый сложный процесс со свойствами самоорганизации, имитируя
хРабота поддержана грантами Президиума РАН (Проект №15.9, 2011), СО РАН (Интеграционный проект ИП-32б 2009), РФФИ (грант № 11-01-0567a).
разного рода автоволновые процессы, процессы образования устойчивых структур и другие диссипативные явления [5].
В настоящее время клеточно-автоматное моделирование реакционно-диффузионных процессов развивается интенсивно [2, 3] как в плане исследований [5], так и в плане применений. При этом для построения КА-модели требуется найти способ перехода от физико-химического описания моделируемого явления к КА-представлению и после получения результата моделирования необходимо его интерпретировать в привычных для исследователя физических или химических терминах. В ряде случаев это удаётся сделать исходя из опыта, иногда приходится ограничиваться качественным сходством с предполагаемым процессом. Хотя некоторые попытки решить эту задачу уже предприняты [6, 7], систематического метода построения КА-моделей по заданным параметрам процесса пока не существует. Чтобы заполнить этот пробел, в работе предлагается новый подход, основанный на понятии инварианта КА-модели.
Инвариантом далее называется безразмерная характеристика процесса, не зависящая от способа его математического представления, которая позволяет определить масштабирующие коэффициенты между реальными и модельными величинами. Поскольку КА моделируют пространственную динамику, то характеризующий КА-модель инвариант связывает величину пространственного (размер клетки в метрах) и временного (длительность итерации в секундах) масштабов. Отсюда следует, что в общем случае эти масштабы независимо выбираться не могут. Кроме того, КА-мо-дели реакционно-диффузионных процессов имеют по крайней мере две составляющие: диффузию и реакцию, каждая из которых характеризуется своим инвариантом. Инварианты таких сложных процессов не всегда возможно выразить в виде функции от инвариантов его составляющих, и масштабирующие функции не всегда возможно выразить аналитически через инвариант. В этих сложных случаях единственным путём определения инвариантов и масштабов остается вычислительный эксперимент и построение связывающих их таблиц и кривых.
1. Клеточный автомат
КА представляет собой множество одинаковых простых вычислителей, называемых клетками, которые задаются парами (и, х), где и Е А — состояние клетки из алфавита А; х Е X — имя клетки, равное вектору координат ^-мерного дискретного пространства X конечных размеров. В пространстве X определены подмножества, называемые шаблонами,
Тк(х) = {х, х + аь ..., х + Вп-1],
где а^, ] = 1,... ,п — 1 —векторы смещения координат х, п = |Т(х)|. Клетки с именами из Тк(х) образуют локальную конфигурацию
Бк (х) = {(ио, х), (иь х + а1),..., (ип-1, х + а^)}.
Множество клеток П = {(и,,х,) : х, € X} называется клеточным массивом, а перечень состояний клеток Па = (и1,и2,... ,и|Х|) — глобальным состоянием КА.
Функционирование КА задается локальным оператором и режимом его применения к клеткам из П. Локальный оператор 0(х) есть композиция более простых локальных операторов: 0(х) = Ф(01(х),... , 0^(х)). При этом 0к(х) может быть элементарной подстановкой в к (х). Подстановка выражается через локальные конфигурации следующим образом:
вк(х) : Бк(х) ^ Бк(х) (1)
где |Бк (х)| ^ |Бк (х)|, причём Тк С Тк и первые п' клеток в локальной конфигурации составляют базу подстановки, а остальные п — п' клеток играют роль контекста.
Подстановка в к (х) применима к клетке (и, х) Е П, если Б к (х) Е П, причём клетка (и, х) с переменным состоянием считается принадлежащей П, если область значений и равна алфавиту А. Применение подстановки в к (х) сводится к замене состояний базовых клеток (uj, х + а^) Е Б' (х) на значения
и'у = fj (и1,...,ип), 3 = 0,..., п', (2)
где fj (и1,... , ип) — функция перехода. Контекстные клетки не меняют своих состояний при применении подстановки (1).
В локальных операторах реакционно-диффузионных КА обычно используются следующие способы композиции подстановок: суперпозиция, случайный выбор одной из подстановок, применение подстановок в случайном порядке [8].
Применение 0(х) ко всем х Е X изменяет глобальное состояние Па (^ на новое Па(£ + 1). Такое изменение составляет итерацию. Итерация может выполняться разными способами, которые называются режимами функционирования КА. Основными из них являются синхронный и асинхронный.
При синхронном режиме на каждой £-й итерации выполняется следующее:
1) для всех (и, х) Е П(£) вычисляются новые состояния и'(х) путём применения к ним функции перехода (2);
2) во всех клетках (и, х) Е П(£) производится замена состояний и(х) на новые и'(х);
3) П(*) ^ П(* + 1).
Синхронные КА-модели реакционно-диффузионных процессов более изучены, чем асинхронные. Однако они отображают моделируемое явление либо косвенным образом, либо путём преобразования дифференциального уравнения в частных производных к клеточно-автоматному виду [9, 10]. Локальные операторы синхронных КА ограничены тем, что их подстановки (1) должны иметь одноклеточную базовую локальную конфигурацию, т. е.
Ув, Е 0(х) |Б'(х)| = 1, (3)
что следует из условий корректности вычислений [11]:
Ух, у Е X Ук,т е{1,..., |0(х)|} Тк (х) П Тт(у) = 0. (4)
Ограничение (3) затрудняет имитацию явления химических или биологических взаимодействий веществ по двум причинам. Во-первых, взаимодействия предполагают участие в них двух или более веществ, что требует изменений состояний двух или более клеток одновременно. Во-вторых, синхронное обновление состояний клеток не естественно, так как в моделируемом микромире нет синхронизации. Поэтому асинхронный режим оказывается часто более подходящим.
Асинхронный режим предполагает следующий порядок применений локального оператора:
1) с вероятностью р = 1/^| выбирается клетка (и, х) Е П;
2) к выбранной клетке применяется локальный оператор 0(х), и состояния клеток базовых локальных конфигураций (и', х) Е Б'(х) сразу же меняются на новые значения;
3) условно принимается, что ^| повторений пп. 1 и 2 составляет одну итерацию, такое соглашение удобно для сравнения синхронного и асинхронного режима и соответствует понятию одного шага метода, называемого кинетическим методом Монте-Карло [12, 13].
Поскольку в асинхронных КА локальный оператор применяется последовательно к выбранным клеткам, условие (4) всегда выполняется, т. е. вычисление эволюции КА всегда корректно. Это значит, что не может произойти потери данных (в одну клетку одновременно записаться два разных значения). Проблема корректности возникает только тогда, когда алгоритм асинхронного КА реализуется параллельно на нескольких процессорах [8, 14].
Режим функционирования КА является его существенным параметром, т. е. если два КА различаются только режимами функционирования, то они являются представлениями двух разных КА. Чтобы определить однозначно КА, надо задать четыре параметра: A, X, 0,р, где р Е {а, а}, а обозначает асинхронный режим, а — синхронный, а КА Ra и Ra обозначают асинхронный и синхронный КА соответственно.
Задание КА в общем случае не определяет характера его поведения. Эволюции одного и того же КА с разными исходными массивами П0 могут кардинально различаться. Ярким примером является всем известный КА «Игра ‘Жизнь’», который при разных начальных массивах порождает самые разные эволюции. Поэтому для определения КА-модели необходимы пять символов:
A = (A, X, 0,р, ОД).
Инвариант является характеристикой КА-модели. В общем случае один и тот же КА при разных начальных условиях имеет разные инварианты.
КА-модели реакционно-диффузионных процессов с явно заданными подстановками диффузии и реакции характеризуются инвариантами, которые зависят от свойств обеих составляющих, и, следовательно, имеет смысл рассмотреть инварианты каждой из них отдельно.
2. Инварианты КА-моделей диффузии
Процесс диффузии имеет два математических представления: дискретное (КА) и непрерывное (уравнение Лапласа), которые характеризуются безразмерным коэффициентом D = d ■ т/h2, где d — коэффициент диффузии, выраженный в м2 ■ с-1; т — шаг по времени в секундах; h — шаг по пространству в метрах. Поскольку коэффициент D не зависит от начальных условий, то КА диффузии можно отождествлять с КА-моделью, инвариант которой равен соответствующему D.
Из содержательного определения коэффициента диффузии как массы вещества, продиффундировавшего за единицу времени через единицу площади при градиенте концентрации равном единице (первый закон Фика для диффузии в изотропной среде), простым рассуждением можно получить значение D следующим образом. В булевом клеточном массиве на границе между областью, в которой все клетки имеют состояние и = 1, и областью, в которой состояния клеток и = 0, градиент концентрации равен 1. Следовательно, масса вещества (количество частиц), продиффундировавшая за At = 1, равна вероятности того, что смежные пограничные клетки обменяются состояниями. Однако, поскольку такое рассуждение нельзя считать доказательством, ниже приводится его подтверждение результатами вычислительных экспериментов.
Наиболее известны два КА с булевым алфавитом, моделирующих процесс диффузии: 1) синхронный КА, названный в [15] КА с окрестностью Марголуса, и 2) асинхронный КА, получивший название наивной диффузии [16].
Синхронный КА диффузии = (A, X, 0), A = {0,1}, X = {x : x = (i, j)}, функционирует следующим образом. На множестве X выделяются два подмножества: X0 = = {(i,j) : i mod 2=0, j mod 2=0} — чётное и X1 = {(i,j) : i mod 2=1, j mod 2=1} —
нечётное. Вводится шаблон T2x2 = {(i, j), (i, j + 1), (i + 1, j + 1), (i + 1, j)}, который определяет локальную конфигурацию:
S ^ j) = {(u0, (i, j)), (u1 (i, j + 1)), (u2, (i + 1, j + 1)), (u3, (i + 1 j))}.
Вероятностный локальный оператор 0(i, j) перемещает состояния клеток таким образом, что новые состояния и'к, k = 0,1, 2, 3, равны
U __ / u(k+1) тоё4> если rand < p
k \ u(k-1) mod 4, если rand ^ (1 - p),
где p — вероятность применения 0(i, j). Итерация подразделена на две стадии: на первой стадии 0(i, j) применяется синхронно к клеткам (i, j) Е X0, на второй — к клеткам (i, j) Е X1. В [15] сделана попытка доказать, что при p = 0,5 безразмерный коэффициент диффузии D2o- = 1,5. Эта величина поставлена под сомнение в [17], где утверждается, что в [15] не учтена разница в вероятности перемещения частицы в противоположных направлениях на чётных и нечётных стадиях, и показано, что при p = 0,5 D2o- = 1, что согласуется с приведённым выше рассуждением. При p < 0,5 D2o- = p.
Одномерный вариант синхронной диффузии функционирует в пространстве X = {i = 0,1, 2,... , N}, подразделённом на три подмножества: X = X0 U X1 U X2, где Xk = {i : i mod 3 = k}, k = 0,1, 2. Локальный оператор 0(i) производит обмен состояниями клетки (u0,i) со случайно выбранным соседом (u1,i — 1) или (u2,i + 1):
J (u0 = u1) & « = u0), если rand < p,
[ (u0 = u2) & (u'2 = u0), если rand ^ (1 — p).
Итерация подразделена на три стадии: на первой стадии 0(i) применяется синхронно
к клеткам i Е X0, на второй — к i Е X1, на третьей — к i Е X2, коэффициент диффузии
D1a = p.
Асинхронный КА диффузии Na = (A, X, 0), A = {0,1}, X = {x : x = (i, j)} (наивная диффузия) предполагает, что выбранная согласно асинхронному режиму равновероятно клетка из X обменивается состоянием с одним из своих соседей, тоже выбранным с равной вероятностью.
В двумерном случае используется 5-точечный шаблон с локальной конфигурацией
S(x) = {(uo,x), (u1,x + a1)), (u2,x + a2), (u3,x + аз), (u4,x + a4)}.
Локальный оператор 0(x) применяется к клеткам (u0, x) и выполняет обмен состояниями с одной из соседних клеток (uk, x + ak) Е S(x), выбранных равновероятно:
(u0 = uk) & (uk = u0), если (k — 1)/4 < rand < k/4) & rand < p, k = 1, 2, 3, 4. (5)
Из (5) следует, что средняя вероятность того, что за одну итерацию частица переместится из клетки (1, x) в одну из соседних клеток (0, x + ak), равна 0,5, т. е. при p = 1 D2a = 0,5, при p < 1 D2a = 0,5 ■ p.
В одномерном варианте X = {i : i = 0,1,..., N}. Случайно выбранная клетка (u0, i) производит обмен состояниями со случайно выбранным соседом (u1, i — 1) или (u2, i + 1) с вероятностью p. Функция переходов и, следовательно, коэффициент диффузии такие же, как в синхронном случае, D1o- = p.
Целочисленные КА-модели диффузии предложены в [18], где они названы «многочастичными». В [18] даны две версии целочисленной КА-модели диффузии, соответствующие синхронной и асинхронной булевым моделям. Состояние клетки выражается в виде суммы двух частей: V = 7 • V + (1 — 7)1», 7 Е (0,1], из которых только одна совершает перемещение в соседнюю клетку. Такой принцип позволяет уменьшить или даже совсем избавиться от автоматного шума и, следовательно, от необходимости осреднения булевого массива.
Синхронная целочисленная 2Б КА-диффузия использует алфавит А = {0,1, 2,... , М}, дискретное пространство X = {х : х = (г,3)}, подразделённое на два подмножества клеток: чётное Xo и нечётное X!, как и в булевом случае. Оператор перехода 0(г,3) выполняет следующую функцию переходов:
vk/
(1 — Y)vk + Y ■ V(k+1) mod 4, если rand < p,
(1 — Y)vk + Y ■ V(k— 1) mod4, если rand ^ (1 — p),
где p — вероятность применения 0(i,j); Y Е (0,1] —перемещаемая часть величины v, k = 0,1, 2, 3. Как и в булевом варианте, итерация подразделена на две синхронные стадии: на первой стадии 0(x) применяется к клеткам x Е X0, на второй — к клеткам x Е X1 .
Одномерный вариант синхронной целочисленной КА диффузии функционирует в пространстве X = {i : i = 0,1, 2,... , N}, подразделённом на три подмножества: X = X0 U X1 U X2, где Xk = {i : i mod 3 = k}, k = 0,1,2. Локальный оператор
0(i) производит обмен частью состояния y ■ v0 со случайно выбранным соседом из
S(i) = {(vo, i), (V1, i — 1), (V2, i + 1)}:
(v0 = (1 — y)v0 + Y ■ v1) & (v^ = (1 — y)v1 + Y ■ v0), если rand < p,
(v0 = (1 — Y)vo + Y ■ V2) & (v2 = (1 — Y)v2 + Y ■ vo), если rand ^ (1 — p).
Как и в синхронном одномерном КА, итерация подразделена на три стадии: на k-й стадии 0(i) применяется к Xk.
Асинхронная целочисленная диффузия работает аналогично наивной диффузии. Отличие как в двумерном, так и в одномерном случаях состоит в том, что межклеточный обмен происходит частью от значения состояния. В двумерном случае функция
перехода переводит состояния (v0,x) и (vk,x + ak) в (v0,x) и (vk,x + ak) следующим образом:
v0 = <; — Y>vo + Y ■ vk' (k — 1)/4 < rand < k/4.
vk = (1 — Y)vk + Y ■ vo,
В одномерном случае применяется та же функция перехода (2), что и в синхронном случае, но без разделения итерации на стадии.
Коэффициенты диффузии целочисленных КА равны соответствующим коэффициентам булевых КА, умноженных на Y.
Знание инварианта позволяет определить масштабы параметров КА-модели, т. е. длину h стороны клетки в метрах и время итерации т в секундах. При этом один из этих масштабов может быть выбран из условий задачи, а второй вычисляется из соотношения
D = ж- <6>
где d — коэффициент диффузии среды в м2■с-1. Например, h может быть выбран равным размеру реальных или абстрактных частиц. Тогда т = (D ■ h2)/d с.
Для подтверждения правильности предположений о значениях коэффициентов диффузии проведены следующие вычислительные эксперименты. Для каждого случая выполнялось моделирование процесса диффузии двумя методами: 1) клеточным автоматом при предполагаемом значении Dca и 2) численным решением уравнения Лапласа при Dl = Dca. Одинаковыми были также размеры I х J клеточного и сеточного пространств, а также начальные условия: (v(0, x)) = u(0, у) для всех x Е X и всех у Е Y, где Y — множество узлов в численном решении уравнения Лапласа. Производилось сравнение значений (v(t, (I/2, j))) с u(t, (I/2, j) для j = 0,... , J и t = 0,..., T, где T = 300 — число итераций, в течение которых процесс практически завершался. Если значения (v(t, (I/2, j))) и u(t, (I/2, j) совпадали (с точностью до ошибки осреднения состояний КА), то считалось, что Dca = Dl — инвариант соответствующей КА-модели.
Результаты вычислительного моделирования подтвердили предполагаемые значения инвариантов для всех КА-моделей диффузии, в которых D = 1 для одномерной и двумерной синхронной КА-модели, а также для одномерной асинхронной; D = 0,5 для двумерной наивной КА-модели. Все коэффициенты диффузии целочисленных КА равны соответствующим коэффициентам булевых КА, умноженных на Y.
3. Инварианты КА-моделей реакции
В КА-моделях химических реакций за инвариант следует принять безразмерный коэффициент R = k ■ т, где k — константа скорости реакции с размерностью с-1, которая для большинства изученных реакций известна; т — длительность итерации клеточного автомата в секундах. Следует отметить, что инвариант КА-модели химической реакции не зависит от размеров клетки и, следовательно, масштабирующие коэффициенты h и т не зависят друг от друга, и h может быть выбран исходя из условий задачи. Физический смысл инварианта R — относительное изменение концентрации (AC/C) реактанта за время одной итерации т. Величина т и является тем значением, которое связывает модельные величины с реальными. Если в модели участвует только одна реакция, то т естественно выбирать равным длительности этой реакции. Если в моделируемом процессе участвует несколько реакций и, возможно, другие действия (адсорбция, диффузия), то соотношения скоростей регулируются вероятностями применения соответствующих локальных операторов, зависящих от всех инвариантов. Так, например, если в процессе участвует n действий, описываемых n подстановками, т. е. 0(x) = {01(x),... , 0n(x)}, то вероятность выполнения каждой #j(x) равна
Rk
pi = n = n . (7)
Е Rj Е kj
j=1 j=1
В этом случае инвариантом процесса следует считать набор R = {R1,... , Rn}, в котором R^ = k^ ■ т, а время итерации т выбирается таким, чтобы самое медленное действие закончилось.
Например, реакция, которая получила название ZGB-реакции по фамилиям авторов [19] и положила начало широкому применению асинхронных вероятностных КА (кинетического метода Монте-Карло), представлена тремя химическими формулами:
0 —COads,
0 + 0 —^ Oads + Oads,
COads + O
ads ^
CO2(gas).
На поверхность катализатора из газа адсорбируются моноокись углерода и диссоциированный кислород, которые, реагируя, образуют углекислый газ. Моноокись углерода диффундирует на пустые клетки, углекислый газ освобождает поверхность катализатора. Значения констант адсорбции (А^) и десорбции (к2) зависят от начальных условий процесса, они пропорциональны парциальным давлениям рН(СО), рН(02) соответствующих газов, значение константы реакции (А3) от начальных условий не зависит. Реакция окисления происходит мгновенно, как только молекулы СО и 02 окажутся в смежных клетках.
КА-модель этого процесса = (А, X, 0, П(0)) определяется значениями А = (0,
О, СО}, X = (х = (г,з)}, П(0) = ((0, х) : х € X}, 0(х) = Ф5(ФГ(^(х),^, (х)),0з(х)), где Ф5, Фг —локальные операторы случайного выбора и суперпозиции подстановок соответственно:
01 (х) #2 (х)
#3(х)
(0,х) —-4 (СО,х),
(0 х^ (0, х + а} ^ {(° (і,з)),(0, (х + з.г)}, (СО, х), (О, (х + аг)} —4 {(0, х), (0, (х + аг)},
Вероятности равны
кі
Рі
рН(СО)
&2
кі + ^2 рН(СО) + рН(02):
Р2
I = 1, 2, 3,4.
рН(О2)
кі + &2 рН(СО) + рН(О2)‘
Инвариант этой модели равен Я = (Я^Я2), где Я1 = А4 • т, Я2 = А2 • т. Масштаб
времени т следует выбрать равным времени адсорбции, так как адсорбция наиболее медленное действие. Масштаб длины к не зависит от инварианта, его можно принять равным размеру молекулы СО. Результаты моделирования легко перевести в физические величины: например, количество образующегося углекислого газа за секунду равно
М(С02)= *<С0*>т(С0*> кг. с-1,
т
где N(С02) —количество применений 03 за одну итерацию; т(С02) —масса молекулы СО2.
4. Инвариант КА-модели «распространение фронта»
К классу реакционно-диффузионных относятся также явления, в которых реакция выражена нелинейной функцией (чаще всего полиномом). В традиционной математике эти явления описываются дифференциальными уравнениями вида
и* = + Я (и), (8)
где правая часть содержит оператор диффузии (Лапласиан) и нелинейную функцию Я (и), отображающую реакционную составляющую процесса.
Типичным примером является уравнение «диффузии, соединенной с возрастанием количества вещества» [20, 21] с нелинейной функцией вида
Я (и) = аи(1 — и)
и начальными условиями
Я(0) = Я(1) = 0, 0 < а ^ 1, 0 ^ и ^ 1.
В химии это уравнение моделирует реакцию превращения вещества B в вещество A с вероятностью, пропорциональной плотности А, при условии, что A+B=const, т. е.
A+B —^ A+A.
В экологии уравнение (8) моделирует распространение эпидемий, сорняков, животных, в социологии — распространение слухов, привычек и т. д. Во всех случаях уравнение вида (8) моделирует пространственную экспансию, и процесс называется распространением фронта. Инвариантом процесса является скорость распространения фронта.
Так как в [20] даны аналитические оценки для одномерного случая, то далее исследуется одномерный аналог уравнения (8) A = (A,X, 0, П(0)), который построен как композиция КА диффузии и КА реакции [22] при A = {0,1}, X = {i :
i = —L,... , 0,..., N}, П(0) = {(1, i) : i < 0; (0, i) : i ^ 0}:
0(i) = $sum(01(i), 02(i)). (9)
Композиционный локальный оператор Ф^т выполняет покомпонентное арифметическое суммирование глобальных состояний Па(01) и Па(02) с преобразованием булева алфавита в вещественный и обратно [23], 01(i) —локальный оператор КА-модели диффузии, 02(i) : (u,i) ^ (u', 1) —одноклеточный бесконтекстный локальный оператор с функцией перехода u' = F(u).
Таким образом, состояния w(i) Е Па^+1), полученные в результате применения (9) к (v,i) Е n(t), равны
w(i) = Disc((v'(i)) + F ((v'(i)))), (10)
где v'(i) — значение функции перехода оператора 01(i) к (v, i) Е n(t); Disc(x) — операция, преобразующая состояния x Е {0,1} в у = (0,1) следующим образом [23]:
_ , , _ J 1, если rand < у,
X DiSC(у) | л
0 в иных случаях.
В [20] аналитически доказано, что скорость распространения фронта при t ^ то зависит от коэффициентов диффузии (D) и скорости реакции (а) следующим образом:
V0 = 2 VD ■ а.
Много позднее в [24, 25] было обнаружено, что дискретный характер компьютерного моделирования влияет на скорость распространения фронта в сторону её уменьшения. Отсюда следует, что полученные аналитическим путём асимптотические оценки не всегда соответствуют результатам клеточно-автоматного моделирования, т. е. значение инварианта необходимо устанавливать путём вычислительных экспериментов.
Моделирование было проведено для КА-модели и численного решения уравнения (8). Численное решение u(t,x) сравнивалось с осреднёнными состояниями (w(t, j)) результатов КА-моделирования (10). Для более точного осреднения в одномерном случае клеточный массив КА целесообразно взять прямоугольным, т. е. X = {(i, j) : i = = 0,... , I; j = 0,..., J}, оставив локальный оператор (9) моделирующим распространение фронта вдоль одной оси j. Такой приём перехода от одномерного случая к ква-зидвумерному позволяет выбирать достаточно большую окрестность осреднения, захватывающую всю ширину фронта Av(j) = {(0, j),... , (I, j)}.
В процессе моделирования вычислялась скорость распространения фронта в зависимости от значений коэффициента диффузии D от 0,2 до 3 при а = 0,5 (рис. 1).
Для преодоления ограничения Куранта при решении уравнения (8) с D > 0,2 каждая итерация составлялась из m = D/0,2 временных шагов, в которых прибавление реакционной части выполнялось с вероятностью p = 1/m. Размер клеточного массива для всех испытаний равен I х J = 200 х 800. Начальные значения клеточного массива: v(i,j) = 1 при 0 < i < 200, j < 40; в остальных случаях v(i,j) = 0. Скорость распространения фронта V вычислялась в процессе моделирования на участке межДУ j 1 = 500 и j2 = 700 следующим образом: V = (j2 — j1)/(t2 — t1) при t2,t1, таких, что u(j2,t2) = u(j1,t1) = 0,5. Эксперименты показали, что скорость фронта перестаёт изменяться уже после t = 200 при 0,1 < а ^ 1. Поэтому измеренные при j > 500 значения можно считать установившимися.
3
0 т--------1-------1-------1-------1--------1-------1----
0,2 0,5 1 1,5 2 2,5 3
Коэфф. диффузии
—♦-КА -0-Диф.ур-ие -А-Теор.оценка
Рис. 1. Полученные путём моделирования зависимости скорости распространения фронта от значений коэффициента диффузии при а = 0,5
5. Инвариант КА-модели диффузионно-ограниченной агрегации
Эта КА-модель (diffusion limited aggregation) [26] отображает целый ряд явлений, связанных с фазовыми переходами. Например, кристаллизация в растворах, образование льда, снежных хлопьев, рост кораллов и др. Наиболее известны асинхронные КА-модели, в которых диффузионная составляющая описывает случайное блуждание частиц в пространстве, а реакционная составляющая преобразует частицу в неподвижную, если она касается другой неподвижной. Процессу соответствует КА-модель Aa = {A, X, 0, П(0)), где A = {a, b, c}, X = {(i, j)}, П(0) —исходный клеточный массив, в котором > 90 % случайно распределённых клеток имеют состояния b, < 10 % — состояния a и несколько клеток-зародышей находятся в состоянии c. Локальный оператор 0(i,j) = (01(i,j),02)), где 01(i,j) —оператор наивной диффузии (5),
#2(i,j) : {c, (i,j)), (a, (i + k,j + I))} {c, (i,j)), (c, (i + k,j + I))}, (11)
где (k,l) E {(0,1), (1,0), (—1, 0), (0, —1))}; p — коэффициент прилипания. В процессе эволюции вокруг зародышей образуется растущая структура (рис. 2).
^ = 2000 ^ = 30000 f = 50000
Рис. 2. Эволюция КА, моделирующего процесс диффузионноограниченной агрегации. Размер клеточного массива 400 х 400. Коэффициент прилипания р = 1, 5 = 1,64
Коэффициент прилипания р определяет однозначно эволюцию КА при заданном П(0), а также её главную характеристику — фрактальную размерность:
6 = 1о^ (Я))/ log('кЯ2),
где Я — радиус круга с центром в зародышевой клетке; N (Я) —число клеток в состоянии с в этом круге.
Инвариантом КА-модели следует считать коэффициент прилипания. Он однозначно определяет главную характеристику реального явления диффузионно-ограниченной агрегации. Аналитическое соответствие между р и 6 пока не установлено. Известны только полученные путем моделирования таблицы [26]. Таким образом, чтобы построить КА-модель Аа процесса агрегации вещества а из раствора Ь с заданной фрактальной размерностью 6, необходимо выполнить следующее:
1) выбрать размер клетки к равным размеру частицы (молекулы, гранулы) агрегируемого вещества и определить исходный клеточный массив;
2) по (6) определить шаг по времени т, исходя из инварианта наивной диффузии Б = 0,5 и физического коэффициента диффузии вещества й;
3) по заданной величине 6 определить по таблице из [26] соответствующее значения коэффициента прилипания р;
4) записать (5) и в2 (11) и задать П(0).
Заключение
Введено понятие инварианта КА-модели, при помощи которого предлагается решить проблему определения масштабирующих коэффициентов, позволяющих построить КА-модель по заданному физическому описанию процесса, а также выразить результаты моделирования в привычных физических понятиях (метрах, секундах, килограммах).
Рассмотрены инварианты для КА-моделей ограниченного класса явлений, называемых реакционно-диффузионными процессами, хотя предлагаемый подход может быть распространён и на КА-модели любого типа. Приведённые примеры показывают, что для КА-моделей, в которых диффузия и реакция заданы явно, инварианты могут быть определены исходя из известных инвариантов диффузии и реакций. В более сложных КА-моделях инварианты должны быть получены путём моделирования и сохраняться в виде таблиц и кривых, подобно тому, как это делалось при определении физических и химических констант в традиционном моделировании.
ЛИТЕРАТУРА
1. Фон Нейман Дж. Теория самовоспроизводящихся автоматов. М.: Мир, 1971. 240 с.
2. Ванаг В. К. Диссипативные структуры в реакционно-диффузионных системах. Эксперимент и теория. Ижевск: ИКИ, 2008. 300 c.
3. Boccara N. Reaction-Diffusion complex systems. Berlin: Springer, 2004. 397 p.
4. Wolfram S. Statistical mechanics in cellular automata // Rev. Mod. Phys. 1983. V. 55. P. 601-642.
5. Ванаг В. К. Исследование пространственно-распределенных динамических систем методами вероятностного клеточного автомата // УФН. 1999. Т. 169. №5. C. 481-504.
6. Bandman O. Mapping physical phenomena onto CA-models // AUT0MATA-2008. Theory and Applications of Cellular Automata. London: Luniver Press, 2008. P. 381-397.
7. Бандман O. Л. Отображение физических процессов на их клеточно-автоматные модели // Вестник Томского госуниверситета. Сер. Математика. Кибернетика. Информатика. 2008. №3. C. 345-356.
8. Kalgin K. Comparative study of Parallel Algorithms for Asynchronous Cellular Automata Simulation on Different computer Architectures // LNSC. 2010. V. 6350. P. 399-408.
9. Weimar J. Cellular Automata for reaction-diffusion systems // Parallel Computing. 1997. V. 23. No. 11. P. 1699-1715.
10. Bandman O. A Hybrid Approach to Reaction-Diffusion Simulation // LNCS. 2001. V. 2127. P. 1-16.
11. Ачасова С. М., Бандман О. Л. Корректность параллельных процессов. Новосибирск: Наука, 1998.
12. Jansen A. P. J. An Introduction to Monte Carlo Simulations of Surface Reactions. arXiv: cond-mat/0303028v1[stst-mech]. 2003. 51 p.
13. Metropolis N. and Ulam S. The Monte Carlo Method // Amer. Stati. Assoc. 1949. V. 44. No. 247. P. 335-341.
14. Bandman O. Parallel Simulation of Asynchronous Cellular Automata Evolution // LNCS. 2006. V. 4173. P. 41-47.
15. Малинецкий Г. Г., Степанцов М. Е. Моделирование диффузионных процессов с помощью клеточных автоматов с окрестностью Марголуса // Журн. вычисл. математики и мат. физики. 1998. Т.38. №6. С. 1017-1031.
16. Toffolli T. and Margolus N. Cellular Automata Machines: A New Environment for Modeling. USA: MIT Press, 1987. 260 p.
17. Калгин К. В. Клеточно-автоматные модели процесса диффузии // Частное сообщение.
18. Medvedev Yu. G. Multiparticle Cellular Automata for Diffusion Simulation // LNCS. 2011. V.6083. P. 204-211.
19. ZiffR.M., GulariE., and Barshad Y. Kinetic phase transition in an irreversible surface-reaction model // Phys. Rev. Lett. 1986. V. 56(24). P. 2553-2556.
20. Колмогоров А. Н., Петровский И. Г., Пискунов И. С. Исследование уравнения диффузии, соединенной с возрастанием количества вещества и его применение к одной биологической проблеме // Бюллетень МГУ, секция А. 1937. Вып. 6. С. 1-25.
21. Fisher R. A. The genetical theory of natural selection. Oxford: Univ. Press, 1930.
22. Бандман O. Л. Методы композиции клеточных автоматов для моделирования пространственной динамики // Вестник Томского госуниверситета. 2004. №9(1). С. 183-193.
23. Bandman O. Synchronous versus Asynchronous Cellular Automata for Simulating NanoKinetics // Bul. Nov. Comp. Center. Iss.21. Novosibirsk: NCC Publisher, 2004. P. 9-21.
24. Warren C., SomfaiE., and Sander L. M. Velocity of front propagation in 1-dimensional autocatalytic reactions // Braz. J. Phys. 2000. V. 30. No. 1. P. 157-162.
25. Brunet E. and Derrida B. Effect of Microscopic Noise on Front Propagation Velocity // J. Stat. Phy. 2001. V. 103. No. 1-2. P. 269-282.
26. Ackland G. J. and Tweedie E. S. Microscopic model of diffusion limited aggregation and electrodeposition in the presence of leveling molecules // Phys. Rev. E. 2006. V. 73. Iss. 1. 011606.
