Pchełki Python: opowieść o najsłynniejszym żuku dwudziestego wieku

https://xpil.eu/aes 📋

Miałem cztery lata, kiedy matematycy po raz pierwszy podali definicję tego zbioru. Dwa lata później udało się zaprząc komputery do wykonania wizualizacji zbioru, a po kolejnych dwóch latach zbiorowi nadano oficjalnie nazwę "zbiór Mandelbrota", na cześć Benoit Mandelbrota, francusko-amerykańskiego matematyka (urodzonego, ciekawostka, w Polsce!). Benoit zajmował się fraktalami i opublikował na ich temat całkiem sporo fachowej literatury.

Z większej odległości zbiór Mandelbrota może się kojarzyć z żukiem. Ma duży tułów (z dwoma półokrągłymi półdupkami), mniejszą głowę oraz mnóstwo nóżek ze wszystkich stron. Trzeba oczywiście przymknąć oko na dodatkowe, coraz mniejsze głowy wystające tu i ówdzie, ale ogólnie podobieństwo do żuka jakieś tam jest.

Od strony matematycznej Żuk reprezentuje obszar stabilności pewnego całkiem prostego przekształcenia zespolonego.

Tak, wiem, ja sam na widok słowa "zespolony" przeważnie uciekam na drugi koniec, nigdy tak do końca nie przetrawiłem tego kawałka matematyki. Ale tu jest naprawdę prosto.

Otóż tak: bierzemy dowolną liczbę zespoloną P.

Następnie bierzemy liczbę zespoloną Z równą 0 (a konkretnie 0 + 0i)

I teraz tak: podnosimy Z do kwadratu, do wyniku dodajemy P i to co wyjdzie zapisujemy jako, uwaga, Z.

Po czym powyższy krok powtarzamy w nieskończoność.

Jeżeli po wykonaniu nieskończonej liczby kroków wyjdzie nam w wyniku jakaś konkretna liczba, to uznajemy punkt P za stabilny.

Punkt? Jaki, kurdę, punkt? Przecież P to liczba jest!

Liczba, owszem, ale zespolona, czyli składająca się z rzeczywistego iksa oraz urojonego igreka. A jak jest iks i igrek to już mamy dwie współrzędne, czyli punkt na płaszczyźnie.

Jeżeli natomiast po wykonaniu nieskończonej liczby kroków dostaniemy w wyniku nieskończoność (dodatnią bądź ujemną, bez znaczenia), uznajemy punkt P za niestabilny.

I teraz tak: jeżeli punkt P jest stabilny, rysujemy go na czarno, a jeżeli nie - na kolorowo. Kolor uzależniamy od stopnia niestabilności, czyli tego, jak szybko (w ilu krokach) dany punkt dotarł do nieskończoności.

Następnie bierzemy kolejny punkt P (jakiś inny) i powtarzamy cały proces od zera (dosłownie!). Czyli ustawiamy Z na 0 + 0i i zaczynamy wyliczać Z^2+P, a to co wyjdzie podstawiamy do Z i zapętlamy do nieskończoności. Powtarzamy powyższe dla wszystkich punktów na płaszczyźnie.

Czujny Czytelnik powinien teraz zakrzyknąć w myślach: "ten gość bredzi!"

Proszę bardzo, krzyczymy: "ten gość bredzi!"

Nie martwcie się. Nie obrażam się za takie coś, zresztą przywykłem już 😉

Po pierwsze, nie da się policzyć nieskończonej ilości punktów. Skończy nam się ilość pamięci w komputerze, potem Słońce zgaśnie i Wogle. No nie da się.

Po drugie nawet dla pojedynczego punktu nie da się wykonać nieskończonej ilości kroków. Czas się skończy, pamięć, Kosmos, takie tam...

Po trzecie, jak to "w ilu krokach punkt dotarł do nieskończoności"?!? No przecież że po nieskończonej, c'nie? Więc jak to jest z tym kolorem?

Skoro nie radzimy sobie za bardzo z nieskończonościami, zamiast tego musimy proces uprościć. Zbudować model.

Jak powiedział kiedyś jeden mądry facet:

Żaden model nie jest dokładny. Niektóre modele są użyteczne.

No dobra. Czyli tak: zamiast dla każdego punktu wyliczać nieskończenie wiele kroków, wyliczymy ich tylko pięćdziesiąt.

Czemu akurat tyle?

Bo jeżeli weźmiemy mniej niż pięćdziesiąt, przybliżenie będzie raczej niedokładne, a jeżeli weźmiemy dużo więcej, będziemy musieli długo czekać na wyniki. Pięćdziesiąt to dobry kompromis między jednym a drugim. Jeżeli ktoś ma szybki komputer (albo dużo czasu), może poeksperymentować z większymi wartościami.

Problem z nieskończoną ilością iteracji mamy załatwiony.

Drugi problem to jak po pięćdziesięciu krokach stwierdzić, czy dany punkt doszedł do nieskończoności czy nie?

Tu w sukurs przychodzi nam całkiem zaawansowana matematyka, której nie będę prezentował w szczegółach, bo sam jej do końca nie rozumiem, ale która mówi, że jeżeli po którymś z kroków |Z^2+P| przekroczy 2, to już będzie rosło w nieskończoność.

A więc zamiast sprawdzać, czy dotarliśmy do nieskończoności, wystarczy sprawdzić czy przekroczyliśmy dwójkę. Jeżeli tak, punkt jest z całą pewnością niestabilny. Jeżeli nie... No cóż. Jeżeli nie, to *raczej* jest stabilny, ale tak naprawdę to nie wiadomo. Czym więcej iteracji, czym większa dokładność, jednak stuprocentowej pewności nigdy mieć nie będziemy. Ale nie szkodzi. To tylko model. Przybliżenie w zupełności wystarczy.

