В своей последней статье я попытался создать фрактал, порожденный простыми числами. Но он меня не очень устроил эстетически. Поэтому я решил воспользоваться zeta функцией Римана для создания фракталов.

Самый популярный фрактал, фрактал Мандельброта, создается итерацией функции:

X_{n+1} = X_n^2 + c

где X - комплексное число, а с - значение, которое для каждого пикселя берется согласно позиции этого пикселя. Если итерация уходит в бесконечность, то мы рисуем цвет согласно числу итераций до того, как функция превысила некоторое значение, а если такого не происходит, то мы рисуем пиксель черным. В качестве начального значения. X0. используется 0

Таким образом, при расчете фрактала Мандельброта для каждого пикселя используется разная функция. В отличие от такого подхода я в итерациях использовал одну и ту же функцию, от пикселя зависело только начальное значение, X0.

Первый эксперимент был с формулой

X_{n+1} = \frac{1}{Zeta(X_n)}

Особенно интересно поведение итераций вблизи нетривиальных корней, например первого (real = 0.5, img = 14.13..):

Желтая стрелка показывает, где находится первый корень. Цвета соответствуют числу итераций до abs(Z)>100. Сделаем zoom:

И еще:

И еще, мы сдвинулись вбок и корень уже ушел из поля зрения:

Мне кажется, этот фрактал можно назвать "бабочки Римана"?

Лирическое отступление о Julia

Я продолжал эксперименты с другими формулами, но процесс генерации миллиона точек (1000x1000) на Python занимал иногда до одного дня. Дело не только в медленности самого Python - подходящий пакет с zeta я нашел только в mpmath - вычисления с произвольной точностью. Это замечательный пакет, но очень медленный.

Я решил перескочить на другой язык, и мое внимание привлек язык Julia. Быстрый и полностью компилируемый, он позволяет вольготно обращаться с типами, почти как Python. Кроме того, он очень хорошо организует multi-threading. Например, перед циклом расчета фрактала, где все точки вычисляются независимо, достаточно написать

Threads.@threads for j in 0:iry
  for i in 0:irx
    ...

и все. Вот она, кнопка "сделать все зашибись". (да, да, не все так просто). В итоге то, что на Python занимало день, считалось несколько секунд. Я вышел на новый уровень и вот некоторые новые вычисления снова уперлись в продолжительность ожидания)

Я перебрал очень много вариантов, некоторые не давали интересных фракталов, зато итерация точки создавала интересные аттракторы:

Нетривиальные нули zeta

Мне все же хотелось найти функцию, которая бы показала все нетривиальные корни - ведь в этом случае картинка содержит бесконечно тонкую структуру всех простых чисел вообще. Такой функцией оказалась

X_{n+1} = X_n + exp(\frac{1}{Zeta(X_n)})

Сам по себе фрактал не был очень красивым, но у одного из его элементов была такая картина:

Как вы видите, расстояние между ножками соответствует расстоянию между корнями. Если увеличить острие хвостика второго корня:

то мы увидим ту же картину в более мелком масштабе.

Видео!

Но самой красивой оказалась формула

X_{n+1} = X_n + \frac{i}{Zeta(X_n)}

(условие выхода: Re(Xn)<-1.0)

Цвета соответствуют ближайшему номеру корня, у которого остановился расчет.

Сделаем zoom:

Можно заметить, что дерево 'растет' за счет все бОльших и бОльших корней. Это можно показать на видео:

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


  1. Matshishkapeu
    04.12.2022 23:24
    +1

    подходящий пакет с zeta я нашел только в mpmath

    Риманова Зета функция, не оно?

    from scipy.special import zeta


    1. Tzimie Автор
      05.12.2022 07:01

      Я знаю, но судя по описанию, это сумма ряда, а не аналитическое расширение.


    1. Tzimie Автор
      05.12.2022 07:30

      >>>zeta(3)
      1.2020569031595942
      >>> zeta(complex(0.5,14.1))
      Traceback (most recent call last):
      File "", line 1, in
      File "C:\...\site-packages\scipy\special_basic.py", line 2572, in zeta
      return ufuncs._riemann_zeta(x, out)
      TypeError: ufunc '_riemann_zeta' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''


  1. belch84
    05.12.2022 14:20
    +1

    Существует универсальное интегральное представление для дзета-функции, оно задается формулой Абеля-Плана
    image
    и работает для всех s, не равных 1

    С его помощью можно построить график для вещественных значений аргумента

    График дзета-функции для вещественных значений аргумента
    image