Зона кода

А давайте немного попрограммируем...

Построение броуновских мостов

C99Графические программы

В своё время мною была написана статья, в которой рассказывалось о моделировании броуновского движения. На очереди — ещё одно явление, связанное с британским ботаником Броуном. На этот раз мы будем рассматривать кривые, называемые броуновскими мостами.

Принцип построения броуновских мостов я взял из книги Седжвика, Уэйна и Дондеро "Программирование на языке Python: учебный курс". Эта книга уже служила ранее источником материалов для написания статей данного сайта; речь идёт об этой и этой статьях.

Статья адресована в первую очередь новичкам в программировании. Мы напишем на языке C99 простую программу, создающую изображения броуновских мостов, в основе которой будет лежать рекурсия. Для построения изображений будем использовать графическую библиотеку pgraph, предполагая, что читатель знаком с ней в общих чертах.

Также будет полезным знакомство читателя с понятиями случайная величина, и, более конкретно, случайная величина, распределённая по нормальному закону. Будем рассчитывать и на то, что читатель знает, что представляет собой генератор псевдослучайных чисел.

Понятие броуновского моста

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

Но нас сейчас не будут интересовать случайные процессы. Мы собираемся лишь рисовать кривые и поэтому примем то определение броуновского моста, которое можно сконструировать, руководствуясь книгой "Программирование на языке Python: учебный курс".

Итак, будем считать, что броуновский мост — это кривая, которая строится по следующему принципу.

Рассмотрим на координатной плоскости две точки с координатами (x1y1) и (x2y2), где x1 < x2.

На первом шаге построим точку (x3y3), где

x3 = (x1 + x2) / 2,
y3 = (y1 + y2) / 2 + δ1.

Таким образом, построенная точка — это середина отрезка, соединяющего точки (x1y1) и (x2y2), смещённая вверх или вниз (в зависимости от знака δ1) по оси ординат.

На втором шаге по тому же принципу строим точки  (x4y4) и (x5y5), где

x4 = (x1 + x3) / 2,
y4 = (y1 + y3) / 2 + δ2,
x5 = (x3 + x2) / 2,
y5 = (y3 + y2) / 2 + δ2.

Точно таким же образом действуем и на последующих шагах. На каждом i-м шаге строим 2i точек, абсциссы которых есть средние арифметические абсцисс двух соседних, уже построенных до этого шага точек, а ординаты отличаются от средних арифметических ординат соседних точек на δi. В результате n шагов таким образом оказываются построенными 2n − 1 точек (назовём их "внутренними").

Поговорим теперь более подробно о числах δi. Эти числа выбираются случайным образом. Каждое такое число представляет собой значение случайной величины, распределённой по нормальному закону, с математическим ожиданием, равным 0, и среднеквадратическим отклонением, равным σ / αi − 1.

Здесь σ — это начальное среднеквадратическое отклонение, применяющееся на первом шаге. А на каждом последующем шаге среднеквадратическое отклонение, как мы видим, получается из предыдущего делением на α. Так что среднеквадратичное отклонение на i-м шаге можно рассматривать как начальное отклонение, умноженное на "масштабирующий коэффициент", равный α1 − i.

Параметр α находится по формуле α = 2H. Здесь H — это так называемый показатель Херста.

Нужно иметь в виду, что для каждого фиксированного i все числа δi (в количестве 2i) — вообще говоря, разные. Я не использовал второй индекс в обозначениях этих чисел лишь для упрощения обозначений.

Итак, в результате выполнения n шагов мы имеем 2 исходные точки и 2n − 1 внутренних точек. Упорядочиваем эти 2n + 1 точек в порядке возрастания их абсцисс, после чего соединяем их последовательно отрезками, двигаясь от первой к последней. В результате получаем ломаную, состоящую из 2n звеньев и 2n + 1 вершин. Эта кривая и есть броуновский мост.

Если начало и конец ломаной фиксированы, то форма её определяется двумя параметрами: σ и H (разумеется, определяется не полностью, поскольку на неё также влияют случайные факторы).

От первого из них зависит, насколько сильно в среднем будет отклоняться каждая вершина ломаной от отрезка, соединяющего две соседние по отношению к ней вершины. Можно сказать, что σ влияет на "амплитуды колебаний".

Второй параметр влияет на "плавность" кривой. Чем больше H, тем более плавной выглядит кривая.

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

Но сначала разберёмся с генерацией псевдослучайных чисел, распределённых по нормальному закону.

О генерации нормально распределённых псевдослучайных чисел

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

Стандартная библиотека языка C99 содержит лишь функцию, генерирующую псевдослучайные целые числа, равномерно распределённые в диапазоне от 0 до значения макроса RAND_MAX. Так что нам придётся решать промежуточную задачу генерации нормально распределённых псевдослучайных чисел.

