Моделирование в Maple: просачивание (перколяция)

Views: 454

В этой заметке я приведу пример использования Maple для математического моделирования на примере моделирования задачи двумерного  просачивания. Здесь мы будем использовать только функциональные средства Maple. В задаче описывается физическое явление, но в тоже время ее решение не требует глубоких знаний физики. Сформулируем эту задачу как физическую, например, в следующем виде.

Имеются два материала в сыпучем виде: идеальный проводник (далее на рисунках — красный) и идеальный диэлектрик (белый).  Мы перемешиваем эти материалы в некоторой пропорции и измеряем проводимость смеси. Нетрудно сообразить, что проводимость  будет зависить от пропорции. Имеется некоторая минимальная доля проводника, ниже которой проводимость смеси равна нолю.

Из формулировки задачи видно, что для ее моделирования достаточно обладать немногими знаниями из школьной геометрии и из теории вероятности.  Мы будем моделировать ее двумерный вариант, предполагая, что элементы материалов имеют форму одинаковых квадратов и при перемешивании прилегают вплотную друг к другу. Впрочем, мы обсудим обобщение нашего описания на трехмерный случай с элементами произвольной формы. Вот такая картинка иллюстрирует перколяционную решетку 19×19:

maple моделирование перколяции
Перколяционная решетка 19×19

Мы будем использовать  функции из двух графических пакетов Maple (plottools и plots), которые необходимо загрузить в ядро Maple дополнительно:

  restart:with(plottools): with(plots):

Для загрузки дополнительных пакетов используется оператор with(имя пакета) с именем пакета в качестве входящего параметра. По умолчанию естественно загружается не полный функционал Maple (для рационального использования памяти компьютера).

Затем мы определяем функцию rnd задающую равномерное распределение на интервале от 0 до 1:

rnd:= rand(0.0 .. 1.0):

Используется модуль Maple rand(a..b) создающий функцию-генератор псевдослучайных  равномерно распределенных по интервалу от a до b чисел. Отметьте что a и b мы записали в виде действительных чисел (с десятичной точкой). Функция  rand(0..1) создает генератор псевдослучайных равновероятных чисел 0 и 1.
Псевдослучайность  генератора rand(a..b) и как с ней бороться мы обсудим в следующей заметке.

Определяем функцую генерирующую квадратную перколяционную решетку, в качестве входных параметров служит размер решетки N и вероятность p того, что ячейка является проводящей. Функция возвращает 2×2 матрицу цветов ячеек. Матрицу здесь мы реализуем как список N списков N элементов.

MY_C_initialize:=proc(N,p) local i,j,r,MY_C:
#Создаем матрицу цветов, задавая всем ячейкам белый цвет.
MY_C:=[seq([seq(white,j=1..N)],i=1..N)]:
#Генерируем случайную перколяционную решетку
for i from 1 to N do
for j from 1 to N do r:=rnd():
if r < p then MY_C[i,j]:=red: else MY_C[i,j]:=white: fi:
od: od:
return(MY_C):
end proc:

Отметте использование знака # для комментария.
Название функции: MY_C_initialize.
Оператор присваивания :=, за которым следует объявление функции proc(N,P) со списком входных параметров (данных). В качестве входых параметров функций вообще говоря могут использоваться любые типы переменных Maple. Естественно, что для нашей функции N
это целое положительное число, а p число из отрезка [0,1].
Параметр local со списком локальных переменных. Использование оператора seq(.., i=n..m) для создание списка и использование квадратных скобок […] для определения упорядоченного множества (списка). Неупорядоченное множество определяется фигурными скобками, например, {1, f, ij,[12,black],{«ase»,1}}. Элементами множеств вообще говоря может быть все, что угодно, в том числе другие множества, функции, модули и т.п. Используется вложенный цикл по номеру ряда и колонки решетки. Оператор return(MY_C) возвращает матрицу соответствующую сгенерированной перколяции. Как и в других языках в функции может быть несколько операторов return(). Если такого оператора нет, или исполняемая часть функции закончилось другим оператором, то функция возвращает результат выполнения этого последнего опреатора. Завершается функция оператором end proc:

