0,2
рейтинг
18 декабря 2013 в 16:01

Разработка → Построение множества Жюлиа

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

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

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


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



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


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


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



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


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

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



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



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



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

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



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



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


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



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



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





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



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



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



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


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



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



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

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



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



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



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



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



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


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



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

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


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



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

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



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



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


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

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

Алгоритм


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


Процедура генерации множества Жюлиа
/// <summary>
/// Generate bmp with Julia set
/// </summary>
/// <param name="c">parameter of square dynamics</param>
/// <param name="w">width of bmp</param>
/// <param name="h">height of bmp</param>
/// <param name="maxIter">max iterations of function</param>
/// <param name="xMin">window in complex plane</param>
/// <param name="yMin">window in complex plane</param>
/// <param name="xMax">window in complex plane</param>
/// <param name="yMax">window in complex plane</param>
/// <returns></returns>
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 IList<ComplexNumber> SqPolyIteration(ComplexNumber z0, ComplexNumber c, int n, double r = 0)
{
    IList<ComplexNumber> res = new List<ComplexNumber>();
    res.Add(z0);
    for (int i = 0; i < n; i++)
    {
        if (r > 0)
        {
            if (res.Last().Mod > r)
            {
                break;
            }
        }
        res.Add(res.Last() * res.Last() + c);
    }
    return res;
}

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) = z2 — 0.74543 + 0.11301i, всё множество содержится в шаре радиуса 1.50197192317588.













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

Павел Нестеров @mephistopheies
карма
97,7
рейтинг 0,2
data scientist

Самое читаемое Разработка

