Fractals met Python
[oOO]
Fractals zijn al meer dan een eeuw geleden voor het eerst beschreven, maar ze zijn pas sinds de komst van de computer buiten de wiskundige wereld bekend geworden. In dit artikel laten we zien hoe je ze zelf met Python kunt maken.
In recente jaargangen van Pythagoras kwamen verschillende fractals aan bod: de Koch kromme, de Peano kromme, de Sierpiński driehoek en de Mandelbrot verzameling. De meeste van die fractals zijn al meer dan een eeuw geleden 'uitgevonden', maar ze zijn pas echt populair geworden — en daarmee bedoel ik dat je ze niet enkel in wiskundeboeken tegenkomt — met de komst van de computer die het toelaat om plaatjes te genereren van die fractals. En al lijkt zo'n fractal zoals de Sierpiński driehoek een lastige klus om te tekenen, voor een computer is het op één, twee, drie klaar. Nou ja, niet echt, want zo'n fractal vraagt, voor een perfect plaatje ervan, een oneindige nauwkeurigheid, en dat kunnen we nu eenmaal niet verwachten van een computer. We laten zien hoe we een goede benadering kunnen tekenen, en gebruiken daarvoor Python.
De driehoek van SierpiNski
Traditioneel wordt de Sierpiński driehoek gedefinieerd als een limiet. Je start met een volle gelijkzijdige driehoek. Zo'n gelijkzijdige driehoek kan je opdelen in vier gelijke maar kleinere gelijkzijdige driehoeken, en daar neem je de middelste van weg. Er blijven er dan drie over, en voor deze drie herhaal je het proces, tot in het oneindige (indien mogelijk). Je ziet het gebeuren in figuur 1, die de eerste vijf stadia toont.
[--- Figuur 1 ---]
Dit is niet zo eenvoudig na te bootsen met een computer. Maar er is gelukkig een betere manier om deze fractal te genereren. Bekijk het even anders: wat we in feite doen in de verschillende stappen is een gelijkzijdige driehoek in vieren delen en het middelste deel ervan weghalen. Nu kan je een transformatie definiëren die precies dat doet: een transformatie die de punten van de volledige driehoek herverdeelt over de drie kleine!
We gebruiken hierbij cartesische coördinaten in het vlak. We voeren daarvoor een assenstelsel in met een horizontale $x$-as en een verticale $y$-as. We tekenen de assen niet in figuur 2 maar we kiezen de oorsprong in het hoekpunt $A$. $A$ heeft dus coördinaten $(x,y) = (0,0)$. Het hoekpunt $B$ geven we coördinaten $(x,y) = (1,0)$, en dan heeft het punt $C$ als coördinaten $(x,y) = \left(\tfrac{1}{2} , \tfrac{\sqrt{3}}{2}\right)$ (waarom?).
[--- Figuur 2 ---]
De punten van de grote driehoek herverdelen we over één van de drie kleine. Dat kan eenvoudig met wat we noemen een lineaire transformatie, en die heeft wiskundig de volgende vorm:
$\begin{align*}
x' &=a \cdot x + b \cdot y + c\\
y' &=d \cdot x + e \cdot y + f
\end{align*}$
Het punt met coördinaten $(x,y)$ wordt hierdoor afgebeeld op het punt met coördinaten $(x’,y’)$. We noteren dit zo: $(x,y) \to (x’,y’)$.
De getallen $a$, $b$, $c$, $d$, $e$ en $f$ kiezen we zodanig dat de hoekpunten van de grote volledige driehoek afgebeeld worden op de overeenkomstige hoekpunten van de kleine.
We hebben duidelijk drie van die transformaties nodig. Om deze te vinden, is het voldoende om enkele punten in te vullen in de transformatieformules om zo de waarde van $a$, $b$, $c$, $d$, $e$ en $f$ te bepalen. Nemen we bijvoorbeeld de grijze driehoek, dan is dit wat de transformatie doet met de hoekpunten van driehoek $ABC$:
$(0,0)\to(-,0)$, $(1,1)\to \left(\tfrac{1}{2},0\right)$, $\left(\tfrac{1}{2}, \tfrac{\sqrt{3}}{2}\right) \to \left(\tfrac{1}{4}, \tfrac{\sqrt{3}}{4}\right)$
In dit geval worden beide coördinaten gewoon gehalveerd! Dit is het resultaat:
|
naar grijs: |
$\begin{align*} x’&= \tfrac{1}{2} \cdot x + 0 \cdot y + 0\\ y’&= 0 \cdot x + \tfrac{1}{2} \cdot y + 0 \end{align*}$ |
Ga dit zelf na door de coördinaten in te vullen in de transformatie. Op dezelfde manier vinden we ook de transformatie naar de groene en de rode driehoek: door de hoekpunten in te vullen. Dit zijn ze:
|
naar groen: |
$\begin{align*} x’&= \tfrac{1}{2} \cdot x + 0 \cdot y + \tfrac{1}{2}\\ y’&= 0 \cdot x + \tfrac{1}{2} \cdot y + 0 \end{align*}$ |
|
naar rood: |
$\begin{align*} x’&= \tfrac{1}{2} \cdot x + 0 \cdot y + \tfrac{1}{4}\\ y’&= 0 \cdot x + \tfrac{1}{2} \cdot y + \tfrac{\sqrt{3}}{4} \end{align*}$ |
Je kan zelf nagaan dat dit klopt!
Met deze transformaties kunnen we de punten die in de volledige driehoek liggen, herverdelen over de drie kleintjes. Daarvoor gebruiken we een van de drie transformaties. Nu blijkt dat we de fractal in kwestie als volgt kunnen opbouwen: start met een willekeurig punt in de driehoek, en pas dan daarop een van de transformaties toe, en herhaal dat tot in het oneindige. En laat daarbij het toeval beslissen welke van de drie transformaties je gebruikt, dus elk met kans $\tfrac{1}{3}$. We kunnen gewoon starten in het punt $(0,0)$. Een Python programma dat precies dit doet, maar niet oneindig doorgaat, staat hiernaast.
In figuur 3 zie je het resultaat na $50\,000$ stappen. Probeer zeker eens het aantal iteraties te verhogen!
[--- Figuur 3 ---]
Besluit: om dit te bereiken, hebben we de oorspronkelijke figuur opgesplitst in kleinere delen die allemaal dezelfde vorm hebben, en we hebben dan transformaties gedefinieerd die het geheel afbeelden op elk van de delen. Die 'zelfgelijkvormigheid' is een typisch kenmerk van fractals. Vele bekende fractals zijn van deze vorm, en je kan er nu zelf mee experimenteren. We geven er hier enkele.
import matplotlib.pyplot as plt
# We definiëren de 3 transfomaties
def f1(x, y):
return 0.5*x , 0.5*y
def f2(x, y):
return 0.5*x+0.5 , 0.5*y
def f3(x, y):
return 0.5*x+0.25 , 0.5*y+0.433
#
fs = [f1, f2, f3]
# We kiezen het startpunt
x, y = 0, 0
# We initialiseren de rij punten
xs = [x]
ys = [y]
# We herhalen het proces een aantal keer
for i in range(50000):
# waarbij we het lot laten beslissen welke
# transformatie we in elke stap toepassen
f = np.random.choice(fs, p=[1/3, 1/3, 1/3])
x, y = f(x, y)
xs.append(x)
ys.append(y)
# We tekenen de punten
plt.scatter(xs, ys, s=0.05, color='green', lw=0)
plt.axis('equal')
plt.axis('off')
plt.show()
Het tapijt van SierpiNski
(of de tweedimensionale spons van Menger) Het gaat hier om figuur 4 (links):
[--- Figuur 4 ---]
Je ziet duidelijk in figuur 4 (rechts) dat er in dit geval 8 transformaties nodig zijn, die telkens het volledige vierkant afbeelden op één van de acht kleinere.
[--- Figuur 5 ---]
Hierbij is het ideaal om het punt $O$ met coördinaten $(0,0)$ in het midden van het vierkant $ABCD$ te leggen (zie figuur 5), en de hoekpunten de volgende coördinaten te geven: $A(-1, -1)$, $B(1,-1)$, $C(1,1)$, $D(-1,1)$. Als je het zo doet, dan heb je voor het kleine lichtblauwe vierkantje links bovenaan:
$(-1,-1)\to \left(-1,\tfrac{1}{3}\right)$, $(1,-1)\to\left(-\tfrac{1}{3}, \tfrac{1}{3}\right)$, $(1,1)\to \left(-\tfrac{1}{3}, 1\right)$, $(-1,1)\to(-1,1)$.
Je hebt aan drie punten genoeg om alle onbekenden in je transformatie te bepalen. Je ziet voor dit geval direct dat $(0,0)\to\left(-\tfrac{2}{3},\tfrac{2}{3}\right)$. Met dit punt vind je zonder rekenwerk de waarden van $c$ en $f$ in de formules. Met twee van de hoekpunten vind je de rest.
Opgave 1Bepaal de acht transformaties. Tip: gebruik de symmetrie van de figuur. Als je bijvoorbeeld de transformatie gevonden hebt voor het vierkant linksboven, dan kan je die voor het vierkant rechtsonder hieruit afleiden door hier en daar van een plusteken een minteken te maken! |
Pas nu het Pythonprogramma aan en teken je resultaat.
De zeshoek van SierpiNski en de pentavlok van Dürer
De fractal in figuur 6 is een variatie op hetzelfde thema. Die in figuur 7 vindt zijn oorsprong in een werk van Albrecht Dürer (1525).
[--- Figuur 6 ---]
[--- Figuur 7 ---]
Opgave 2Je weet wat je moet doen! Figuur 7 is een hele uitdaging! |
Een fraCtal uit het boek FraCtals Everywhere van MiChael Barnsley
In dit boek uit 1988 kan je lezen waarom de bovenstaande methode werkt. De wiskundige stelling hierachter wordt de Collagestelling van Barnsley genoemd.
Opgave 3Schrijf een Pythonprogramma dat de fractal in figuur 8 genereert. Goed kijken is hier de boodschap! |
[--- Figuur 8 ---]
Opgave 4Vervang in je Python programma de transformaties door deze:
def f1(x, y):
return 0, 0.16*y
def f2(x, y):
return 0.85*x + 0.04*y, -0.04*x + 0.85*y + 1.6
def f3(x, y):
return 0.2*x - 0.26*y, 0.23*x + 0.22*y + 1.6
def f4(x, y):
return -0.15*x + 0.28*y, 0.26*x + 0.24*y + 0.44
en de kansen door deze: p=[0.01, 0.85, 0.07, 0.07] en bekijk wat dit als resultaat geeft! |
Meer hierover in het volgende nummer van Pythagoras