Apollonische Netze mit Python

Wed 2014-08-27

Translations: en

Apollonisches Netz mit D3-Symmetrie.

Abbildung 1: Apollonisches Netz mit D3-Symmetrie.

Apollonische Netze sind Fraktale, die aus drei sich gegenseitig berührenden (»apollonischen«) Kreisen berechnet werden. Aus diesen ersten drei Kreisen lassen sich rekursiv weitere Kreise erzeugen, die einen umschließenden Kreis ausfüllen.

Den Algorithmus dazu habe ich Python als ein Kommandozeilenprogramm, dass die Kreise als svg-Datei ausgibt, implementiert. Im Folgenden erläutere ich die Mathematik dahinter.

Eine Online-Version zum selbst spielen gibt es auch!

Da dieser Artikel etwas länger geworden ist, gibt es hier zunächst ein

1   Wie es funktioniert

Um ein apollonisches Netz zu erzeugen, können wir wie folgt vorgehen:

Am Anfang stehen drei sich gegenseitig von außen berührende Kreise. Jetzt suchen wir einen vierten Kreis, der die drei gegebenen berührt. Man sieht leicht ein, dass es dazu genau zwei Möglichkeiten gibt: Entweder liegt der vierte Kreis im Zwickel zwischen den dreien und berührt sie außen, oder er umschließt die drei und berührt sie von innen.

Drei sich berührende Kreise, zwei Möglichkeiten für einen vierten.

Abbildung 2: Drei und vier sich gegenseitig berührende Kreise

Links: Drei sich paarweise berührende Kreise mit unterschiedlichen Radien.

Rechts: Zwei Möglichkeiten für einen vierten Kreis, der alle drei gegebenen berührt. Einschließend (grün) oder ausschließend (pink).

Zunächst benötigen wir den einschließenden Kreis.

Es existieren zwei Gleichungen, die die Krümmungen und Mittelpunkte von vier paarweise berührenden Kreisen in Beziehung zueinander setzen. Das bedeutet: Kennt man schon drei Kreise, so ist der vierte leicht durch Einsetzen auszurechnen.

Das ist aber noch nicht alles: Die besagten Gleichungen sind quadratisch, d.h. man kann den Satz von Vieta auf sie anwenden. Das ist recht nützlich, da quadratische Gleichungen meist zwei Lösungen haben (in unserem Fall den einschließenden und den ausschließenden Kreis) und man mit dem Satz von Vieta ohne viel Aufwand die zweite Lösung ausrechnen kann, wenn man die erste schon kennt.

Auf unsere vier Kreise angewandt bedeutet das:

  1. Hält man die drei inneren Kreise fest und betrachtet den einschließenden als die erste Lösung, so liefert Vieta den ausschließenden, kleinen Kreis zwischen den dreien.
  2. Betrachtet man der Reihe nach einen der drei inneren Kreise als die erste Lösung und hält die anderen drei fest, so erhält man je einen weiteren neuen Kreis.
Vier Möglichkeiten, den Satz von Vieta anzuwenden

Abbildung 3: Vier Möglichkeiten, den Satz von Vieta anzuwenden

Die grauen Kreise werden festgehalten, der grüne ist die erste, bereits bekannte Lösung.

Der pinke Kreis ist die zweite Lösung, die der Satz von Vieta liefert.

Insgesamt erhalten wir so vier neue Kreise.

Vier neue Kreise wurden zum Apollonischen Netz hinzugefügt.

Abbildung 4: Vier neue Kreise wurden zum Apollonischen Netz hinzugefügt.

Dieser Prozess kann jetzt auch mit den neu hinzugekommenen Kreisen fortgeführt werden: Finde vier paarweise berührende, betrachte der Reihe nach jeden als die erste Lösung und berechne mit Vieta die zweite.

Neue Kreise nach dem zweiten Schritt der Iteration.

Abbildung 5: Neue Kreise nach dem zweiten Schritt der Iteration.

Wenn uns das apollonische Netz vollständig genug erscheint brechen wir die Berechnung ab, legen die Farbe eines Kreises anhand seines Radius fest und speichern als svg.

Soweit ein grober Überblick. Jetzt kommen die Details. Wer lieber direkt selbst ausprobieren möchte: Ausprobieren!.

