Linien

Anschauliche Mathematik: Die Schmetterlingskurve

Seit ich Ende der 1980er Jahre mit meinem damals hochmodernen Atari Mega ST erste Schritte mit einem graphikfähigen Personalcomputer unternommen hatte, habe ich die Schmetterlingskurve immer wieder als Test für die Graphikfähigkeit und Schnelligkeit von Programmiersprachen und Rechnern benutzt. Sie wird in Polarkoordinaten beschrieben und ihre Formel ist

$$ \rho=e^{\cos(\theta)}-2\cdot \cos(4\cdot \theta)+\sin(\tfrac{\theta}{12})^5 $$

oder in Python-Code:

r = exp(cos(theta)) - 2*cos(4*theta) + (sin(theta/12))**5

Die Gleichung ist hinreichend kompliziert um selbst in C geschriebene Routinen auf meinen damals unglaubliche 8 MegaBit schnellen Atari alt aussehen zu lassen. Rechenzeiten von 10 - 20 Minuten waren keine Seltenheit. Heute dagegen muß man den Rechner schon künstlich verlangsamen, damit man sieht, wie sich die Kurve aufbaut. Denn sonst erscheint sofort die fertige Kurve, um die sinnliche Erfahrung, wie diese entsteht, wird man betrogen. Daher habe ich sie in Processing.py innerhalb der draw()-Schleife zeichnen lassen, wobei die Schleifenvariable theta bei jedem Durchlauf um 0.02 erhöht wurde.

Der Code ist – dank Processing.py – wieder von erfrischender Einfachheit und Kürze:

def setup():
    global theta, xOld, yOld
    theta = xOld = yOld = 0.0
    size(600, 600)
    background(100, 100, 100)
    colorMode(HSB, 100)

def draw():
    global theta, xOld, yOld
    strokeWeight(2)
    stroke(theta, 100 - theta, 100)
    r = exp(cos(theta)) - 2*cos(4*theta) + (sin(theta/12))**5
    # aus Polarkoordinaten konvertieren
    x = r*cos(theta)
    y = r*sin(theta)
    # auf Fenstergröße skalieren
    xx = (x*60) + 300
    yy = (y*60) + 300
    if (theta == 0.0):
        point(xx, yy)
    else:
        line(xOld, yOld, xx, yy)
    xOld = xx
    yOld = yy
    theta += 0.02
    if (theta > 75.39):
        print("I did it, Babe!")
        noLoop()

In setup() ist eigentlich nur bemerkenswert, daß ich nach der Festlegung des grauen Hintergrunds (noch als RGB), den colorMode auf HSB geändert habe. Damit lassen sich nämlich recht einfach diverse Farbeffekte erzielen. Ich habe dabei den Hue-Wert in Abhängigkeit von theta gesetzt, die Sättigung auf 100 - theta und die Brightness konstant bei 100 belassen. Da theta nie größer als \(75,39\) wird, wird es also auch nie größer als 100 und damit sind diese Umrechnungen gefahrlos.

Damit erreicht man, daß zu Beginn, wo die Sättigung noch ziemlich voll ist, die Zeichnung mit einem satten rot beginnt, während im Laufe der Iteration die weiteren Farben immer blasser werden. Ich fand dies das ästhetisch anspruchvollste Ergebnis, aber um das selber nachvollziehen zu können, sollten Sie ruhig damit experimentieren, zum Beispiel mit stroke(theta, 100, 100) oder stroke(100-theta, theta, 100) oder was immer Sie wollen.

Sie bekommen so diesen wunderschönen Schmetterling auf den Monitor gezeichnet:

Um die Entstehung der Kurve zu verstehen, empfiehlt Stan Wagon¹, nacheinander folgende Formlen plotten zu lassen:

In Polarkoordinaten:

r = exp(cos(theta))  # ergibt eine Art Kreis
r = -2*cos(4*theta)  # ergibt eine Art Blume
r = exp(cos(theta)) - 2*cos(4*theta)  # ergibt einen sehr einfachen Schmetterling

Dann in kartesischen Koordinaten:

x = -2*cos(4*theta)
y = -sin(theta/12)**5

Und dann ruhig auch noch einmal (wieder in Polarkoordinaten):

r = exp(cos(theta)) - 2*cos(4*theta) - (sin(theta/12))**5

Sie sehen dann, daß es eigentlich unerheblich ist, ob Sie den Störungsgteil der Formel addieren oder subtrahieren: Der Schmetterling ist nahezu identisch, lediglich an der anderen Farbgebung erkennen Sie, daß es zwei verschiedene Formeln sind.

Die Schmetterlingskurve und ähnliche Kurven wurden von Temple Fay² an der Universität von Southern Mississpi entwickelt. Sie eignen sich vorzüglich zum Experimentieren. So weist Pickover³ darauf hin, daß die Kurve

r = exp(cos(theta)) - 2.1*cos(6*theta) + sin(theta/30)**7

eine bedeutend größere Wiederholungsperiode besitzt. Sie sollten sich auch das ruhig einmal ansehen. Interessante Vergleiche mit der Originalschmetterlingskurve können Sie auch ziehen, wenn Sie mit

r = exp(cos(2*theta)) - 1.5*cos(4*theta)

