Моделирование броуновского движения на C99
Представляю вашему вниманию пятидесятую, юбилейную статью сайта. Этот юбилей я решил отметить не формулами и математическими выкладками, а картинками и видеороликом.
Мы будем заниматься моделированием броуновского движения программным способом.
Броуновское движение — это достаточно известное явление, которое изучается, если я не ошибаюсь, ещё в школе на уроках физики. Однако оно представляет интерес и с математической точки зрения. Недаром выдающийся русский советский математик Андрей Николаевич Колмогоров посвятил этому явлению статью "Броуновское движение и задача о блуждании по плоскости". Эта статья вошла в виде главы в книгу "Математика — наука и профессия", представляющую собой сборник статей, написанных Колмогоровым преимущественно для школьников и учителей математики и вышедшую в печать уже после смерти автора.
Вот что Андрей Николаевич рассказывает о броуновском движении в данной статье:
... в 1827г. ботаник Р. Броун открыл явление, которое по его имени названо броуновским движением. Броун наблюдал под микроскопом взвешенную в воде цветочную пыльцу. К своему удивлению, он обнаружил, что взвешенные в воде частицы пыльцы находятся в непрерывном беспорядочном движении, которое не удаётся прекратить при самом тщательном старании устранить какие-либо внешние воздействия, способные эти движения поддержать (например, вызвать движение самой воды под влиянием неравномерности температуры и т. п.). Вскоре было обнаружено, что это общее свойство любых достаточно мелких частиц, взвешенных в жидкости. Его интенсивность зависит только от температуры и вязкости жидкости и от размеров частиц (движение тем интенсивнее, чем температура выше, вязкость меньше, а частицы мельче). Каждая частица движется по своей собственной траектории, не похожей на траектории соседних частиц, так что близкие вначале частицы очень быстро становятся удалёнными (хотя могут иногда случайно вновь встретиться).
В этой же статье автор предлагает рассмотреть упрощённую математическую модель броуновского движения, которую и я собираюсь использовать. Вот что он пишет:
Основные черты броуновского движения частицы можно наблюдать уже на упрощённой модели блуждания частицы по плоскости, разделённой на квадратики. К таким упрощённым моделям при изучении более сложных явлений прибегают и в серьёзных научных исследованиях.
Будем считать, что наша частица перемещается отдельными шагами и за один шаг переходит из квадратика, в котором она находится в начале шага, в один из четырёх соседних квадратиков.
Далее Колмогоров исследует движение одной частицы в рамках предложенной им модели с точки зрения теории вероятностей и приходит к достаточно интересным выводам. Рекомендую читателю прочитать эту статью. Кстати, она является частью главы, посвящённой теории вероятности, которая является отличным введением в данный раздел математики для новичков.
Ну а мы будем использовать предложенную Андреем Николаевичем модель для создания ролика, демонстрирующего броуновское движение большого количества частиц.
Модель весьма простая; как видите, она изложена всего лишь в одном предложении; простой будет и программа, поэтому данную статью можно рассматривать как адресованную новичкам.
Предварительные замечания
Наше поле, состоящее из квадратиков (в дальнейшем мы будем называть их клетками), будет не бесконечным, а квадратным. В начальный момент времени некоторые клетки будут заполнены частицами. На каждом шаге каждая из частиц будет перемещаться на одну из четырёх соседних клеток. Состоянием поля будем называть распределение частиц по его клеткам в данный момент времени.
Начальное и последующие состояния поля будут преобразовываться в квадратные изображения. Каждой клетке будет соответствовать свой пиксель изображении. Он будет иметь белый цвет, если соответствующая ему клетка пуста. Если же клетка содержит одну или более частиц, то соответствующий ей пиксель будет чёрным.
Для построения изображений будем использовать библиотеку pgraph (предполагаем, что читатель с ней знаком хотя бы в общих чертах). Созданные изображения будем рассматривать как кадры анимации, которую мы сконструируем посредством программы VirtualDub.
Частицы неразличимы между собой, и перемещение каждой из них на каком-либо шаге не зависит от того, каким образом частица перемещалась ранее, а зависит только от её местонахождения непосредственно перед данным шагом. Поэтому мы можем не хранить информацию о каждой частице в отдельности. Для формирования состояния поля после некоторого шага нам нужно лишь знать исходное состояние (перед шагом). Поэтому будем фиксировать только состояния полей.
Для хранения состояний квадратных полей будем использовать двухмерные массивы. Каждый элемент массива будет соответствовать одной клетке поля и содержать количество частиц, находящихся в данной клетке. Нам понадобятся 2 таких массива. На каждом шаге из "старого" состояния поля, хранящегося в одном из массивов, мы будем формировать "новое" состояние, помещая его в другой массив. После завершения данного процесса на следующем шаге массивы уже поменяются ролями.
Для выбора направлений перемещений частиц на каждом шаге будем использовать генератор псевдослучайных чисел, реализованной в стандартной библиотеке языке C. Об этом мы ещё подробно поговорим в одном из последующих разделов.
Начальное состояние поля я предполагал сформировать, поместив большое количество частиц в одну-единственную клетку, чтобы наблюдать расширение минимально возможной области, занятой частицами. Однако я отказался от этого плана, в связи с одной особенностью нашей модели, которая заключается в следующем.
Предположим, что все клетки раскрашены точно так же, как поля шахматной доски. Пусть, для определённости, изначально частица находится в клетке чёрного цвета. Тогда после чётного количества шагов она всегда будет находиться в чёрной клетке, а после нечётного — в белой. А если сразу все частицы изначально располагаются в одной клетке, то в результате каждого шага все клетки, занятые частицами перед этим шагом, окажутся свободными от частиц. Таким образом, область, содержащая большое число частиц, будет выглядеть как своеобразное "решето". А в ходе шагов решёта будут переходить одно в другое, и всё это будет вызывать эффект мерцания.
Полагаю, что такое мерцание будет не очень хорошо смотреться, поэтому я решил разместить все частицы в 5 клетках, расположенных "крестом". Крест выбран из соображений симметрии. А чтобы после каждого шага количество частиц, расположенных в "белых" клетках, совпадало с количеством частиц, расположенных в "чёрных", изначально частицы предлагается разместить так, чтобы в центральной клетке их было столько же, сколько и в остальных 4-ёх в сумме. При этом частицы между 4-мя данными клетками должны быть размещены равномерно.
Пусть, например, общее число частиц равно N, где N делится на 8. Тогда частицы внутри креста будут распределены так, как показано на рисунке.
Этот "крест" будет располагаться примерно в центре квадратного поля.
А что будем делать с частицами, достигшими границ квадратного поля, и пытающимися выйти за его пределы? Пусть себе выходят. Вернуться в поле они, разумеется, уже не смогут. Будем считать, что поле — это квадратная столешница стола, а частицы, вышедшие за его пределы, просто упали на пол.
Структура программы
Программа состоит из главного файла brownian_motion.c и библиотечных файлов, входящих в библиотеку pgraph. В состав файла brownian_motion.c входят 3 функции — основная функция main() и 2 вспомогательные: rand4() и create_frame(). Первая из них генерирует псевдослучайные числа от 0 до 3, а вторая преобразовывает состояния полей в изображения и сохраняет последние в графических файлах.
Файл brownian_motion.c начинается с директив препроцессора #include и #define:
#include "pgraph.h" #define N 16000 //Количество частиц #define S 360 //Сторона квадрата #define F 4000 //Количество кадров
Как мы видим, к файлу brownian_motion.c подключается лишь один заголовочный файл, принадлежащий графической библиотеке. Макрос S содержит длину стороны квадратного поля, выраженную в длинах сторон клеток, его образующих. Разумеется, это значение совпадает с длиной и шириной соответствующего полю изображения в пикселях. Значение 360 взято, "с прицелом" на построение видеоролика: данная высота кадра является одной из стандартных.
Значения макросов N и F, смысл которых понятен из комментариев, подбираются эмпирически. Читатель может поэкспериментировать, переопределяя данные макросы.
Генерация псевдослучайных чисел — функция rand4()
В ходе выполнения программы постоянно будут генерироваться направления перемещений частиц. Всего возможны 4 варианта таких направлений, поэтому нам понадобится генератор псевдослучайных чисел, возвращающий при каждом запуске одно из чисел 0, 1, 2, или 3. Псевдослучайные числа из данного набора и будут возвращаться функцией rand4().
В основе этой функции будет лежать стандартная библиотечная функция rand(), объявленная в заголовочном файле <stdlib.h>. Эта функция генерирует целые равномерно распределённые псевдослучайные числа в диапазоне от 0 до RAND_MAX.
Компилятор, которым я пользуюсь, присваивает макросу RAND_MAX значение 32767. Таким образом, функцией rand(), фактически, генерируются 15-битовые псевдослучайные числа. Но нам для реализации одного перемещения одной частицы достаточно всего лишь 2-битового псевдослучайного числа.
Поступим мы следующим образом: будем расходовать числа, возвращаемые rand(), максимально экономно. Для этого мы будем использовать эти числа парами, "склеивая" каждый раз два 15-битовых числа в одно 30-битовое. Из полученного числа будем извлекать биты парами. Таким образом, мы сможем получить из одного такого числа 15 чисел от 0 до 3.
А теперь приведём код функции.
1.char rand4() 2.{ 3. static int count = 0; 4. static int random; 5. if (!count) 6. random = rand() << 15 | rand(); 7. char result = (random >> count) % 4; 8. if (count == 28) 9. count = 0; 10. else 11. count += 2; 12. return result; 13.}
Функция rand4() при каждом вызове возвращает одно псевдослучайное число от 0 до 3.
В строках 3 и 4 объявляются 2 статические переменные — count и random, которые сохраняют свои значения между вызовами функциями. Первая из них приобретает первоначальное значение 0 и увеличивается на 2 перед окончанием работы функции, за исключением случая, когда к окончанию работы её значение равно 28; тогда значение count обнуляется. Вторая переменная предназначена для хранения 30-битового псевдослучайного числа, получаемого "склейкой" двух чисел, генерируемых rand().
Как мы видим, новое значение переменная random приобретает только в случае, если значение count — нулевое (см. стр. 5, 6). Таким образом, random обновляется каждые 15 вызовов функции rand4(), начиная с первого.
В строке 7 содержимое random сдвигается на count разрядов влево, после чего 2 последних бита результата (они и образуют остаток от деления на 4) записываются в переменную result. Таким образом, в этой переменной оказываются биты числа, содержащегося в random с номерами count и count + 1 (биты нумеруем от младшего к старшему; нумерация начинается с 0).
В строках 8-11 значение count обновляется уже описанным выше способом, а в строке 12 из функции возвращается результат.
Заметим, что генератор псевдослучайных чисел от 0 до 3 получился весьма неплохим. Эксперимент показал, что при 10000 запусков функция rand4() возвращает числа 0, 1, 2 и 3 в 25,19%, 24,99%, 24,77% и 25,05% случаев соответственно. Для 100000 запусков данные показатели будут равны уже 24,85%, 25,01%, 25,01% и 25,13%. Наконец, если функцию запустить миллион раз, то показатели будут такими: 24,97%, 24,99%, 25,01% и 25,03%.
Видно, что в среднем доли выпадений 4-х чисел с увеличением количества вызовов функции rand4() приближаются к 25%. Математические расчёты это подтверждают: средние квадратические отклонения от 25% в рассмотренных трёх случаях составлют 0,151%, 0,098% и 0,023%. Они, как и следовало ожидать, уменьшаются с увеличением числа вызовов.
Создание кадра — функция create_frame()
Следующая вспомогательная функция — create_frame(). Вот её код:
1.void create_frame(int *p, image *img, int n) 2.{ 3. char filename[14]; 4. for (int i = 0, f = 0; i < S; i++, f += S) 5. for (int j = 0; j < S; j++) 6. set_color(img, i, j, p[f + j] ? BLACK : WHITE); 7. sprintf(filename, "results\\%04d.bmp", n); 8. save_to_file(img, filename); 9.}
Функция create_frame() принимает в качестве параметров адрес первого элемента двухмерного массива, содержащего информацию о состоянии поля, адрес объекта типа image и целое число. Функция записывает в объект типа image графическое представление состояния поля, после чего содержимое данного объекта сохраняется в графическом файле, который в дальнейшем становится одним из кадров анимации.
Код функции достаточно прост. Как мы видим, в двух циклах for, один из которых вложен в другой (см. стр. 4, 5), перебираются все элементы двухмерного массива. Если содержимое текущего элемента отлично от нуля (что означает наличие частиц в клетке, информация о которой хранится в элементе), то пиксель изображения, соответствующий элементу, окрашивается в чёрный цвет; если же значение элемента равно нулю (т. е. в случае отсутствия частиц в клетке), то пиксель окрашивается в белый цвет (стр. 6).
Обратите внимание на то, что соответствие между элементами массива устанавливается очень просто: индексы элемента, описывающего содержимое клетки, совпадают с координатами пикселя, отображающего данную клетку.
Далее генерируется имя графического файла с расширением bmp, в который записывается изображение (стр. 7). Часть имени, предшествующая расширению, представляет собой число, переданное в функцию, при необходимости дополненное нулями слева так, чтобы данная часть имени всегда состояла из 4-х символов. Такой способ наименования файлов обусловлен нашими дальнейшими планами по сбору всех файлов в ролик посредством VirtualDub.
Наконец, последняя операция, выполняемая функцией, — непосредственное сохранение изображения в файле (стр. 8), который создаётся в директории results, находящейся в корневом каталоге исполняемого файла (к моменту выполнения программы эта директория уже должна существовать).
Функция main()
Переходим к рассмотрению функции, выполняющей основную "работу" по моделированию броуновского движения. Ниже помещён её код.
1.int main() 2.{ 3. int square1[S][S], square2[S][S]; 4. int *p1 = (int *) square1, *p2 = (int *) square2; 5. for (int i = 0; i < S * S; i++) 6. p1[i] = p2[i] = 0; 7. square1[S / 2][S / 2] = N / 2; 8. square1[S / 2 + 1][S / 2] = square1[S / 2][S / 2 + 1] = N / 8; 9. square1[S / 2 - 1][S / 2] = square1[S / 2][S / 2 - 1] = N / 8; 10. image *img = create_image(S, S); 11. create_frame(p1, img, 1); 12. for (int k = 2; k <= F; k++) 13. { 14. for (int i = 0, f = 0; i < S; i++, f += S) 15. for (int j = 0; j < S; j++) 16. { 17. int ind = f + j, elem = p1[ind]; 18. if (elem) 19. { 20. int left = 0, right = 0, up = 0, down = 0; 21. while (elem--) 22. switch (rand4()) 23. { 24. case 0: 25. right++; 26. break; 27. case 1: 28. left++; 29. break; 30. case 2: 31. up++; 32. break; 33. case 3: 34. down++; 35. break; 36. } 37. if (j < S - 1) 38. p2[ind + 1] += right; 39. if (j > 1) 40. p2[ind - 1] += left; 41. if (i < S - 1) 42. p2[ind + S] += up; 43. if (i > 1) 44. p2[ind - S] += down; 45. p1[ind] = 0; 46. } 47. } 48. create_frame(p2, img, k); 49. int *t = p1; 50. p1 = p2; 51. p2 = t; 52. } 53. free(img); 54. return 0; 55.}
В строке 3 объявляем два двухмерных массива, каждый из которых состоит из S строк и S столбцов. В этих массивах будут храниться состояния полей. Сразу же объявляем 2 указателя и присваиваем им адреса данных массивов (стр. 4). Работать с массивами мы будем, в основном, через эти указатели. Во-первых, это позволит нам легко реализовать ситуацию, в которой массивы меняются ролями при переходе от одного шага к другому, поскольку такой обмен ролями будет эквивалентен обмену указателей значениями. Во-вторых, это ускорит обращения к элементам массивов. Данной теме, в частности, будет посвящена, статья о двухмерных массивах, которую я собираюсь написать.
Присваиваем всем элементам обоих массивов нулевые начальные значения (стр. 5, 6). Формируем в массиве square1 начальное состояние поля, распределяя все частицы между 5-ю клетками, образующими "крест" (стр. 7-9). Создаём динамически изображение, в которое будем последовательно помещать графические представления всех состояний поля (стр. 10) и тут же преобразуем начальное состояние поля в файл 0001.bmp (стр. 11) — первый кадр нашей будущей анимации.
Теперь в цикле for (стр. 12) формируем оставшиеся N - 1 кадров.
На каждой его итерации в двух циклах for, один из которых вложен в другой (см. стр. 14, 15) перебираем все элементы массива, адресуемого указателем p1. Значение каждого элемента сохраняем в переменной elem (стр. 17).
Далее вызываем функцию rand4() столько раз, сколько частиц находится в клетке, соответствующей текущему элементу массива (т. е. количество вызовов равно значению elem) и после каждого вызова увеличиваем значение одной из переменных left, right, up или down в зависимости от того, какое из 4-х чисел возвращено функцией (см. стр. 21-36). Эти переменные, как несложно догадаться из их названий, являются счётчиками, которые в ходе цикла while "подсчитывают" частицы, "перешедшие" на текущем шаге в каждую из четырёх соседних клеток (начальные значения этих переменных — нулевые (см. стр. 20)).
По окончании цикла while сумма значений, содержащихся в переменных left, right, up и down равна начальному значению переменной elem. Теперь мы увеличиваем на величины, равные значениям этих переменных, содержимое элементов массива, адресуемого указателем p2, соответствующих четырём соседним клеткам (см. стр. 37-44). При этом мы каждый раз удостоверяемся в том, что соседняя клетка принадлежит нашему квадратному полю (см. стр. 37, 39, 41 и 43). Если мы уверены в том, что частицы не выйдут за границы поля, то эти 4 строки, содержащие инструкции if, можно убрать, тем самым немного уменьшив время работы программы.
Мы не стали сразу в ходе цикла while модернизировать содержимое элементов массива, адресуемого указателем p2, а прибегли к использованию четырёх промежуточных переменных всё с той же целью — чуть-чуть сократить время работы программы, избавив компьютер от необходимости многократно вычислять адреса одних и тех же четырёх элементов массива. Экономия времени, скорее всего, весьма незначительна, но пусть она будет.
Далее обнуляем содержимое текущего элемента массива, адресуемого указателем p1 (стр. 45).
После обработки всех элементов массива, адресуемого p1, в массиве, адресуемом p2, содержится новое состояние поля. При этом все элементы первого из этих массивов — нулевые. Перед переходом к новой итерации самого внешнего цикла нам остаётся лишь сохранить изображение, соответствующее текущему состоянию поля и содержащееся во втором массиве, в графическом файле (стр. 48), после чего переключить указатель p1 на второй массив, а p2 — на первый.
А после построения всех графических файлов мы не забываем уничтожить изображение, созданное динамически в начале функции (стр. 53).
Результаты выполнения программы
После выполнения программы в директории results, расположенной в том же каталоге, что и исполняемый файл, появляются 4 тысячи файлов 0001.bmp, 0002.bmp, …, 4000.bmp, содержащих кадры нашей будущей анимации.
В первом кадре, как и ожидалось, изображён "крестик", состоящий и 5 пикселей. Далее, с увеличением номеров кадров, он начинает расти, "превращаясь" в почти идеальный ромб:
Далее ромб плавно переходит в круг, от которого отделяются отдельные частицы:
Постепенно круг увеличивается в размерах, теряя свою круглую форму, а окружающих его частиц становится всё больше и больше:
Этот процесс продолжается: пятно расплывается всё сильнее и сильнее, а внутри него начинает появляться всё больше и больше пустых областей. Мне это напоминает рой пчёл:
Наконец, пятно расплывается ещё больше, а некоторые частицы уже начинают выходить за пределы нашего поля:
Теперь переходим к созданию самой анимации. Как уже было сказано, для её построения я использую программу VirtualDub. Процесс создания анимации из отдельных графических файлов с помощью данного приложения уже был рассмотрен нами в предыдущих статьях (в этой и в этой), поэтому повторяться не буду. Скажу лишь, что я установил частоту кадров, равную 25 кадр./с., а для сжатия применил кодек Huffyuv, сжимающий без потерь. В итоге получился видеофайл длительностью 2мин. 40с., размером 376МБ.
Заранее было известно, кадры видеоролика, содержащие большое количество отдельных пикселей, цвета которых резко отличны от фона, будут испоганены Ютьюбом по причине высокой степени сжатия. Несмотря на это, я загрузил итоговый ролик на данный видеохостинг. Вот что из этого получилось:
Интересно, что первые секунд 15 сжатие вообще не заметно (если сильно не присматриваться). Видимо, на первых порах битрейта хватает, чтобы сжимать кадры с хорошим качеством. Дальше уже начинают появляться артефакты сжатия, но пока они ещё не мешают смотреть. Далее (примерно после 23-ей секунды) вокруг пятна появляется мерцающая "пелена", пока ещё терпимая. Но где-то уже к 35-й секунде эта пелена начинает сильно мешать. А потом всё уже идёт "вразнос". Один из последних кадров на Ютьюбе выглядит примерно так:
Так что лучше просматривать видеоролик в исходном качестве.
Заключение
Ну, вот и всё. Мы убедились в том, что с помощью достаточно простой программы можно смоделировать, пусть и очень грубо, весьма интересный физический процесс.