2   Die Mathematik

Bevor es an die Implementierung geht, müssen mehrere kleinere Probleme gelöst werden.

2.1   Warum es funktioniert: Vieta

Für alle, die nicht mehr so fit in der Mittelstufenmathematik sind, hier eine kurze Wiederholung.

Betrachtet man die quadratische Gleichung \(f(x)=x^2 + px + q\) mit den Lösungen \(x_1\) und \(x_2\) von \(f(x)=0\), dann besagt der Satz von Vieta:

\begin{equation} \label{eq:vietaplus} x_1 + x_2 = -p \end{equation}

\begin{equation} x_1\cdot x_2 = q \end{equation}

Wenn wir \(x_1\) schon kennen, so können wir \(x_2\) trivialerweise berechnen als \(x_2 = -p - x_1\).

Beachte, dass wir nicht \(x_2=q/x_1\) benutzen, da \(x_1\) auch gut \(0\) sein könnte und uns der ganze Krempel dann um die Ohren fliegt.

2.2   Darstellung von Kreisen

Es gibt viele Möglichkeiten, Kreise zu repräsentieren. Zwei davon sind

  1. Mittelpunkt und Radius
  2. Mittelpunkt und Krümmung

Der Mittelpunkt ist ein Punkt in der Ebene oder eine komplexe Zahl. Unter der Krümmung versteht man den Kehrwert des Radius.

Ein Kreis mit unendlich großem Radius hat eine Krümmung von 0 und ist eine Gerade. Andersherum hat ein Kreis mit unendlicher Krümmung den Radius 0 -- ein Punkt.

Für unsere Anwendung ist es praktikabel, sich einen Kreis als durch eine Komplexe Zahl für den Mittelpunkt und einer reellen Zahl für die Krümmung festgelegt vorzustellen.

2.3   Des Satz von Descartes

Obwohl das Problem der sich berührenden Kreise erstmalig um 200 v.Chr. von Apollonius von Perge beschrieben und gelöst wurde, so ist eine Lösung von Hand jedoch eher umständlich.

Einen ersten Schritt in Richtung einer rechnerischen Lösung kommen wir mit dem von René Descartes im Jahr 1643 aufgestellten und nach ihm benannten Satz. Dieser besagt, dass für die Krümmungen \(b_1,b_2,b_3,b_4\) von vier sich paarweise berührenden Kreisen gilt:

\begin{equation} \label{eq:descartes} b_1^2 + b_2^2 + b_3^2 + b_4^2 = \frac{1}{2}(b_1+b_2+b_3+b_4)^2 \end{equation}

Dazu muss ein einschließender Kreis eine negative Krümmung aufweisen.

Kennen wir nun schon die ersten drei Krümmungen und stellen nach \(b_4\) um, so erhalten wir

\begin{equation} \label{eq:descartesb4} b_4 = b_1+b_2+b_3\pm 2\sqrt{b_1b_2 + b_2b_3 + b_3b_1} \end{equation}

Damit können wir aber zunächst nur die Krümmungen berechnen. Glücklicherweise haben Wilks et al im Jahr 2001 eine sehr ähnliche Gleichung für die Mittelpunkte von apollonischen Kreisen gefunden.

Für vier paarweise sich berührende Kreise mit den Krümmungen \(b_1, b_2, b_3, b_4\) und den Mittelpunkten \(z_1, z_2, z_3, z_4\) (als komplexen Zahlen) gilt:

\begin{equation} \label{eq:wilks} (b_1z_1)^2 + (b_2z_2)^2 + (b_3z_3)^2 + (b_4z_4)^2 = \frac{1}{2} (b_1z_1 + b_2z_2 + b_3z_3 + b_4z_4)^2 \end{equation}

Analog zum Satz von Descartes können wir wieder umstellen und damit \(z_4\) berechnen.

\begin{equation} \label{eq:wilksb4z4} b_4z_4 = b_1z_1+b_2z_2+b_3z_3\pm 2\sqrt{b_1b_2z_1z_2 + b_2b_3z_2z_3 + b_3b_1z_3z_1} \end{equation}

Wirklich elegant geht anders, aber diese Gleichungen werden auch nur ein einziges mal benötigt, nämlich wenn zu Beginn der einschließende Kreis für die drei Startkreise berechnet werden soll.