Комментарии (42)

  • +4
    Стоит еще добавить тег «фракталы».
    Спасибо Вам за эту красоту!
    • +1
      точно -)
      • +11
        чорт вы успели отредактировать, а я нет, отвечу тогда вторым коментом

        Спасибо Вам за эту красоту!


        спасибо Эйлеру, Риману, Коши, Вейрштрассу, Фату, Жюлиа, Мандельброту -)
  • +7
    Как бы ненароком между вступлением к статье и тегами в конце читателю предоставляется возможность окунуться в красоту мира математики.
    Спасибо, действительно интересно!
  • +9
    Не занудства ради, но всё же: не «корень из минус единицы равен i», а «квадрат i равен -1».
    За статью, большое спасибо.
    • +3
      будем считать, что до углубления в ТФКП я этого не знал, а теперь вот знаю что именно i2 = -1 -)
    • +9
      Более интригующе
      корень из минус единицы равен плюс-минус i.

      И еще хлеще
      корень из минус единицы равен нулю плюс-минус i.
  • +5
    Меня всегда удивляли такие люди, в хорошем смысле. Теги ударные!
  • –20
    пост странное впечатление производит — как-будто хорошая книжка измазана говном
    я не понимаю, почему у многих людей так тяжело укладывается в голове, что правила нужно соблюдать? к чему эта позиция мелкого пакостника?
    • +3
      глубоко -)

      так а какие правила?
      • –6
        habrahabr.ru/info/help/rules/

        «Хабр — не для политики. На сайте крайне не приветствуются дискуссии на политические темы в любом их проявлении.»
        • +15
          Хабр — не для правонарушителей… Говоря проще — запрещено обсуждать правила сайта...

          habrahabr.ru/info/help/rules/

          воу воу парень палехче, ты что правила сайта обсуждаешь? да ты анархист
        • +4
          Вы удивитесь, но политику много лет никто не трогал, пока политика не начала трогать нас.
          • –9
            по-моему, в интернетах достаточно профильных площадок, чтобы поделиться своими печалями
            • +4
              так скоро не останется ни одной -)
            • +8
              «Площадок»? Да мне плевать на площадку. Интересующую тему, если понадобится, я могу обсуждать хоть в фейсбуке, хоть на форуме, хоть вообще по емылу. Ценность представляет не площадка, а сообщество.
  • +1
    не туда, нада сюда было habrahabr.ru/post/206516/#comment_7114136
  • +3
    Можно просто траурную ленточку в углу всех картинок и разъяснение позиции про скорбь о свободе слова и информации в тегах.
  • +6
    Впечатление от поста испортил тег. Ну е-мое, ну как будто на заборе «ху*й» написали, как дети.
    • НЛО прилетело и опубликовало эту надпись здесь
  • +3
    Спасибо! Отличный пост.
  • +3
    Бассейном аттрактора к бесконечности A(∞) будем называть множество тех точек z, для которых траектории уходят в бесконечность, т.е. другими словами не сходятся к какой то неподвижной точке.

    Эти понятия не эквивалентны. В общем случае бывает, что последовательность не стремится ни к бесконечности, ни к конкретной точке (а, например, осциллирует между двумя точками). Возможно, конкретно в данном случае траектория fn(z) всегда имеет предел (конечный или бесконечный), но это нужно отдельно доказывать.
    То же самое касается выкладок в разделе «построение множества Жюлиа».
    • 0
      да точно, про осцилляторы совсем забыл упомянуть, спасибо за верное замечание, сейчас добавлю это уточнение в определение A(∞)
      • 0
        Так ведь случаи не исчерпываются тремя — «сходится к точке», «осциллирует между точками», «уходит в бесконечность».
        Что, если последовательность осциллирует между неподвижной точкой и уходящей в бесконечность траекторией?
        • 0
          т.е. дойдет до бесконечности и вернется обратно? -)

          «сходится к точке», «осциллирует между точками»

          на самом деле суть одна, все это периодические точки, а в частном случае неподвижная точка — это периодическая с периодом равным единице
          • 0
            Художник бы смонтировал анимацию, а я для пояснения могу предложить
            только это
            время
            |  *
            |   *
            |  *
            V     *
               *
                    *
               *
                      *
               *
                        *
               *
                          *
            положение --->
            

            Т.е. ни с какого момента траектория не выйдет из окрестности неподвижной точки; но и не сходится к ней.
  • +1
    Легко заметить,..

    Как говорил мой научный руководитель, подобные обороты — явный признак того, что читатель «хрен повторит выводы автора». Но не в этом случае, т.к… и предыдущие выкладки я пониманию условно. С уважением отношусь к понимающим.
    • 0
      ну там где я это упомянул действительно легко заметить, что первые два члена это из общей квадратичной формы, а три слагаемых являются свободными
      image

      а вот тут я чот не рискнул писать, «легко заметить», что "Множество A(∞) (в вышеописанном контексте) является открытым, связным и неограниченным.".

      я в последнем случае даже не пытался найти доказательство =)
  • 0
    Круто!
    Кстати, меня всегда терзала мысль, что фракталы придумали производители обоев, а не Бенуа Мандельброт :)
    • 0
      Хочу фрактальные обои, которые можно бесконечно увеличивать.
      • +2
        Фрактальные обои — беда для аутиста: компьютер хотя бы выключить можно, а от обоев не спрячешься.
        И смерть для параноика. Из каждой розетки за ним будет следить бесконечное множество розеток, из которых за ним будет следить…
  • +1
    Мы в школе каждые полгода-год изучали новый язык. На каждом языке писали одни и те же программы. Программа для навигации по множеству Жюлиа или Мандельброта была одной из них.
    • +6
      крутая школа у вас была, мы в школе каждые пол года изучали новую карточную игру и новый пивной ларек -)
      • +1
        Одно другому не мешает :)
      • 0
        Сим-сити, что ли? Там и карты, и ларьки, и стройка, по которой можно мышкой полазить
  • 0
    Fractal Orgy — www.youtube.com/watch?v=zvHhnGjSmr8 и исходный код www.shadertoy.com/view/4dsGWB
  • 0
    Спасибо за статью.

    Она мне очень кстати пришлась, поскольку «легла» на мои нынешние рассуждения над серьезными наборами экспериментальных данных.

    Думаю, может в окрестностях этого поста я найду кого-нить, интересующегося изложенным выше с плавным переходом к «E.M.»-алгоритму и его применению к реально большим данным?..

    А вдруг!

    • 0
      ух ты, интересна связь этого поста и ем алгоритма =) мои контакты тут pavelnesterov.info/ пишите, пообщаемся
  • 0
    Связь, как я ее вижу:

     f(z)=z^2+c,

    А теперь представьте, что c=N(0,s^2) (хотя бы для примера, или что-нить случайное, но дискретное).

    А вместо вашего «i» — давайте «по взрослому» не будем замыкаться в плоскости, а возьмем многомерное (бесконечномерное?...) пространство.

    И скажите мне — не будет ли это все похоже на траекторию из оценок е.м.? Или траектория оценок е.м. похожа на «это»…
    Я уже пол-года или больше «болею» этим вопросом.
    Сам я его разрешить не в состоянии.

    И при каких s они будут куда-нить плавненько сходиться?
    А есть ли каки-нить хитрые условия, когда хоть парочка уровней фрактала" будет просматриваться невооруженным глазом.

    Еще раз говорю — я могу заблуждаться.
    Тогда просто считайте — это преддверие НГ.
    А вообще у меня здесь лето и последний час последнего рабочего дня перед отпуском :)

    Всех хабровчан с Новым Годом!
  • 0
    Эх, на первых курсах универа писал программу, где строились разные варианты множеств Мандельброта, Жюлиа и Фату. С разными модификациями. И даже на университетской конференции про них рассказывал… В общем, спасибо за ностальгию! =)

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