Ostatni, trzeci problem to jak policzyć to draństwo dla wszystkich punktów na płaszczyźnie? Przecież jest ich całkiem sporo (nieskończenie wiele!), w dodatku są cholernie gęsto upakowane (nieskończenie gęsto!). Hmmm.

Rozwiązaniem tego problemu jest (a) wziąć pod uwagę tylko niektóre punkty (jeden piksel na ekranie to jeden piksel, gęściej nie narysujemy) oraz (b) wybrać sobie tylko fragment płaszczyzny. Ten najbardziej interesujący rozciąga się w poziomie od -2 do 0.8 oraz w pionie od -1.3 do 1.3.

Jeżeli do tego założymy, że nasz obrazek ma mieć rozdzielczość, dajmy na to, 600x600 pikseli, to będziemy musieli wykonać obliczenia dla 360 tysięcy pikseli, dla każdego z nich maksymalnie 50 iteracji, co daje razem... niecałe 18 milionów operacji.

Na papierze, z ołówkiem byłoby to raczej kłopotliwe. Ale komputery radzą sobie z takimi liczbami całkiem sprawnie.

Czas na kod:

from PIL import Image, ImageDraw
LimitIteracji = 50

# wymiary obrazka w pikselach
PikselePoziomo = 600
PikselePionowo = 558

# współrzędne kartezjańskie - granice obszaru
xStart = -2
xKoniec = 0.8
yStart = -1.3
yKoniec = 1.3


def mandelbrot(punkt):
    z, kroki = 0, 0
    while abs(z) < = 2 and kroki < LimitIteracji:
        z = z*z + punkt
        kroki += 1

    return kroki


# wybieramy HSV zamiast RGB żeby łatwiej było sterować kolorem
# # (w HSV o kolorze decyduje wyłącznie H,
# natomiast S i V mówią jak bardzo kolor jest nasycony i jasny)
obrazek = Image.new('HSV', (PikselePoziomo, PikselePionowo), (0, 0, 0))
kredka = ImageDraw.Draw(obrazek)

for x in range(0, PikselePoziomo):
    for y in range(0, PikselePionowo):
        # Skalujemy współrzędne pikseli ekranowych do współrzędnych rzeczywistych
        # oraz zamieniamy współrzędne rzeczywiste na liczbę zespoloną
        c = complex(xStart + (x / PikselePoziomo) * (xKoniec - xStart),
                    yStart + (y / PikselePionowo) * (yKoniec - yStart))
        # Wyliczamy liczbę iteracji
        m = mandelbrot(c)
        # kolor piksela zależy od liczby iteracji
        odcien = int(255 * m / LimitIteracji)
        # nasycenie ustawiamy na stałe na maksa
        nasycenie = 255
        # jasność maksymalna jeżeli piksel niestabilny, w przeciwnym razie zero (czarny)
        jasnosc = 255 if m < LimitIteracji else 0
        # stawiamy kropkę
        kredka.point([x, y], (odcien, nasycenie, jasnosc))

# po narysowaniu wszystkich pikseli, zapisujemy obrazek na dysk
obrazek.convert('RGB').save('output.png', 'PNG')

Powyższy kod uruchamiamy (może zajść potrzeba uprzedniego zainstalowania biblioteki Pillow - po szczegóły odsyłam do Google) a następnie zaglądamy do folderu ze skryptem - powinien tam pojawić się plik o nazwie output.png, z takim oto obrazkiem:

Piękny Żuk, prawda?

Po zwiększeniu limitu iteracji do 100 oraz podwojeniu wymiarów w pionie i w poziomie czeka się ciut dłużej, ale efekt jest jeszcze piękniejszy:

Przy 500 iteracjach oraz kolejnych zwiększeniu wymiarów 2x w każdym wymiarze efekt jest taki:

Z bliska wygląda to tak:

Można teraz, manipulując współrzędnymi oraz ilością iteracji, robić sobie dowolne zbliżenia - efekty są całkiem interesujące, ale o tym może już kiedy indziej.

Czytasz dalej?

Kolejna ciekawostka dotycząca Żuka jest taka, że są w nim "zakodowane" dwie bardzo słynne stałe matematyczne: Pi oraz Fi.

Pi zakodowana jest w sposób dość niezwykły: otóż jeżeli ustawimy się w punkcie (-3/4, 0) oraz zrobimy "zoom" o powiększeniu 10^x - krotnym (dajmy na to, dziesięć milionów razy, czyli 10^7), żeby narysować każdy piksel w tym powiększeniu, będziemy musieli wykonać Pi * 10^x iteracji (w tym przypadku: 31415928 iteracji).

Jeśli natomiast chodzi o Fi (złoty podział, te sprawy), to okazuje się, że sumy częściowe ilości "antenek" na okrągłych "głowach" Żuka kolejnych poziomów układają się w ciąg Fibonacciego.

Na koniec kilka obrazków z naprawdę dużym zoomem (pobrane z Wikipedii, a wygenerowane przez Wolfganga Beyera za pomocą programu UltraFractal 3):

Disklajmer: kod do dzisiejszej Pchełki ukradłem stąd. Trochę go spolszczyłem dla czytelności.

Dopisane półtora roku później: Żuk Mandelbrota jest też głównym tematem hipotezy MLC (skrót od Mandelbrot Locally Connected) mówiącej, że jest on lokalnie spójny (na nasze, w uproszeniu, że nie jest porozrywany na mniejsze kawałki, czyli z każdego miejsca na Żuku da się dotrzeć w każde inne miejsce bez opuszczania Żuka).

https://xpil.eu/aes 📋