2.4   Vieta rettet den Tag

Mit unseren nun bekannten vier Kreisen können wir anfangen, den Satz von Vieta (\(\ref{eq:vietaplus}\)) anzuwenden.

Zuerst wieder für die Krümmungen:

Angenommen, wir kennen schon \(b_i, i=1\dots4\) und wollen die zweite Lösung für \(b_1\) ausrechnen, d.h. die Krümmung des anderen Kreises, der die Kreise 2, 3 und 4 berührt. Dazu stellen wir zunächst den Satz von Descartes ein wenig um.

\begin{eqnarray*} 2(b_1^2 + b_2^2 + b_3^2 + b_4^2) &=& (b_1+b_2+b_3+b_4)^2 \\ &=& b_1^2+b_2^2+b_3^2+b_4^2 +2b_1b_2 +2b_1b_3 +2b_1b_4 +2b_2b_3+2b_2b_4 +2b_3b_4 \\ &=& (b_1^2 + b_2^2 + b_3^2 + b_4^2) + 2b_1(b_2+b_3+b_4) +2(b_2b_3+b_2b_4 +b_3b_4) \end{eqnarray*}

Jetzt schaffen wir einiges Zeug auf die linke Seite

\begin{eqnarray*} b_1^2 + b_2^2 + b_3^2 + b_4^2 -2b_1(b_2+b_3+b_4)-2(b_2b_3+b_2b_4+b3_b4)&=& 0 \\ b_1^2 +\underbrace{(-2(b_2+b_3+b_4))}_{p}b_1 + \underbrace{(-2(b_2b_3+b_2b_4+b3_b4) + b_2^2 + b_3^2 + b_4^2)}_q&=&0 \end{eqnarray*}

Und können die Koeffizienten für den Satz von Vieta ablesen!

Eingesetzt in (\(\ref{eq:vietaplus}\)) erhalten wir \(b_{1,alt} + b_{1,neu}=-p1\).

\begin{equation} \label{eq:descartesvieta} b_{1,neu}=-p - b_{1,alt} = 2(b_2 + b_3 + b_4) - b_{1,alt} \end{equation}

Da wären wir. Mit dieser Formel können wir sehr elegant die Krümmung des zweiten Kreises aus den vier bekannten ausrechnen.

Für die zweite Gleichung und die Mittelpunkte der Kreise funktioniert es analog, es muss nur jeweils \(b_i\) durch \(b_iz_i\) ersetzt werden. Am Ende bekommt man

\begin{equation} \label{eq:wilksvieta} z_{1,neu}=\frac{2(b_2z_2 + b_3z_3 + b_4z_4) - b_{1,alt}z_{1,alt}}{b_{1,neu}} \end{equation}

2.5   Auffinden von Startkreisen

Jetzt, da wir wissen, wie wir aus bekannten Kreisen neue ausrechen können, brauchen wir nur noch eine Methode um an die ersten drei Kreise heranzukommen.

Sonderlich schwierig ist auch das nicht. Am Ende wollen wir ein Verfahren haben, dass uns zu drei beliebigen Krümmungen \(b_1, b_2, b_3\) die Mittelpunkte zugehöriger apollonischer Kreise und den einschließenden Kreis liefert.

  1. Für die Konstruktion brauchen wir die Radien und nicht die Krümmungen der Kreise. Setze also \(r_i=1/b_i\).

  2. Der erste Kreis ist einfach, wir packen ihn einfach in den Ursprung \((0,0)\)

    Erster Kreis im Ursprung
  3. Der zweite Kreis ist auch nicht viel schwieriger, wir setzen ihn rechts neben den ersten nach \((r_1+r_2,0)\).

    Zweiter Kreis rechts vom ersten.
  4. Der dritte ist schon etwas schwieriger. Wie kommen wir an die Koordinaten \((m_x,m_y)\)?

    Dritter Kreis unter den ersten beiden.

    Schauen wir uns einmal nur das Dreieck an, dass durch die drei Radien definiert wird.

    Dreieck aus den Radien.

    Man kann sagen, dass \((m_x,m_y)=(x,y)\) und weiterhin nach dem Satz des Pythagoras

    \begin{eqnarray*} y^2 + x^2&=& (r_1+r_3)^2 \\ y^2 + (r_1 + r_2 -x)^2 &=& (r_2+r_3)^2 \end{eqnarray*}

    Nach \(x\) und \(y\) aufgelöst erhalten wir

    \begin{eqnarray} \label{eq:thirdx} x = \frac{r_1^2 + r_1r_3 + r_1r_2 - r_2r_3}{r_1+r_2} \\ \label{eq:thirdy} y = \sqrt{(r_1+r_3)^2-x^2} \end{eqnarray}

  5. Der vierte Kreis ist wieder leicht, wir füttern einfach unsere Gleichungen (\(\ref{eq:descartesb4}\)) und (\(\ref{eq:wilksb4z4}\)) mit den drei schon festgelegten Kreisen.

