Построение множества Жюлиа

Привет. Кипят страсти, конец года, сессии, дедлайны, новый год, а так же цензура проникает во все слои интернетов, что не может не печалить. Хабр уже не торт ^[1]. Просто хотелось написать, что я не согласен с таким подходом, но тогда бы меня просто забанили. Так что придется написать интересный контент. Хотя если забанят из-за предисловия к посту о множестве Жюлиа, ну что, тогда остатки торта стухли и шансов нет.

Итак, вернемся к теме поста. Я давно хотел немного больше узнать о комплексных числах, а не только то, что корень из минус единицы равен i. Особенно вызывали интерес фигуры имеющие фрактальную структуру, хотелось понять, что это значит, и как сделать такую визуализацию. Где то на полке стояла книжка по ТФКП ^[2], а так же закончился курс по комплексному анализу на курсере ^[3], и появилось немного свободного от работы времени. Приступим.

Итерация функции

Перед тем как написать алгоритм, придется изучить несколько базовых понятий. Начнем мы с понятия итерация функции f. Введем следующее обозначение n-ой итерации функции ^[4] f:

• где знак ○ — композиция функций ^[5]

За нулевую итерацию принимается тождественное отображение ^[6]:

Существует множество теорем о неподвижных точках ^[7] определенных типов отображений, вспомним что такое неподвижная точка некоторого отображения g, это решение следующего уравнения:

Аналогично, для итерации функции вводится притягивающая точка ^[8] — такая точка из области определения функции f, что последовательность значение итераций функции f сходится к некоторой неподвижной точке отображения f:

• где y — некоторая точка достаточно близкая к точке x

Множество всех y, удовлетворяющее предыдущему условию — называется предельным множеством точки x под действием итерации функции f.

Как правило, интерес представляет изучение последовательности функций образованных итерациями (подобную идею можно увидеть в методе Ньютона ^[9], искомое решение оказывается притягивающей неподвижной точкой построенного отображения, и потому может быть найдено как предел последовательности итераций.)

а так же последовательность значений вычисленных в некоторой начальной точке

Для комплексного мира визуализировать итеративную сходимость точки y к x под действием итерации функции f немного проблематично, все таки 4d, но для действительного мира картина всего 2d (красная траектория расходится, в то время как синяя сходится к неподвижной точке).

Итак, первый камушек, необходимый для понимания множества Жюлиа заложен, перейдем к следующему.

Квадратичный комплексный полином

Вообще множество Жюлиа строится для функции вида

• где p и q это комплексные полиномы

Мы же упростим задачу, и будем строить множество Жюлиа только для квадратичного полинома. Квадратичный комплексный полином в общем виде выглядит следующим образом:

Так же существует понятие унитарной (a = 1) центрированной (b = 0) формы квадратичного полинома

Введем следующее отображение φ, и найдем его обратное отображение:

Рассмотрим следующее выражение

Легко заметить, что последнее выражение напоминает квадратичный полином в общем виде p(z), давайте найдем такое c, что бы привести полученное выражение к квадратичному полиному общего вида

Таким образом мы показали, что при вышеприведенной замене c, полином общего вида можно записать как

Поведение квадратичного полинома под действием итераций

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

Не трудно показать, что

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

Хаотичное и нормальное поведение

Уже скоро можно будет дать определение множества Жюлиа. Мы будем рассматривать скорее интуитивные определения, нежели строго математические. Но сперва, рассмотрим понятие устойчивости. Решение некоторого дифференциального уравнения называется устойчивым, если поведение решения с условиями, близкими к начальным, не сильно отличается от поведения исходного решения. Не трудно догадаться, что вся соль во фразе не сильно отличается. Вообще выделяют различные типы устойчивости, например в определении множества Жюлиа участвует устойчивость по Ляпунову (это более ясное определение вытекает из теоремы Монтеля о компактном семействе функций ^[10]), но мы даже не будем вдаваться в суть различий устойчивостей, тут главное понять идею. Попробуем это проиллюстрировать. Для начала взглянем на интуитивное восприятие устойчивости, тут синяя и желтые траектории точек при небольшом изменении ведут себя почти одинаково, в то время как для красной траектории, небольшое изменение ведет совершенно в другую сторону.