eine ganz simple Form des Schmetterlings zeichnen lassen. Denn die heutigen Rechner sind schließlich hinreichend schnell, daß Sie nicht mehr minuten- oder gar stundenlang auf ein Ergebnis warten müssen und zum anderen lädt die Möglichkeit des schnellen Skizzierens mit der Processing-IDE geradezu zu eigenen Experimenten ein.

Der Lorenz-Attraktor, eine weitere Ikone der Chaos-Theorie

Nachdem ich im letzten Abschnitt die Schmetterlingskurve mit Processing.py gezeichnet hatte, wollte ich nun darauf aufbauen und eine Ikone der Chaos-Forschung, den Lorenz-Attraktor damit zeichnen. Ich hatte das ja auch schon einmal mit R getan – dort finden Sie auch weitere Hintergrundinformationen zu diesem Attraktor –, aber mit R wurde nur das fertige Ergebnis visualisiert. Hier kommt es mir aber wieder darauf an, die Entstehung der Kurve verfolgen zu können und dafür ist, wie schon bei der Schmetterlingskurve, Processing gut geeignet:

Als einer der ersten hatte 1961 Edward N. Lorenz, ein Meteorologe am Massachusetts Institute of Technology (MIT), erkannt, daß Iteration Chaos erzeugt. Er benutzte dort einen Computer, um ein einfaches nichtlineares Gleichungssystem zu lösen, das ein simples Modell der Luftströmungen in der Erdatmosphäre simulieren sollte. Dazu benutzte er ein System von sieben Differentialgleichungen, das Barry Saltzman im gleichen Jahr aus den Navier-Stokes-Gleichungen ⁴ hergeleitet hatte. Für dieses System existierte keine analytische Lösung, also mußte es numerisch (d.h. wie damals und auch heute noch vielfach üblich in FORTRAN) gelöst werden. Lorenz hatte entdeckt, daß bei nichtperiodischen Lösungen der Gleichungen vier der sieben Variablen gegen Null strebten. Daher konnte er das System auf drei Gleichungen reduzieren:

$$ \begin{align} \frac{dx}{dt} & = -\sigma (y - z) \\ \frac{dy}{dt} & = (\rho - z)x - y \\ \frac{dz}{dt} & = xy - \gamma z \end{align} $$

Dabei sind \(\sigma = -10\), \(\rho = 40\) und \(\gamma = - \frac{8}{3}\). Die Parameter der Gleichung habe ich [Herm1994] entnommen, [Stew1993] gibt \(\rho = 28\) an, aber der Wert ändert nichts an dem Verhalten der Kurve und \(\rho = 40\) füllt das Fenster einfach besser aus. 😜

Processing.py besitzt im Gegensatz zu R oder NumPy kein Modul zur numerischen Lösung von Differentialgleichungen und so habe ich das einfache Eulersche Poligonzugverfahren zur numerischen Berechnung benutzt

    dx = -sigma*(x - y)*dt
    dy = (x*(r - z) - y)*dt
    dz = (x*y - b*z)*dt
    x += dx
    y += dy
    z += dz

und dabei konstant dt = 0.01 gesetzt. Das benötigt natürlich mehr Rechenkapazität, als sie Lorenz je zur Verfügung standen, aber trotz der größeren Genauigkeit ändert sich nichts am chaotischen Verhalten der Kurve. Für die Farbberechnugn habe ich dieses mal nur den Farbwert (Hue) bei jeder Iteration geändert, Sättigung (Saturation) und Helligkeit (Brightness) bleiben konstant auf dem höchsten Wert. Das ergibt kräftige Farben, die von Rot über Orange nach Gelb und dann nach Grün, Blau und Violett wandern. So kann man schön erkennen, daß die beiden »Flügel« des Attraktors immer wieder, aber für uns unvorhersehbar, durchlaufen werden.

Der Quellcode

Hier nun der vollständige Quellcode des Skripts. Er ist kurz und selbsterklärend und folgt weitestgehend dem Pascal-Programm aus [Herm1994], Seiten 80ff.

b = 8.0/3
r = 40.0
sigma = 10.0
dt = 0.01
x = y = 0.01
z = t = 0.0
xOld = zOld = 0.0
first = True

def setup():
    size(640, 480)
    background(0)
    colorMode(HSB, 100)

def draw():
    global x, y, z, t, xOld, zOld
    global first
    strokeWeight(1)
    stroke(t, 100, 100)
    dx = -sigma*(x - y)*dt
    dy = (x*(r - z) - y)*dt
    dz = (x*y - b*z)*dt
    x += dx
    y += dy
    z += dz
    # auf Fenstergröße skalieren
    xx = (x*8) + 320
    zz = 470 - (z*5.5)
    if first:
        point(xx, zz)
    else:
       line(xOld, zOld, xx, zz)
    xOld = xx
    zOld = zz
    first = False 
    t = t + dt
    if ( t >= 75.0):
       print("I did it, Babe!")
       noLoop()

Links

• Der Lorenz Attractor auf Wolfram MathWorld

Literatur

• [Herm1994] Dieter Hermann: Algorithmen für Chaos und Fraktale, Bonn (Addison-Wesley) 1994, S. 80ff.
• [Pief1991] Frank Piefke: Simulationen mit dem Personalcomputer, Heidelberg (Hüthig) 1991
• [Stew1993] Ian Stewart: Spielt Gott Roulette?, Frankfurt (Insel TB) 1993