2.6   Wann es nicht funktioniert

Es ist leider nicht möglich, aus allen Kombinationen von Startkreisen ein Apollonisches Netz zu erzeugen. Dies ist zum Beispiel der Fall, wenn es keinen einschließenden vierten Kreis gibt, da dieser einen unendlichen Radius haben müsste.

Eine Kombination, bei der das auftritt ist \(b_1=1,b_2=1,b_3=4\), oder allgemein jede Kombination, die \(b_3=2\sqrt{b_1\cdot b_2}+b_1+b_2\) erfüllt. Hier erhalten wir in Gleichung (\(\ref{eq:descartesb4}\)) den Wert \(b_4=0\).

3   Der Code

Der Quelltext ist auf github zu finden.

Zunächst brauchen wir eine Möglichkeit, die Informationen über einen Kreis (Position und Radius) zu speichern. Dies besorgt die Klasse Circle.

 1 class Circle(object):
 2     """
 3     A circle represented by center point as complex number and radius.
 4     """
 5     def __init__ ( self, mx, my, r ):
 6         """
 7         @param mx: x center coordinate
 8         @param my: y center coordinate
 9         @param r: radius
10         """
11         self.r = r
12         self.m = (mx +my*1j)
13 
14     def curvature (self):
15         return 1/self.r

Als nächstes brauchen wir den einschließenden Kreis für die drei Startkreise, wir erreichen das durch eine Implementierung der Formeln (\(\ref{eq:descartesb4}\)) und (\(\ref{eq:wilksb4z4}\)).

 1 def outerTangentCircle( circle1, circle2, circle3 ):
 2     """
 3     Takes three externally tangent circles and calculates the fourth one enclosing them.
 4     """
 5     cur1 = circle1.curvature()
 6     cur2 = circle2.curvature()
 7     cur3 = circle3.curvature()
 8     m1 = circle1.m
 9     m2 = circle2.m
10     m3 = circle3.m
11     cur4 = -2 * sqrt( cur1*cur2 + cur2*cur3 + cur1 * cur3 ) + cur1 + cur2 + cur3
12     m4 = ( -2 * sqrt( cur1*m1*cur2*m2 + cur2*m2*cur3*m3 + cur1*m1*cur3*m3 ) + cur1*m1 + cur2*m2 + cur3*m3 ) /  cur4
13     circle4 = Circle( m4.real, m4.imag, 1/cur4 )
14 
15     return circle4

Diese Funktion können wir jetzt direkt nutzen, wenn wir die ersten vier Kreise aus drei gegebenen Radien berechnen wollen. (Die Nummerierung der Kreise im Code unterscheidet sich von der Nummerierung in der Erklärung weiter oben).

 1 def tangentCirclesFromRadii( r2, r3, r4 ):
 2     """
 3     Takes three radii and calculates the corresponding externally
 4     tangent circles as well as a fourth one enclosing them. The
 5     enclosing circle is the first one.
 6     """
 7 
 8     circle2 = Circle( 0, 0, r2 )
 9     circle3 = Circle( r2 + r3, 0, r3 )
10     m4x = (r2*r2 + r2*r4 + r2*r3 - r3*r4) / (r2 + r3)
11     m4y = sqrt( (r2 + r4) * (r2 + r4) - m4x*m4x )
12     circle4 = Circle( m4x, m4y, r4 )
13     circle1 = outerTangentCircle( circle2, circle3, circle4 )
14     return ( circle1, circle2, circle3, circle4 )