Давайте рассмотрим следующий полином f(z) = z² — 1 + 0.213i, и построим траекторию для 100 итераций при z₀ = 0.3 + 0.2i и z₀ = 0.31 + 0.2i

В этом примере вы видите, что небольшое изменение в начальном условии вылилось в расхождение траекторий, такое поведение мы назовем хаотичным в некоторой окрестности начальной точки, т.е. поведение сильно зависит от небольших изменений в начальных условиях. В то время как поведение при котором сохраняется устойчивость, т.е. небольшие изменения в начальных условиях не сильно влияют на поведение в целом, мы будем называть нормальным, как в следующем примере. Рассмотрим следующий полином f(z) = z² — 1 + 0.2i (всего лишь чуть изменили константу c), и построим траекторию для 100 итераций при z₀ = 0.3 + 0.2i и z₀ = 0.31 + 0.2i

Множества Фату, Жюлиа и Мандельброта

• Множеством Фату полинома f(z) = z² + c называется такое подмножество множества комплексных чисел, для каждой точки которого, поведение функции под действием итераций является нормальным, т.е. траектория точек генерируемая итерациями, не сильно изменяется при изменении начальных условий в некоторой небольшой окрестности начальной точки. Обозначается F(f).
• Множеством Жюлиа полинома f(z) = z² + c соответственно называется такое подмножество множества комплексных чисел, для каждой точки которого, поведение функции под действием итераций является хаотичным, т.е. небольшие изменения в начальных условиях в некоторой небольшой окрестности начальной точки, значительно влияют на траекторию. Обозначается J(f)
• Множество Мандельброта — это такое множество параметров c полинома f(z) = z² + c, для которого множества Жюлиа связно ^[11].

Бассейн аттрактора

Итак мы узнали смысл множества Жюлиа, можно приступить к алгоритму построения такого множества. Но и для этого, придется сперва ввести некоторые другие определения.

• Аттрактором ^[12] называется множество состояний динамической системы, в направлении которых развивается система с течением времени. В нашем случае можно понимать под аттрактором подмножество неподвижных точек к которым сходятся последовательности итераций функции независимо от начальных условий.
• Бассейном аттрактора точки z функции f называется подмножество точек из окрестности z (обозначаемая как A(z)), такое что любая траектория начатая в одной из этих точек сходится к точке z
• Бассейном аттрактора к бесконечности A(∞) будем называть множество тех точек z, для которых траектории уходят в бесконечность, т.е. другими словами не сходятся к какой то неподвижной точке.

Рассмотрим A(∞) в контексте наших рассуждений для некоторого полинома f(z) = z² + c

Оказывается существует следующая теорема. Множество A(∞) (в вышеописанном контексте) является открытым ^[13], связным ^[11]и неограниченным ^[14]. Оно полностью содержится во множестве Фату. Множество Жюлиа совпадает с границей ^[15]A(∞), которая является замкнутым ^[16]и ограниченным ^[14] подмножеством всех комплексных чисел.

Построение множества Жюлиа

Напоминаю, что мы напишем алгоритм для построения множества Жюлиа квадратичной динамики, а не любой функции. Итак мы узнали что множество Жюлиа совпадает с границей ^[15]A(∞), которая является замкнутым ^[16]и ограниченным ^[14] подмножеством всех комплексных чисел, т.е. множество Жюлиа тоже замкнуто и ограничено. Обозначим следующим образом то подмножество комплексных чисел, которое генерируется итерациями функции f и остается ограниченным. Назовем такое множество заполненным множеством Жюлиа:

И последняя теорема которая нам понадобится звучит следующим образом:

• пусть дан некоторый квадратичный полином вида f(z) = z² + c
• обозначим
• тогда для некоторой точки и n > 0 выполняется следующие условие
• 
• т.е. если модуль n-ой итерации больше R, то итерации функции f расходятся, это эквивалентно тому, что точка z₀ входит в бассейн аттрактора к бесконечности, откуда следует, что z₀ не является точкой заполненного множества Жюлиа

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