К счастью, такую задачу мы уже решили ранее. Все подробности можно узнать в предыдущей статье. Говоря коротко, мы создали функцию get_norm_rand_val(), при каждом вызове возвращающую псевдослучайное число. При многократных вызовах этой функции генерируемые ею числа оказываются распределёнными по стандартному нормальному закону, т. е. нормальному закону с параметрами 0 (математическое ожидание) и 1 (среднеквадратическое отклонение).

В свою очередь, из функции get_norm_rand_val() вызывается вспомогательная функция f_1(), подробно рассмотренная в упомянутой выше статье. Для корректной работы этой функции требуется ряд макросов и один константный глобальный массив.

Я решил, что разумно будет эти две функции, макросы и массив объединить в библиотеку norm-rand, которую в дальнейшем мы будем подключать к нашему коду. Она будет состоять из двух файлов: заголовочного norm-rand.h и файла с реализацией функций norm-rand.с.

Вот содержимое файла norm-rand.h:

#include <stdlib.h>
#include <math.h>

double get_norm_rand_val();

А вот — файла norm-rand.с:

#include "norm-rand.h"

#define sqrt_2 1.414213562373        //sqrt(2)
#define sqrt_pi_2 1.253314137316     //sqrt(pi / 2)
#define eps 1e-5

static const double val[] = {0.0,             //f_1(0.5)
                             0.253347,        //f_1(0.6)
                             0.524401,        //f_1(0.7)
                             0.841621,        //f_1(0.8)
                             1.281552,        //f_1(0.9)
                             1.880794};       //f_1(0.97)

double f_1(double x)
{
    if (x == 0.5)
        return 0.0;
    if (x < eps)
        return -4.3;
    if (x > 1 - eps)
        return 4.3;
    int sign = 1;
    if (x < 0.5)
    {
        x = 1 - x;
        sign = -1;
    }
    double xk1 = val[(int) (round(10 *x) - 5)], xk;
    double b = 1 - 2 * x;
    do
    {
        xk = xk1;
        double t = xk / sqrt_2;
        xk1 = xk - sqrt_pi_2 * exp(t * t) * (erf(t) + b);
    }
    while (fabs(xk - xk1) >= eps);
    return sign == 1 ? xk1 : -xk1;
}

double get_norm_rand_val()
{
    return f_1(rand() / ((double) RAND_MAX));
}

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

Структура программы

Программа будет состоять из файла main.c, а также подключаемых к этому файлу библиотек pgraph и norm-rand. Также нам понадобится стандартная библиотека stdio.h. Так что данный файл будет начинаться следующим образом:

#include <stdio.h>
#include "pgraph.h"
#include "norm-rand.h"

Далее поместим в файл main.c определения трёх макросов-констант:

#define sigma 30
#define H 0.8
#define N 9

Смысл макросов следующий:

  • sigma — начальное значение среднеквадратического отклонения (мы обозначали его σ);
  • H — показатель Херста (мы обозначали его H);
  • N — количество шагов, задействованных при построении броуновского моста (мы обозначали его n).

Затем объявим в файле глобальную переменную ordinates и глобальный массив deviations:

double *ordinates;
double deviations[N];

Вершины создаваемой нами ломаной будут представлять собой целые числа в диапазоне от 0 до 2n. А для хранения ординат будет динамически создан массив размера 2n + 1. Элемент массива с индексом i будет содержать ординату вершины, абсцисса которой совпадает с i. Глобальная переменная ordinates — это указатель на начальный элемент этого массива.

Для хранения среднеквадратических отклонений, задействованных на каждом шаге, будем использовать глобальный массив deviations. Элемент массива с индексом i будет содержать среднеквадратическое отклонение, использующееся при вычислении ординат вершин на шаге с номером n − i .

Далее в файле main.c будут находиться определения функций create_bridge() и main().

Построение вершин броуновского моста — функция create_bridge()

Рекурсивная функция create_bridge() занимается вычислением ординат вершин броуновского моста. Её код приведён ниже.

 1.void create_bridge(int i, int j, int n)
 2.{
 3.    int s = (i + j) / 2;
 4.    ordinates[s] = (ordinates[i] + ordinates[j]) / 2.0 + deviations[--n] * get_norm_rand_val();
 5.    if (n)
 6.    {
 7.        create_bridge(i, s, n);
 8.        create_bridge(s, j, n);
 9.    }
10.}

В качестве исходных данных функция получает координаты начала и конца броуновского моста, а также количество шагов, за которые нужно построить мост. Формальные параметры i и j содержат целочисленные абсциссы первой и последней вершин создаваемой ломаной. Ординаты этих точек содержатся в соответствующих элементах массива, адресуемого указателем ordinates (т. е. в элементах ordinates[i] и ordinates[j] соответственно). Количество шагов задаётся значением параметра n.

Таким образом, задачей нашей функции является вычисление ординат промежуточных вершин и заполнение ими элементов массива с индексами от i + 1 до j − 1 включительно. Будем предполагать, что значения i и j чётны и различны.