Wir nähern uns dem interessanteren Teil. Hier springen wir mit Hilfe von Vieta von der bekannten ersten Lösung zur zweiten.

 1 def secondSolution( fixed, c1, c2, c3 ):
 2     """
 3     If given four tangent circles, calculate the other one that is tangent
 4     to the last three.
 5 
 6     @param fixed: The fixed circle touches the other three, but not
 7     the one to be calculated.
 8 
 9     @param c1, c2, c3: Three circles to which the other tangent circle
10     is to be calculated.
11     """
12 
13     curf = fixed.curvature()
14     cur1 = c1.curvature()
15     cur2 = c2.curvature()
16     cur3 = c3.curvature()
17 
18     curn = 2 * (cur1 + cur2 + cur3) - curf
19     mn = (2 * (cur1*c1.m + cur2*c2.m + cur3*c3.m) - curf*fixed.m ) / curn
20     return Circle( mn.real, mn.imag, 1/curn )

Das Apollonische Netz ist in der Klasse ApollonianGasket verpackt.

1 class ApollonianGasket(object):
2     """
3     Container for an Apollonian Gasket.
4     """
5     def __init__(self, c1, c2, c3):
6         self.start = tangentCirclesFromRadii( 1/c1, 1/c2, 1/c3 )
7         self.genCircles = list(self.start)

Und zum guten Schluss die Methode, die rekursiv die neuen Kreise aus den Startkreisen berechnet.

 1 def recurse(self, circles, depth, maxDepth):
 2     """Recursively calculate the smaller circles of the AG up to the
 3     given depth. Note that for depth n we get 2*3^{n+1} circles.
 4 
 5     @param maxDepth: Maximal depth of the recursion.
 6     @type maxDepth: int
 7 
 8     @param circles: 4-Tuple of circles for which the second
 9     solutions are calculated
10     @type circles: (L{Circle}, L{Circle}, L{Circle}, L{Circle})
11 
12     @param depth: Current depth
13     @type depth: int
14     """
15     if( depth == maxDepth ):
16         return
17     (c1, c2, c3, c4) = circles
18     if( depth == 0 ):
19         # First recursive step, this is the only time we need to
20         # calculate 4 new circles.
21         del self.genCircles[4:]
22         cspecial = secondSolution( c1, c2, c3, c4 )
23         self.genCircles.append( cspecial )
24         self.recurse( (cspecial, c2, c3, c4), 1, maxDepth )
25 
26     cn2 = secondSolution( c2, c1, c3, c4 )
27     self.genCircles.append( cn2 )
28     cn3 = secondSolution( c3, c1, c2, c4 )
29     self.genCircles.append( cn3 )
30     cn4 = secondSolution( c4, c1, c2, c3 )
31     self.genCircles.append( cn4 )
32 
33     self.recurse( (cn2, c1, c3, c4), depth+1, maxDepth )
34     self.recurse( (cn3, c1, c2, c4), depth+1, maxDepth )
35     self.recurse( (cn4, c1, c2, c3), depth+1, maxDepth )

Für jeden Kreis konstruieren wir mit den Satz von Vieta die zweite Lösung und steigen, mit dem betreffenden Kreis durch die zweite Lösung ersetzt, eine Ebene in die Rekursion ab. Beachte, dass wir den neuen Kreis an erster Stelle in das circles Quadrupel stecken und für eben diesen ersten Kreis (außer im ersten Durchlauf) nicht in die Rekursion absteigen. Damit vermeiden wir, dass wir die zweite Lösung der zweiten Lösung ausrechen, die ja wieder der ursprüngliche Kreis wäre.

Der restliche Code verpackt das Ganze noch in ein Kommandozeilenprogramm und enthält ein paar Hilfsfunktionen für die Farben und die Ausgabe als svg.

4   Ausprobieren!

Und noch ein apollonisches Netz.

Abbildung 6: Dasselbe in grün.

Neugierig geworden?

Zusätzlich zur Version für die Kommandozeile habe ich eine interaktive Online-Version geschrieben. Diese ist unter http://lsandig.org/cgi-bin/apollon/index.cgi zu finden.

Tags:

This text by Ludger Sandig is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.