CC BY
73
Паньгин Александр Викторович
МОДЕЛИРОВАНИЕ В MS EXCEL ТЕЧЕНИЯ ЖИДКОСТИ С ПОМОЩЬЮ КЛЕТОЧНОГО АВТОМАТА
ВВЕДЕНИЕ
Предлагается разработка в электронной таблице MS Excel приложения для моделирования процесса течения жидкости в канале с наличием препятствий различной формы. Использование простейшей техники клеточного автомата (КА) позволяет наглядно представить процессы, для описания которых требуется громоздкий аппарат вычислительных методов решения сложной системы дифференциальных уравнений.
Двумерный КА можно представить в виде прямоугольной сетки m X n, в которой за квантованный промежуток времени At числовые значения в каждой клетке изменяются по некоторому заданному закону, в зависимости от значений (как правило) соседних клеток на предыдущем временном шаге.
Наиболее известным вариантом КА является игра «Жизнь» Дж. Конвэя. Клеточные автоматы являются прообразами современного нейросетевого моделирования.
При изучении «протекания» сложных процессов в исследуемой системе (объекте) их не всегда удается описать аналитическими зависимостями, однако применение стохастических (вероятностных) методов на уровне локальных взаимодействий способно правильно отразить поведение системы в целом [1].
Табличное представление данных в MS Excel, встроенный язык VBA являются
удобными и понятными средствами для моделирования и визуализации работы КА.
В качестве самостоятельной работы с помощью КА предлагается смоделировать и отыскать решение одной олимпиадной задачи.
ПОСТРОЕНИЕ МОДЕЛИ
Определим канал, по которому «течет жидкость», как поле 11X19, состоящее из набора A2:K20 ячеек Excel. Отформатируем его надлежащим образом (рисунок 1), изменив ширину столбцов и прорисовав границы ячеек рабочего поля. Жидкость представляем как «количество молекул» в заданной ячейке. Первоначально поле пустое. Зададим препятствия в поле (удобно значениями -100, при малом размере ячеек «препятствие» отобразится, как ##). Предполагаем, что за временной шаг каждая молекула двигается вниз на следующий ряд по некоторым предопределенным правилам.
ПРАВИЛА КЛЕТОЧНОГО АВТОМАТА
Обозначим, используя обозначения [2], через uij(tk) - число молекул в клетке поля (i - строка, j - столбец) в момент времени tk. Предполагаем, что время изменяется дискретно tk+1 = tk + At, (k = 0, 1, ...), в один и тот же момент совершается движение молекул для каждого ряда сетки и для каждой заня-
тот же момент времени). Ограничение d по пункту 2 снимается при смещении вдоль ряда, в котором находится молекула. На одной из сторон сетки (во «входящем» ряду в первой строке Excel) предполагается присутствие за единицу времени потока новых молекул, движение которых регулируется теми же самыми законами.
ЗАДАНИЕ ИСХОДНЫХ ДАННЫХ
Число молекул каждой ячейки входящего ряда A1:K1 задается пользователем, что удобно для «экспериментирования». В ячейке M8 хранится число d, а в ячейке M5 - число итераций (tstep) по времени. Цветовая палитра определяет раскраску ячейки поля канала в зависимости от содержимого значения (числа молекул). Диапазон ячеек M12:M19 определяет цвет, а ячеек N12:N19 - верхние границы числовых значений для соответствующего цвета (изменяются автоматически при задании числа d). Цвет и числовые границы палитры могут быть переопределены пользователем для достижения определенных эффектов исследования. Ячейки препятствия пользователь образует (или стирает) двойным кликом мыши (или функциональными возможностями Excel).
При нажатии на кнопку «Сделать итерации» можно наблюдать изменение состояния рабочего поля за tstep временных шагов, при этом в ячейке N5 накапливается суммарное число проведенных итераций. Порядное видоизменение ячеек дает зрительный эффект динамики течения жидкости. В ячейке N8 отражается максимальное число молекул, достигнутое в некоторых ячейках поля в ходе «эксперимента». В течение эксперимента пользователь может произвольно изменять входные данные.
При нажатии на кнопку «Очистить поле» эксперимент приводится в исходное начальное состояние.
ЭКСПЕРИМЕНТЫ
1) Модели движения жидкости, описанные с помощью КА по изложенным прави-
Рисунок 1.
той позиции. Правила, по которым формируется модель используемого клеточного автомата, схематично описываются на рисунке 2.
1. Если uij(tk) = p, то каждая молекула в позиции (i, j) (на рисунке 2 ячейка светлосерого - на экране голубого - цвета) может занять новую позицию (i + 1, j + r) в момент времени tk+1, где r - целое случайное число из значений -1, 0 и 1.
2. Если Mi+1 j+r (tk+x) > d, где d - заданное максимальное число молекул, которое может поместиться в одну отдельную ячейку (более темная тональность серого - на экране синего - на рисунке 2) при продвижении в следующий ряд, то молекула остается в i - ряду, и она может двигаться в позицию (i, j+r);
3. Если uj+l j+r(tk+1) < 0, то есть в ячейке (/ + 1, j + г) препятствие, то молекула не может двигаться в эту позицию, она остается в i-ряду и может двигаться в позицию (i, j + r).
Смещение молекул i-ряда (то есть расположенных в одной строке Excel) вычисляется после обработки молекул i + 1 ряда (для каждого ряда поля, начиная с нижнего, в один и
Рисунок 2.
лам, являются слишком простыми. В этом их недостаток, так как не учитываются межмолекулярные взаимодействия, эффекты турбулентности (возможность образования вихрей, движения в обратном направлении). Но даже такие простые модели дают качественную картину происходящих процессов. Например, если увеличивать плотность потока жидкости в ячейке входного ряда, то, как видно из рисунка 1, конфигурация препятствий не может обеспечить направление потока в нужное «русло», образуется «перелив» в нежелательный вторичный «рукав».
На рисунке 1 также виден эффект накапливания молекул на границе препятствий. Плотность заполнения ячеек визуально оценивается градацией темных тонов синего цвета, что позволяет оценить изгиб линий потока.
2) Можно использовать модель движения жидкости для других исследований. Создадим из препятствий некое «сооружение» для изучения свойств треугольника Паскаля [1] (рисунок 3). Будем бросать из некоторой точки (ячейка Р1) по 100 шаров-молекул за один шаг по времени (число итераций можно задать 10). По мере разброса большого числа шаров (к примеру, 5000) значение в ячейке Р1 обнуляется, и итерации продолжаются, пока все шары не распределятся по «лузам» (итоговые количества указаны в ячейках Р2:Р7). В столбце О приведены нормированные значения столбца Р по количеству шаров в одной
плмиЛра onfeefeM&effi- раасраасу $лейки каЛамл б ул&исимлс&и cofeft^uwotoo fftaieftufr (гисма мдмещм).
из крайних луз. После выбора числового целочисленного формата ячеек столбца, можно отметить, что значения совпадают с биномиальными коэффициентами Ск (для п = 5, к = 0, ..., 5).
Упражнение. Можно ли построить преграду для распределения потока жидкости из одного «сопла» по трем руслам в количественной пропорции
a) 1:1:1
b) 1:2:3
ОПИСАНИЕ ПОСТРОЕНИЯ МАКРОСОВ
1. Определим переменные т и }т как число строк и столбцов поля (канала), мас-
Mo^fto ми nocffykouffai преграду распределения no&ootcA фсирыос&и Uf offtotoo «мама» no Лрем русмам...
Рисунок 3.
сив carry для регистрации молекул текущего ряда, перенос которых должен произойти в горизонтальном направлении (вдоль ряда их месторасположения в текущий момент времени), массив carry1 для регистрации молекул, переносимых в следующий ряд. Переменная d определяет ограничение числа молекул в ячейке для предотвращения переноса в следующий ряд.
2. Откроем в меню Вид ® Панель инструментов ® Формы, создадим кнопку с названием «Сделать итерации» и назначим ей имя макроса Iteration.
3. Создадим в окне редактора VBA подпрограмму, в которой по числу итераций по времени (переменная цикла tstep) для каждого ряда (переменная цикла i), начиная с максимального, выполняется
- обнуление массива carry, carry 1;
- для каждой ячейки данного ряда (переменная цикла j) перенос (процедура uv) числа частиц в ячейки последующего или текущего ряда;
- раскраска ячеек поля (процедура Paint) (см. листинг 1).
Листинг 1
Dim im, jm, d, carry, carryl Sub Iteration() im = 20 jm = 11
d = Range("M8") For tstep = 1 To Range("M5") Application.ScreenUpdating = False For i = im To 1 Step -1 ReDim carry(jm) ReDim carry1(jm) For j = 1 To jm
Call uv(i, j) Next j
For j = 1 To jm
Cells(i, j) = Cells(i, j) + carry Cells(i + 1, j) = Cells (i + 1, j) Next j Next i
Range("N5") = Range("N5") + 1 Call Paint(im, jm) Application.ScreenUpdating = True Next tstep End SubM
4. Процедура uv(i,j) осуществляет перенос молекул из позиции (i, j) согласно описанным правилам. Процедура frnd генерирует случайное число из значений {-1, 0, 1} (см. листинг 2).
ВСПОМОГАТЕЛЬНЫЕ МАКРОСЫ
5. Макрос Paint раскрашивает ячейки рабочего поля согласно заданной палитре цветов и определяет максимальное значение «плотности» (количество молекул в ячейке) в течение эксперимента (см. листинг 3).
Образуем, макрос Сеаг, коийорий агшцл&й ейки рабогего
(j)
+ carryl(j)
6. Образуем макрос Clear, который очищает ячейки рабочего поля, не занятые препятствиями, для подготовки к проведению последующего эксперимента. Откроем в рабочем листе Excel меню Вид ® Панель инструментов ® Формы, создадим кнопку с названием «Очистить поле» и назначим ей имя макроса Clear (см листинг 4).
7. В редакторе VBA в окне VBAProject щелкнем дважды кнопкой мыши на Excel объекте «Лист1» и в открывшемся окне создадим макрос, который стандартно вызывается в случае двойного клика кнопкой мыши на ячейке таблицы. В теле макроса зададим операторы, которые в случае, если целевая ячейка (Target) принадлежит ре-
Листинг 2
Sub uv(i, j)
If Cells(i, j) = -100 Then Exit Sub 'препятствие
km = Cells(i, j) 'количество молекул в ячейке
If i > 1 Then Cells(i, j) = 0
If i = im Then Exit Sub 'последний ряд
For k = 1 To km
r = frnd(j)
If i = 1 Then r = 0 'вариант перехода с входящего ряда
If Cells(i + 1, j + r) < d And Cells(i + 1, j + r) <> -100 Then
carry1(j + r) = carry1(j + r) + 1
ElseIf i > 1 Then
If Cells(i, j + r) = -100 Then r = 0
carry(j + r) = carry(j + r) + 1
End If
Next k
End Sub
Function frnd(j)
frnd = Int(3 * Rnd) - 1
If j = 1 And frnd = -1 Then frnd = 0
If j = jm And frnd = 1 Then frnd = 0
End Function
Листинг 3
Sub Paint(im, jm) vMax = 0 For i = 2 To im For j = 1 To jm If Cells(i, j) <> -100 Then c = Cells(19, 13).Interior.ColorIndex For k = 12 To 18
If Cells (i, j) <= Cells(k, 14) Then c = Cells(k, 13).Interior.ColorIndex Exit For End If Next k
Cells(i, j).Interior.ColorIndex = c If vMax < Cells(i, j) Then vMax = Cells(i, j) End If Next j Next i
If Range("N8") < vMax Then Range("N8") = vMax End Sub
гиону ячеек рабочего поля (канала), формируют (или очищают) ячейку препятствия (см. листинг 5).
8. Следующий стандартный макрос (в том же окне VBA) для отслеживания изменений на листе Excel формирует по заданному числу d верхние границы значений цветовой палитры и/или раскрашивает ячейки канала в соответствии с изменением границы диапазона или конкретного цвета (листинг 6).
Листинг 4
Sub Clear()
im = 20
jm = 11
For i = 2 To im
For j = 1 To jm
If Cells(i, j) <> -100 Then
Cells(i, j) = ""
Cells(i, j).Interior ColorIndex = xlNone
End If
Next j
Next i
Range("N5:N8") = ""
End Sub
Листинг 5
Sub Worksheet_BeforeDoubleClick(ByVal Target As Range, Cancel As Boolean) If Not Intersect(Target, Range("A2:K20")) Is Nothing Then Select Case Target Case -100
Target.ClearContents Case Else
Target = -100 End Select Cancel = True End If End Sub
Листинг 6
Private Sub Worksheet_Change(ByVal Target As Range) If Target.Address = "$M$8" Then ColMax = Range("M8") + 5 For i = 0 To 6
Cells(i + 12, 14) = Int(ColMax * i / 6) Next End If
If Not Intersect(Target, Range("M12:N19")) Is Nothing Then
Call Paint(20, 11) End If End Sub
И последние «штрихи к портрету» модели.
Чтобы скрыть отображение нулевых значений на рабочем поле, выделите ячейки А1:К20. В меню Формат ® Ячейки на вкладке Число в списке Числовые форматы выберите пункт Все форматы, а в поле Тип введите 0;-0;;@.
Для выделения «препятствия» в меню Формат ® Условное форматирование выберите в условии параметр Значение, опера-
цию сравнения Равно, значение «-100», а для формата ячеек выберите вкладку Вид и нужный цвет для заливки ячейки препятствия.
Модель готова к проведению экспериментов.
В приложении на CD представлена готовая модель.
Задания для самостоятельной работы
Смоделировать в Excel одномерный КА распределения нечетных коэффициентов в «треугольнике Паскаля». Выполнить задания (по степени сложности).
a) Сформулировать закон для данного КА. Определить количество нечетных чисел в указанной в задании строке m, считая вершину треугольника нулевой строкой.
b) Вывести общую формулу для количества k нечетных чисел в строке m «треугольника Паскаля» (условие олимпиадной задачи).
c) В предположении, что данный одномерный КА - круговой, определить, при каких длинах цепочка со временем «опустеет».
Указания к решению.
a) Пусть начальное состояние КА определяется первой строкой Excel. Следует выбрать значительное по протяженности (в ширину) поле ячеек в Excel, задать значение «1» примерно в центре входящего (первого) ряда поля. В условиях задачи нижеследующие строки Excel будут отражать изменение (в единицу «времени») значений ячеек предыдущей строки по следующему правилу (закону) КА:
ui+i, j = (\ j-1 + uh j+i) mod 2
Здесь индексы i, j - соответственно, номера строк и столбцов Excel, индекс i отражает временное изменение состояния КА. Аналогом сравнения по модулю 2 в Excel служит математическая функция ОСТАТ с делителем 2. Для удобства визуализации следует выполнить автоподбор ширины столбцов, форматирование ячеек поля, как описано выше.
b) Данный КА разворачивает во «времени» узор так называемой «салфетки
"DaHHuic "Xrf ра^барлчибаеЛ. бх «бремени» Лак набиваемой «самре&ки ферпиАсклы»...
Серпинского». Подсчитайте в свободном столбце Excel количество (сумму единиц) нечетных коэффициентов в каждой строке поля. Наглядно видно, что цепочка чисел результирующих сумм описывается рекуррентной формулой C(2)n+1 = C(2)n & 2*C(2)n (обозначения в условной записи: C(k) определяет цепочку из k чисел; знак & означает объединение (конкатенацию) цепочек; умножение на число 2 проводится для цепочки чисел поэлементно; n = 0, 1, 2, 3, ...).
Для заданной строки m, проведя ретроспективный анализ удаленности числа m от меньшего числа, являющегося степенью 2, получим общую формулу для задания b):
k = 2d , где d - количество единиц в двоичном разложении числа m.
с) Пусть в Excel в ячейке A1 задается длина КА (на рисунке 4 круговой одномерный КА длины 16 определяется ячейками C1:R1). При этом предполагаем, что стол-
и- ™ Я- ■ELn^:TOFEC;
ь. IF п|с J 1 J к H H I I r Ц •fi V *
1 i: 1 1
1 з ■ .........T-Ulll 1 1 ■ I 1 1 J
J 1 ■ —
t 1 i n 1
1 [ 1 1 ■
1 1 1 1 1 1
в г| ■ 1 1 1 1 4 1
?
1
ih 1 1
1 1
■1
J 1
■ I ■ ■ Шп i ■ i : 1 4
ГВ1В-В
Рисунок 4.
Рисунок 5.
бец B зациклен с конечным столбцом R (в примечании к ячейке указана соответствующая формула).
Для выделения совпадающих столбцов при различных значениях ячейки A1 (при «сворачивании» КА в кольцо) зададим в меню Условное форматирование заливку соответствующих ячеек, удовлетворяющих условиям, как на рисунке 5 (обратите внимание, что в задании формулы форматирования нет кавычек, которые Excel может поставить самостоятельно). Остальные значения ячеек, начиная со строки 2 и столбца 3 Excel, определяются строкой формул как для выделенной ячейки на рисунке 4. Та-
кое задание обходит создание циклических ссылок. Экспериментируя с различными значениями А1 и начальными единицами во входной строке, сделайте предположения относительно решения задания.
ВМЕСТО ЗАКЛЮЧЕНИЯ
Тенденция современного школьного образования состоит отчасти в том, чтобы наряду с преподаванием базовых знаний уменьшить временной разрыв между вновь созданными и используемыми в учебном процессе новыми технологиями обработки информации, научными воззрениями и исследованиями. Такая направленность способствует развитию широкого кругозора и интересов, формированию ранней профессиональной ориентированности учащихся, «мягкому» переходу к обучению в высших учебных заведениях, использованию возможностей распространенных пакетов прикладных программ.
Данная методическая разработка и предполагает изложение новых методов исследования сложных процессов в наглядной и доступной форме представления результатов исследования.
Литература
1. Паньгина Н.Н., Паньгин А.А. Статистическое моделирование: метод Монте-Карло // Компью терные инструменты в образовании, № 5, 2002. С. 30-43.
2. Gianluca Argentini. A first approach for a possible cellular automaton model of fluids dynamics. 2003, http://arxiv.org/ftp/cs/papers/0303/0303003.pdf
© Наши авторы. 2006 Ourauthors, 2006.
Паньгин Александр Викторович, инженер Центра информационных технологий, г. Сосновыш Бор.