Определим функцию, которая создает матрицу чисел в соответствии с цветом соответствующей ячейки (в качестве входного параметра служит матрица цветов, сгенерированная например предыдущей функцией):

numberOFcolor:=proc(MY_C) local MY_N,N,i,j:
# Извлекаем размер решетки
N:=nops(op(1,MY_C));
# Создаем матрицу чисел первоначально
# отождествляя ее с матрицей 
цветов.
MY_N:=MY_C:
# Красному цвету будет соответствовать число 1,
# Белому — 0, и синему — 5.
# Двойной вложенный цикл по рядам и колонкам
# решетки (матрицы)

for i from 1 to N do
for j from 1 to N do
if MY_C[i,j]=red then MY_N[i,j]:=1: else MY_N[i,j]:=0: fi:
if MY_C[i,j]=blue then MY_N[i,j]:=5: fi:
od: od:
return(MY_N):
end proc:

Следующая функция обратна предыдущей, т.е. создает матрицу цветов в соответствии с числом соответствующей ячейки. Работает аналогично предыдущей функции.

colorOFnumber:=proc(MY_N) local MY_C,i,j,N:
N:=nops(op(1,MY_C));
MY_C:=MY_N:
for i from 1 to N do
for j from 1 to N do
if MY_N[i,j]=1 then MY_C[i,j]:=red: else MY_C[i,j]:=white: fi:
if MY_N[i,j]=5 then MY_C[i,j]:=blue: fi:
od: od:
return(MY_C):
end proc:

Наконец определяем главную функцию, которая ищет проводящий кластер. Используется следующий алгоритм:

  1. Присваивается число 5 всем проводящим ячейкам нижнего ряда.
  2. Движемся слева направо и (вложено) снизу вверх и присваиваем число 5 всем ячейкам, которые являются проводящими и имеют соседнюю ячейку с числом 5. Условия проверки различны для внутренней части, внешних колонок и верхних угловых ячеек.

Change_color:=proc() global MY_N,MY_C: local N,NN,NN2,i,j,MY_N_OLD,MY_F: N:=nops(op(1,MY_C)); MY_N_OLD:=MY_N:
# Рассчитываем сумму всех чисел в ячейках. Она понадобится для
# проверки роста проводящего кластера
NN:=add(add(MY_N[i,j],i=1..N),j=1..N):
# Затем следует двойной вложенный цикл
# в котором движемся слева направо и снизу вверх от второго ряда
for i from 1 to N do
if MY_N_OLD[i,1]=1 then MY_N[i,1]:=5: fi:
od:
for j from 2 to N-1 do for i from 2 to N-1 do
# Проверяем наличие соседних проводящих ячеек для ячейки
# внутренней части решетки, если такие есть
# присваиваем ячейке число 5
if (MY_N_OLD[i,j-1]+ MY_N_OLD[i,j+1]+MY_N_OLD[i-1,j]+MY_N_OLD[i+1,j] >= 5 and MY_N_OLD[i,j]=1) then MY_N[i,j]:=5: fi:
# Аналогичная проверка и присваивание
# для левой и  правой колонки

if (MY_N_OLD[1,j-1]+ MY_N_OLD[1,j+1]+MY_N_OLD[2,j] >= 5 and MY_N_OLD[1,j]=1) then MY_N[1,j]:=5: fi:
if (MY_N_OLD[N,j-1]+ MY_N_OLD[N,j+1]+MY_N_OLD[N-1,j] >= 5 and MY_N_OLD[N,j]=1) then MY_N[N,j]:=5: fi:
if (MY_N_OLD[i-1,N]+ MY_N_OLD[i+1,N]+MY_N_OLD[i,N-1] >= 5 and MY_N_OLD[i,N]=1) then MY_N[i,N]:=5: fi:
# Аналогичная проверка и присваивание
# для верхних угловых ячеек