Код функции достаточно прост. Сначала вычисляем абсциссу "серединной" точки, т. е. точки, проекция которой на ось абсцисс равноудалена от проекций начальной и конечной вершин (стр. 3). Далее находим её ординату и помещаем в соответствующий элемент массива, адресуемого ordinates (стр. 4). Предполагается, что все элементы массива deviations к моменту вызова функции create_bridge() уже инициализированы. Обратите внимание, что в ходе вычисления ординаты мы уменьшаем на единицу значение переменной n.

Если после такого уменьшения значение n обнулилось, что означает, что текущий шаг был последним, то завершаем работу. В противном случае (см. стр. 5) с помощью двух рекурсивных вызовов вычисляем ординаты вершин, находящихся между первой вершиной и серединной точкой (стр. 7), а также между серединной точкой и последней вершиной (стр. 8).

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

Итак, создав предварительно массив, адресуемый указателем ordinates, состоящий из 2n + 1 элементов, заполнив первый и последний его элементы и вызвав функцию create_bridge() с фактическими параметрами 0, 2n и n, мы получаем в качестве результата её работы ординаты всех промежуточных точек, содержащиеся в элементах массива с индексами от 1 до 2n − 1 включительно.

Управление процессом построения броуновского моста — функция main()

Процессом построения броуновского моста будет заниматься функция main(). Она подготовит исходные данные для функции create_bridge(), вызовет её, после чего по полученным координатам вершин создаст само изображение броуновского моста.

Код функции main() приведён ниже.

 1.int main()
 2.{
 3.    const double alpha = pow(2, H);
 4.    deviations[N - 1] = sigma;
 5.    for (int i = N - 1; i; i--)
 6.        deviations[i - 1] = deviations [i] / alpha;
 7.    const int pow2N = pow(2, N), size = pow2N + 1;
 8.    ordinates = malloc(size * sizeof(double));
 9.    ordinates[0] = ordinates[pow2N] = pow2N / 4;
10.    create_bridge(0, pow2N, N);
11.    image *img = create_image(size, pow2N / 2);
12.    img->cur_pnt = (point) {0, ordinates[0]};
13.    for (int i = 1; i <= pow2N; i++)
14.        line_to(img, i, round(ordinates[i]));
15.    free(ordinates);
16.    save_to_file(img, "result.bmp");
17.    free(img);
18.    return 0;
19.}

Вычисляем значение константной переменной alpha, имеющей тот же смысл, что и параметр α (стр. 3). Заполняем среднеквадратическими отклонениями глобальный массив deviations (стр. 4-6). Инициализируем константные переменные pow2N и size, значения которых, равные 2n и 2n + 1 соответственно, нам пригодятся в дальнейшем (стр. 7). Динамически создаём массив, предназначенный для хранения ординат вершин броуновского моста, и присваиваем адрес его начального элемента указателю ordinates (стр. 8).

Инициализируем первый и последний элементы созданного массива (стр. 9). Обратите внимание на то, что в результате инициализации их значения оказываются равными. Таким образом, начало и конец нашего броуновского моста находятся "на одном уровне", т. е. равноудалены от оси абсцисс.

Что ж, оба массива уже готовы к вызову функции create_image(). Этот вызов мы и осуществляем в строке 10. Теперь массив, адресуемый ordinates, содержит ординаты вершин создаваемого броуновского моста (напомню, что абсциссами являются индексы элементов массива). Остаётся, используя эти данные, "нарисовать" сам мост.

Для этого создаём холст для рисования, т. е. объект типа image и присваиваем его адрес переменной img (стр. 11). Обратите внимание на то, что ширина холста для рисования совпадает с количеством вершин броуновского моста.

Устанавливаем начальную вершину моста в качестве текущей точки (стр. 12) изображения, адресуемого img, после чего рисуем на холсте саму ломаную (стр. 13-14). Теперь массив ординат, адресуемый ordinates, больше не нужен, и мы его удаляем (стр. 15).

Осталось записать изображение в файл result.bmp (стр. 16) и удалить его (стр. 17).

Результаты работы программы

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

Изображение броуновского моста, соответствующего показателю Херста, равному 0,8

Броуновский мост (показатель Херста 0,8)

Это изображение броуновского моста соответствует показателю Херста 0,8. Уменьшим его теперь в 10 раз, изменив соответствующий макрос в файле main.c следующим образом:

#define H 0.08

В результате компиляции и выполнения программы получаем следующую картинку:

Изображение броуновского моста, соответствующего показателю Херста, равному 0,08

Броуновский мост (показатель Херста 0,08)

Как мы видим, колебания стали более резкими. Авторы книги "Программирование на языке Python: учебный курс" сравнивают изображения подобного рода с "графиками биржевых курсов".

А теперь увеличим показатель Херста до 1:

#define H 1

Получаем следующий результат:

Изображение броуновского моста, соответствующего показателю Херста, равному 1

Броуновский мост (показатель Херста 1)

Теперь, наоборот, колебания становятся более плавными, кривая как бы "сглаживается". Авторам книги такого рода линия напоминает "горизонт в гористой местности".

Заключение

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

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

Скачать исходный код