PS: примерно тут любители дизайна и золотых сечений, мистики и оккультики, заметят, что при f(z) = z² + 1, порог R равен золотому сечению о_0

Алгоритм

1. Выбираем с для задания полинома f(z) = z² + с
2. Вычисляем R для заданного полинома f(z) = z² + с
3. Выбираем параметр maxIter для обозначения максимальной итерации, очевидно чем он выше, тем лучше точность, тем медленнее алгоритм
4. Генерируем массив цветов, всего maxIter штук, скажем от менее яркого к более яркому, мы будем обозначать цветом, то на сколько далеко точка расположена от множества Жюлиа
5. Для каждой точки вычисляем является ли она частью заполненного множества Жюлиа или нет, а так же номер итерации на которой порог был превышен
6. если |z| > R то используем первый цвет, далее используем тот свет, на каком номере итерации был превышен порог

private static XBitmap PlotJuliaSet(ComplexNumber c, int w, int h, int maxIter,
    double xMin = Double.NaN, double yMin = Double.NaN, double xMax = Double.NaN, double yMax = Double.NaN)
{
    double r = CalculateR(c);
    if (Double.IsNaN(xMin) || Double.IsNaN(xMax) || Double.IsNaN(yMin) || Double.IsNaN(yMax))
    {
        xMin = -r;
        yMin = -r;
        xMax = r;
        yMax = r;
    }            
    Logger.Instance.Log("R = " + r);
    double xStep = Math.Abs(xMax - xMin) / w;
    double yStep = Math.Abs(yMax - yMin) / h;
    XBitmap bmp = new XBitmap(w, h);

    IDictionary<int, IDictionary<int, int>> xyIdx = new Dictionary<int, IDictionary<int, int>>();
    int maxIdx = 0;
    for (int i = 0; i < w; i++)
    {
        xyIdx.Add(i, new Dictionary<int, int>());
        for (int j = 0; j < h; j++)
        {
            double x = xMin + i * xStep;
            double y = yMin + j * yStep;
            ComplexNumber z = new ComplexNumber(x, y);
            IList<ComplexNumber> zIter = SqPolyIteration(z, c, maxIter, r);
            int idx = zIter.Count - 1;
            if (maxIdx < idx)
            {
                maxIdx = idx;
            }
            xyIdx[i].Add(j, idx);
        }
    }
    for (int i = 0; i < w; i++)
    {
        for (int j = 0; j < h; j++)
        {
            int idx = xyIdx[i][j];
            double x = xMin + i * xStep;
            double y = yMin + j * yStep;
            ComplexNumber z = new ComplexNumber(x, y);
            bmp.SetPixel(w - i - 1, j, ComplexHeatMap(idx, 0, maxIdx, z, r));
        }
    }
    return bmp;
}

private static double CalculateR(ComplexNumber c)
{
    return (1 + Math.Sqrt(1 + 4*c.Mod))/2;
}

public static Color ComplexHeatMap(decimal value, decimal min, decimal max, ComplexNumber z, double r)
{
    decimal val = (value - min) / (max - min);
    return Color.FromArgb(
        255,
        Convert.ToByte(255 * val),
        Convert.ToByte(255 * (1 - val)),
        Convert.ToByte(255 * (z.Mod / r > 1 ? 1 : z.Mod / r))
    );
}

Результаты

Вот несколько в картинок 5к на 5к, сделанных с помощью этой функции

• f(z) = z² — 0.8 + 0.156i ^[17]
• f(z) = z² + 0.285 + 0.01i ^[18]
• f(z) = z² — 0.0085 + 0.71i ^[19]

А теперь давайте позумим множество Жюлиа для полинома f(z) = z² — 0.74543 + 0.11301i, всё множество содержится в шаре радиуса 1.50197192317588.

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

Автор: mephistopheies

Источник ^[20]