if (MY_N_OLD[2,N]+ MY_N_OLD[1,N-1] >= 5 and MY_N_OLD[1,N]=1) then MY_N[1,N]:=5: fi:
if (MY_N_OLD[N,N-1]+ MY_N_OLD[N-1,N] >= 5 and MY_N_OLD[N,N]=1) then MY_N[N,N]:=5: fi:
od: od:
# Закончился двойной вложенный цикл.
# Теперь проверяем вырос ли в результате выполнения цикла
# проводящий кластер (рассчитывая сумму чисел ячеек NN2),
# если он вырос, то мы рекурсивно вызываем
# эту же функцию. Если рост остановился, тогда мы проверяем
# дорос ли кластер до верхнего ряда.
# Если дорос, то выводим на консоль «Yes», в противном случае
# выводим «No».
# Функция возвращает матрицу цветов, в которой проводящий
# кластер обозначен голубым цветом
NN2:=add(add(MY_N[i,j],i=1..N),j=1..N):
if NN2>NN then Change_color() else for i from 1 to N do
if MY_N_OLD[i,N]=5 then print(YES): MY_F:=colorOFnumber(MY_N): return(MY_F) fi:
od: print(NO): MY_F:=colorOFnumber(MY_N): return(MY_F) fi:
end proc:

Обратите внимание на то, что мы используем рекурсивный (вложенный) вызов функции. Например, факториал от целого n можно запрограмировать так:  fac:=proc(n)  if n=0 then return(1) else fac(n-1)*n; fi: end proc:  Впрочем, в Maple есть встроенная операция факториала в обычной нотации: n!

В случае использования рекурсивного вызова следует предотвращать возможную бесконечную рекурсию. В случае функции fac(n) мы ее предотвращаем задав значение fac(0)=1. В нашей функции Change_color() рекурсия прекращается как только кластер прекращает расти. Очевидно его рост ограничен размерами решетки.

Проверяем работу алгоритма, задав решетку размером N=19 и вероятность проводящей ячейки p=0.55. x1 и y1 — упорядоченные списки координат ячеек, используются для рисования решетки. У меня все выполняется очень быстро. Но лучше начать с небольшого значения N и ниже цикл генерации увеличивайте постепенно (начиная, например, с одного, а не с пяти как у меня). Для ускорения расчетов можно закоментировать вывод картинок или символом: #print(display(Sqs)); или взятием в кавычки: «print(display(Sqs))»;  

N:=19: p:=0.55: x1:=[seq(i-1,i=1..N+1)]:y1:=[seq(i-1,i=1..N+1)]:

В цикле пять раз генерируем перколяционную решетку MY_C, рисуем ее. Затем  используя процедуру (функцию) Change_color(MY_C) проверяем имеется ли проводящий кластер (снизу вверх) и рисуем максимальный проводящий кластер. Для рисования используем опрератор  display в качестве исходных данных которого служит Sqs — список списков (матрица) квадратов  (ячеек) задаваемых координатами левого нижнего и правого верхнего угла с соответствующими цветами.

for i from 1 to 5 do MY_C:=MY_C_initialize(N,p): MY_N:=numberOFcolor(MY_C): Sqs:=seq(seq(rectangle([x1[i],y1[j]],[x1[i+1],y1[j+1]],color=MY_C[i,j]),i=1..N),j=1..N):
# рисуем исходную решетку
print(display(Sqs)):
MY_F:=Change_color(MY_N,MY_C): Sqs:=seq(seq(rectangle([x1[i],y1[j]],[x1[i+1],y1[j+1]],color=MY_F[i,j]),i=1..N),j=1..N):
# рисуем решетку с максимальным проводящим кластером
print(display(Sqs));
od:

Заметьте использование двоеточия после od для подавления вывода промежуточных результатов и опреаторы print для принудительного вывода рисунков на консоль.

Результатом выполнения будет последовательность рисунков исходных решеток (см. выше) и результата поиска проводящего кластера (привожу для той же решетки, что изображена выше):

Maple, перколяция, percolation, моделирование
Проводящие снизу вверх кластера — синего цвета, проводящие ячейки не принадлежащие кластеру — краснного цвета, непроводящие ячейки — белого цвета. Имеется один кластер соединяющий нижний и верхний ряд решетки.

 

 

 

 

 

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *