XHTML SVG CSS PHP

Dr. O. Hoffmann

SVG- und PHP-Labor

Das übt ja auch...
Olaf Hoffmann

Numerische Lösung von Differentialgleichungen und Bewegungsgleichungen

Die meisten Prozesse, denen unser Interesse bei einer Animation gelten wird, sind analytisch nicht oder nur mit großem Aufwand lösbar. So kommt ein Autor einer Animation oft entweder mit Numerik schneller voran oder aber auch auf gar keinem anderen Wege, darum möchte ich hier einen kurzen Einblick in die Numerik von Differentialgleichungen geben, der einem zumindest bei einfachen Problemen dieser Art schon deutlich weiterhelfen sollte.
Nachteil numerischer Lösungen gegenüber analytischen ist natürlich immer die prinzipiell begrenzte Genauigkeit. Bei Anfangswertproblemen kumulieren natürlich die Fehler jedes einzelnen Schrittes. Bei chaotischen Systemen ist es sogar so, daß das Ergebnis nach längerer Integration einer Differentialgleichung komplett vom Verfahren oder der Rechenschrittgröße abhängen können - das Ergebnis repräsentiert dann keine Lösung des Anfangswertproblems mehr, sondern nur noch ein mehr oder weniger typisches Beispiel für das jeweilige System.

Für die numerische Lösung von Differentialgleichungen, zu denen auch die Bewegungsgleichungen gehören, gibt es zahlreiche Verfahren, hier werde ich nur ein einfaches Einschrittverfahren nach Runge und Kutta verwenden beziehungsweise für die Bewegungsgleichungen eines nach Runge, Kutta und Nyström.
Wie die genau funktionieren, wie genau das ist und welche anderen Verfahren es gibt, ist der einschlägigen Fachliteratur zum Thema zu entnehmen.
Ausgegangen wird hier von klassischen mechanischen Systemen, die sich auf Massenpunkte reduzieren lassen oder die sich durch wenige verallgemeinerte Koordinaten beschreiben lassen.

Da die Numerik hier ausschließlich mit PHP erfolgt, liegt von vorne herein fest, wie lang eine Animation sein wird. Die Dauer der Animation wird in einzelne Teilschritte zerlegt, deren Wert dann berechnet wird. Das kann manchmal etwas lästig sein, wenn zu Beginn nicht bekannt ist, wie lang die Animation sinnvoller Weise sein soll. Soll eine Trajektorie der Bewegung dargestellt werden, so ergibt sich daraus ohnehin ganz von selbst eine endliche und bekannte Länge der Animation.
Bei einer diskreten Skriptanimation ohne Trajektorienanspruch könnte auch anders vorgegangen werden und nur der jeweils nächste Schritt berechnet werden, ohne eine bestimmte Gesamtdauer festzulegen. Alternativ können vermutlich auch die jeweiligen Attribute der Animation neu berechnet werden, während die Animation läuft und dann ausgetauscht werden.
Bei der Verwendung von PHP kann hingegen so vorgegangen werden, daß das Endergebnis der Animation dem Skript wieder als Anfangswertproblem vorgesetzt wird. Dies passiert entweder manuell durch den Betrachter der Animation oder automatisch durch Senden eines geeigneten 'headers' mit einer verzögerten automatischen Weiterleitung. Da dies aber nicht direkt im HTTP festgelegt ist, sondern nur gängige Praxis bei den meisten Darstellungsprogrammen, sollte der letztere Fall nur alternativ zu einer manuellen Weiterleitungsmöglichkeit genutzt werden. Auch kann da ja gleich per Fallunterscheidung Vorsorge getroffen werden, wenn das Darstellungsprogramm mit automatischer Weiterleitung nicht klar kommt.

Runge-Kutta: Erster Test

Um sich mit dem Werkzeug erst einmal etwas vertraut zu machen, ist es in der Regel vorteilhaft, erstmal ein Problem zu lösen, dessen analytische Lösung bekannt und einfach darstellbar ist. Damit wird dann das numerische Ergebnis mit dem analytischen verglichen, um Fehlern auf die Spur zu kommen.
Einfach ist etwa eine Differentialgleichung vom Typ:
dx/dt = -kx
Ein Anfangswertproblem wird daraus, wenn festgelegt wird, welchen Wert x zum Zeitpunkt 0 hat, also x(0). Allgemein als Lösung bekannt ist:
x(t) = x(t=0) exp(-k t)
Solche Prozesse begegnen einem häufig, etwa das Entladen eines Kondensators, Abkühlen von Tee oder Kaffee, Auslaufen eines Behälters mit Flüssigkeit, biologische, chemische oder physikalische Populationsdynamik, etwa radioaktiver Zerfall oder natürliche Lebensdauer von angeregten elektronischen Zuständen von Atomen oder Molekülen und vieles mehr.

Exponentieller Abfall als Test für Runge-Kutta,
(Quelltext zum Test):
Eine Exponentialfunktion als Lösung einer numerischen Berechnung einer Differentialgleichung dx/dt = -k x.
k wird zufällig gewählt.
Es handelt sich um ein Anfangswertproblem, welches mit dem Runge-Kutta-Verfahren bearbeitet wird. Die exakte Lösung ist rot unterlegt. Es werden verschiedene Schrittweiten gewählt, um direkt zu sehen, wie die Qualität der numerischen Lösung von der Schrittweite abhängt. Nach welchen Kriterien die Schrittweite gewählt werden sollte, um eine sinnvolle Lösung zu bekommen, ist etwas kniffliger, es gibt dazu Bewertungszahlen und Kontrollmöglichkeiten, die hier nicht eingebaut sind, aber der einschlägigen Literatur zu entnehmen sind. Die Schrittweite kann auch während der Rechnung variiert werden, allerdings sind dann keyTimes zu verwenden und die Animationsdauer dur ist bei der Rechnung als Summe aller Schrittweiten zu bestimmen, der jeweilige Wert für keyTimes ist offenbar die Summe bis zum aktuellen Schritt. Sinnvoll ist eine Schrittweitenänderung, wenn sich die durch die Differentialgleichung gegebene Änderung der Größe stark ändert. Beispielsweise kann man die gleiche Rechnung mit halber Schrittweite und doppelter Schrittzahl durchführen, um zu verifizieren, ob die gefundene Lösung stabil ist.
Mittels der Knöpfe rechts unten (oder der entsprechenden GET-Parameter) kann Einfluß drauf genommen werden, welche Schrittweite verwendet wird und ob automatisch nach Ablauf einer Animation eine neue dargestellt werden soll oder nicht. Die Dauer der Animation kann nur ueber den GET-Parameter dur beeinflußt werden.

Ziemlich unabhängig von der Schrittweite sind hier die eigentlich berechneten Werte allerdings immer sehr ähnlich, verwendet man bei der großen Schrittweite die Ableitung, um mit einer kubischen Kurve zu interpolieren und nicht nur mit geraden Linien, wird man kaum einen Unterschied zur kleinen Schrittweite sehen. Dies ist also eine sehr stabile Differentialgleichung. Bei chaotischen Systemen hat die Schrittweite viel mehr Einfluß auf das numerische Ergebnis.

Differentialgleichung erster Ordnung

Gucken wir uns jetzt zwei Beispiele an, die nicht mehr ganz so einfach sind, zudem enthalten die Gleichungen auch noch eine explizite Zeitabhängigkeit, was die Runge-Kutta-Formel geringfügig umfangreicher macht.

DGL erster Ordnung 2,
(Quelltext zur DGL erster Ordnung 2):
Eine Differentialgleichung vom Typ
dx/dt= t/(k(1+x*x/l/l))
wird gelöst. k und l werden zufällig gewählt.
Mittels der Knöpfe rechts unten (oder der entsprechenden GET-Parameter) kann Einfluß drauf genommen werden, welche Schrittweite verwendet wird und ob automatisch nach Ablauf einer Animation eine neue dargestellt werden soll oder nicht. Die Dauer der Animation kann nur ueber den GET-Parameter dur beeinflußt werden.

DGL erster Ordnung 3,
(Quelltext zur DGL erster Ordnung 3):
Eine Differentialgleichung vom Typ
dx/dt= m sin(l t) -k x
wird gelöst. k, l und m werden zufällig gewählt. Mittels der Knöpfe rechts unten (oder der entsprechenden GET-Parameter) kann Einfluß drauf genommen werden, welche Schrittweite verwendet wird und ob automatisch nach Ablauf einer Animation eine neue dargestellt werden soll oder nicht.
Die Dauer (dur) der Animation und k, l und m können nur über den jeweiligen GET-Parameter beeinflußt werden.

Differentialgleichung zweiter Ordnung, Bewegungsgleichung

Bei einer Bewegungsgleichung handelt es sich um eine Differentialgleichung zweiter Ordnung, die allerdings auch in zwei gekoppelte Gleichungen erster Ordnung transformiert werden können, die ebenfalls mit Runge-Kutta behandelt werden können. Technisch erweist es sich als etwas eleganter, einfach das Runge-Kutta-Nyström-Verfahren anzuwenden.
Wieder gucken wir erstmal etwas an, was wir kennen, um grobe Fehler auszuschließen. Dazu verwenden wir einen harmonischen Oszillator, dessen Bewegungsmuster uns bereits aus anderen Abschnitten dieses Projektes und allgemein unseres bisherigen Lebens bekannt sein sollte. Die Differentialgleichung ist:
d2x/dt2 = -k x
oder eben
dv/dt = -k x
dx/dt = v

Harmonischer Oszillator numerisch,
(Quelltext zum harmonischen Oszillator numerisch):
Die horizonale Auslenkung eines Kreises stellt die harmonische Oszillation dar, die vertikale Richtung ist die Zeitachse. Die hellblaue Kurve ist x(t) und sollte eine trigonometrische Funktion (sinus, cosinus) darstellen. Die Qualität der Kurve kann natürlich stark verbessert werden, wenn statt eines Linienzuges eine kubische Bézierkurve verwendet wird. Prinzipiell kann man mittels animateMotion auch das bewegte Objekt einer solchen Kurve folgen lassen und dabei die Qualität der Animation mit dem Modus calcMode spline verbessern. Geschwindigkeit und Beschleunigung sind ohnehin bekannt, insofern ist das nur wenig mehr Aufwand, der hier dem geneigten Leser aber zur Übung überlassen bleibt. Verwendet man wirklich eine solche Kurve als Pfad für animateMotion, ist auch garantiert, daß die Pfadlänge eine monotone Funktion ist, weswegen man in dem Falle die numerischen Probleme mit calcMode spline elegant umgehen kann.

Gucken wir uns nun kurz etwas Komplizierteres an:
Getriebener Van der Poolscher Oszillator,
(Quelltext zum getriebenen Van der Poolschen Oszillator).
Folgende Bewegungsgleichung wird numerisch gelöst:
dv/dt=k (1 - l x x) y - x + m cos(n t)
dx/dt = v
Die Dauer der Animation dur, k, l, m und n sind als GET-Parameter verfügbar.
Beim Van der Poolschen Oszillator wirkt Reibung bei großen Auslenkungen und eine Beschleunigung bei kleinen. Zusätzlich kommt hier noch ein periodischer Antrieb hinzu.

Pendel ohne und mit Reibung

Beim Faden- Stangen- oder mathematischem Pendel treten Abweichungen von der harmonischen Schwingung auf, weswegen es sich lohnt, dieses an sich einfache System numerisch anzugehen. Die Periodendauer eines solchen Pendels etwa kann nur über ein elliptisches Integral angegeben werden, was es bei einem konkreten Beispiel etwas schwierig macht, die Dauer der Animation so genau zu bestimmen, daß man die Periodizität allein mit Wiederholung in der Animation erreichen kann. Bei diesem Beispiel wurde die Dauer aber lediglich über die ersten Glieder einer Potenzreihenentwicklung abgeschätzt. Es läßt sich da mit etwas mehr Aufwand und geschickten Tricks selbst mit weniger Daten in der SVG-Datei sicher noch deutlich mehr herausholen (Symmetrieeigenschaften ausnutzen):
Pendel ohne Reibung,
(Quelltext zum Pendel ohne Reibung).

Reibung ist meist eine recht komplizierte Wechselwirkung eines Objektes mit seiner Umgebung, bei der kinetische Energie in Wärme umgewandelt wird. Reibungskräfte sind stets der Bewegung entgegengerichtet und werden je nach umgebendem Medium mit verschiedenen Ansätzen in die Bewegungsgleichung eingebracht. Stokessche Reibung ist der Geschwindigkeit des Objektes proportional und beschreibt Reibung besonders bei hohen Geschwindigkeiten oder in zähen Flüssigkeiten recht gut. Newtonsche Reibung ist proportional zum Quadrat der Geschwindigkeit und wird eher bei langsamen Bewegungen zur Beschreibung von Reibung verwendet.

Pendel mit Stokesscher Reibung,
(Quelltext zum Pendel mit Stokesscher Reibung).
Ein Pendel mit Stokesscher Reibung als Differentialgleichung wird numerisch gelöst.
dv/dt=-k sin(x) -l v
dx/dt = v
Die Dauer der Animation dur, k und l sind als GET-Parameter verfügbar.

Pendel mit Newtonscher Reibung,
(Quelltext zum Pendel mit Newtonscher Reibung).
Ein Pendel mit Newtonscher Reibung als Differentialgleichung wird numerisch gelöst.
dv/dt=-k sin(x) -l v |v|
dx/dt = v
Die Dauer der Animation dur, k und l sind als GET-Parameter verfügbar.

Wird zu einem Pendel mit Reibung wiederum Energie hinzugefügt, etwa in periodischer Form, so kommt es leicht zu chaotischen Bewegungen, wenn die Ausschläge groß genug werden. Bei kleinen Ausschlägen ist das Verhalten entgegen ähnlich wie beim getriebenenen harmonischen Oszillator.
Getriebener Pendel mit Reibung,
(Quelltext zum getriebenen Pendel mit Reibung).
Ein getriebener Pendel mit Reibung als Differentialgleichung wird numerisch gelöst.
dv/dt=-k sin(x) - l v abs(v) - m sin(n t)
dx/dt = v
Die Dauer der Animation dur, k, l, n und m sind als GET-Parameter verfügbar.

Kombination aus Pendel und harmonischen Oszillator mit Reibung,
(Quelltext zur Kombination aus Pendel und harmonischen Oszillator mit Reibung).
Folgende Bewegungsgleichung wird numerisch gelöst: dv/dt=- k sin(x - m) - l x - n v
dx/dt = v
Die Dauer der Animation dur, k, l, n und m sind als GET-Parameter verfügbar.

Differentialgleichung mit mehreren Veränderlichen, Lorenz-Attraktor etc

Oft sind die Systeme natürlich komplizierter (und chaotischer) als ein Pendel oder ein harmonischer Oszillator. Der erste Ansatz ist immer, ein kompliziertes System in einfachere Einzelsysteme zu zerlegen. Stellt sich etwa heraus, daß es eine veränderliche Größe gibt, für die sich deren Differentialgleichung vom Restsystem separieren läßt, die Gleichung enthält also keine der anderen veränderlichen Größen, so ist es in der Regel vorteilhaft, diese Gleichung auch getrennt vom Rest zu lösen, wenn möglich natürlich exakt. Was übrig bleibt, kann wieder mit dem Runge-Kutta-Verfahren gelöst werden, wobei man versuchen kann, das problemorientiert etwas zu optimieren. Die Einzelschritte sind dann natürlich jeweils für alle Veränderliche durchzuführen, bevor der nächte Schritt vollzogen wird.

Nehmen wir als Beispiel für ein nicht auf Einzelprobleme reduzierbares System eines, welches in der Chaostheorie als einfaches System mit drei Veränderlichen sehr bekannt geworden ist:
Lorenz-Attraktor,
(Quelltext zum Lorenz-Attraktor).
Mittels der Knöpfe rechts unten (oder der entsprechenden GET-Parameter) kann Einfluß drauf genommen werden, ob die Animation fortgesetzt, automatisch fortgesetzt oder neu gestartet werden soll. Mittels der GET-Parameter kann man im Bedarfsfalle auch die Startposition und die Parameter selbst festlegen. s ist ein Skalierungsparameter.
Differentialgleichung Lorenz-Attraktor:
dx/dt = k(y-x)
dy/dt = l x - y - x z
dz/dt = x y - m z

Da kann die z-Komponente natürlich nicht eingesehen werden, daher kann man das Problem geeignet drehen, bevor es dargestellt wird:
Lorenz-Attraktor gedreht,
(Quelltext zum Lorenz-Attraktor gedreht).
Um die Richtung (x, y, z) = (p, q, r) wird das Problem um den Winkel o (in Grad) gedreht. dur ist die Animationsdauer in Sekunden (optional anz die Anzahl der Animationspunkte, per Voreinstellung 30*dur).

Das gleiche mit animateMotion, gzippt und im Profil tiny
Lorenz-Attraktor mit animateMotion,
(Quelltext zum Lorenz-Attraktor mit animateMotion).
Mittels der GET-Parameter tr und kr (0 oder 1) kann im Bedarfsfalle auch noch die Trajektorien-Malanimation beziehungsweise die animateMotion des Kreises (de)aktiviert werden, Voreinstellung ist 1, also anzeigen.

Um zu zeigen, wie langsam oder schnell das System chaotisch wird, wird zum einen die Schrittweite etwas vergrößert, dann aber auch geringfügig verschiedene Schrittweiten miteinander verglichen:
Lorenz-Attraktoren mit unterschiedlichen Schrittweiten,
(Quelltext zu Lorenz-Attraktoren mit unterschiedlichen Schrittweiten).
Bei gleicher Anfangsbedingung haben der grüne, der blaue und der rote Punkt jeweils eine geringfügig andere Schrittweite. Anfangs zeigt sich noch eine gute Übereinstimmung. Irgendwann gehen die Kreise völlig unterschiedliche Richtungen.

Eine andere Methode, das chaotische am System sichtbar zu machen, ist eine kleine Variantion der Anfangsbedingungen:
Lorenz-Attraktoren mit leicht unterschiedlichen Anfangsbedingungen,
(Quelltext zu Lorenz-Attraktoren mit leicht unterschiedlichen Anfangsbedingungen).
Der grüne, der blaue und der rote Punkt haben jeweils geringfügig andere Anfangsbedingungen. Anfangs zeigt sich noch eine gute Übereinstimmung. Irgendwann gehen die Kreise völlig unterschiedliche Richtungen.

Durch wenige Modifikationen kann das natürlich auch für andere Systeme verwendet werden. Viele davon haben einen realen naturwissenschaftlichen Hintergrund, sind in der gemeinhin diskutierten Form aber immer stark abstrahiert. Unter den angegebenen Namen sollte es recht einfach weitere inhaltliche Informationen zu den Systemen geben, auf die ich hier nicht weiter eingehen werde.
Es kann von Vorteil sein, die Schrittweite anzupassen und eine Tanslation vor der Drehung anzuwenden, die den Mittelwert der Fixpunkte ungefähr auf den Drehnullpunkt verschiebt.

Rößler-Attraktor,
(Quelltext zum Rößler-Attraktor).
Differentialgleichung:
dx/dt = -y -z
dy/dt = x + k y
dz/dt = l + z (x- m)

Rikitake-Attraktor,
(Quelltext zum Rikitake-Attraktor).
Differentialgleichung:
dx/dt = -lx + zy
dy/dt = -ly + (z - k)x
dz/dt = 1 - xy

Lotka-Volterra-Attraktor,
(Quelltext zum Lotka-Volterra-Attraktor).
Differentialgleichung:
dx/dt = x - x y +m x x - k z x x
dy/dt = x y - y
dz/dt = k z x x - l z

Getriebener Chua-Attraktor,
(Quelltext zum getriebenen Chua-Attraktor).
Differentialgleichung:
dx/dt = k (y - m x - n x x x/3)
dy/dt = x - y + z
dz/dt = - l y sin(j t)

Rabinovich-Fabrikant-Attraktor,
(Quelltext zum Rabinovich-Fabrikant-Attraktor).
Differentialgleichung:
dx/dt = y (z - 1 + x x) + l x
dy/dt = x (3 z + 1 -x x) + l y
dz/dt = -2 z (k + x y)

Moore-Spiegel-Attraktor,
(Quelltext zum Moore-Spiegel-Attraktor).
Differentialgleichung:
dx/dt = y
dy/dt = z
dz/dt = - z - (l -k + k x x)y - l x

Lorenz84-Attraktor,
(Quelltext zum Lorenz84-Attraktor).
Differentialgleichung:
dx/dt = -k x - y y - z z + k m
dy/dt = - y + x y - l x z + n
dz/dt = - z + l x y + x z

Kleinere Modifikationen zum vorherigen:
AttrH01-Attraktor,
(Quelltext zum AttrH01-Attraktor).
Differentialgleichung:
dx/dt =-y y -z z -k x +k m
dy/dt = x y +z x l -x +n
dz/dt =- l y x + x z -z

Kleinere Modifikationen zum vorherigen:
AttrH02-Attraktor,
(Quelltext zum AttrH02-Attraktor).
Differentialgleichung:
dx/dt =-y y -z z -k y +n
dy/dt = x y +z x m -x +n
dz/dt =- l y x + x z -z

Bewegungsgleichung mit mehreren Veränderlichen

Fangen wir wieder vorsichtig an mit einer Differentialgleichung, deren Lösung wir schon kennen, um zu kontrollieren, ob alles richtig funktioniert. Dazu verwenden wir gekoppelte harmonsiche Oszillatoren, die Lösung ist bekannt und auch im Abschnitt Bewegung dargestellt:
Gekoppelte Oszillatoren,
(Quelltext zu gekoppelten Oszillatoren).
Differentialgleichung:
dv1/dt=- k x1 + l(x2 - x1)
dv2/dt=- m x2 + n(x1 - x2)
dx1/dt = v1
dx2/dt = v2
Die Dauer der Animation dur, k, l, m und n sind als GET-Parameter verfügbar.
Gekoppelte Oszillatoren, alternative Ansicht,
(Quelltext zu gekoppelte Oszillatoren, alternative Ansicht).

Bei einem Oszillator sieht ja noch alles schön regelmäßig aus, ganz anders kann das Bild bei gekoppelten Pendeln aussehen.
Gekoppelte Pendel,
(Quelltext zu gekoppelte Pendel).
Differentialgleichung:
dv1/dt=- k sin(π x1) + l(x2 - x1)
dv2/dt=- m sin(π x2) + n(x1 - x2)
dx1/dt = v1
dx2/dt = v2
Die Dauer der Animation dur, k, l, m und n sind als GET-Parameter verfügbar.
Gekoppelte Pendel, alternative Ansicht,
(Quelltext zu gekoppelte Pendel, alternative Ansicht).

Wie die Kopplung genau aussieht, hängt natürlich vom Problem ab oder wie die Kopplung praktisch realisiert ist.
Gekoppelte Pendel 2,
(Quelltext zu gekoppelte Pendel 2).
Differentialgleichung:
dv1/dt=- k sin(π x1) + l sin(π(x2 - x1))
dv2/dt=- m sin(π x2) + n sin(π(x1 - x2))
dx1/dt = v1
dx2/dt = v2
Die Dauer der Animation dur, k, l, m und n sind als GET-Parameter verfügbar.

Doppelpendel mit gleichen Massen und gleichen Pendellängen,
(Quelltext zum Doppelpendel).
Folgende Bewegungsgleichung wird numerisch gelöst:
Hilfsgrößen:
c=cos(x1-x2)
s=sin(x1-x2)
f=(2 - c c)/2
Differentialgleichung:
dv1/dt=-0.5/f (v2 v2 s +v1 v1 s c + 2 k sin(x1) -k c sin(x2))
dv2/dt=1/f (v2 v2 s c/2 +v1 v1 s +k c sin(x1) -k sin(x2))
dx1/dt = v1
dx2/dt = v2
Die xi sind die Winkel relativ zur Vertikalen, die vi die Winkelgeschwindigkeiten. Die Dauer der Animation dur, Anfangswerte und k sind als GET-Parameter verfügbar.

Doppelpendel mit unterschiedlichen Massen und Pendellängen,
(Quelltext zum Doppelpendel mit unterschiedlichen Massen und Pendellängen).
Folgende Bewegungsgleichung wird numerisch gelöst:
Hilfsgrößen:
c=cos(x2-x1)
s=sin(x2-x1)
f1=l1 - m2 l1 c c
f2=l2 - m2 l2 c c
Differentialgleichung:
dv1/dt=(m2 l2 v2 v2 s +m2 l1 v1 v1 s c -g sin(x1)+m2 g c sin(x2))/f1
dv2/dt=(-m2 l2 v2 v2 s c -l1 v1 v1 s +g c sin(x1)- g sin(x2))/f2
dx1/dt = v1
dx2/dt = v2
Die xi sind die Winkel relativ zur Vertikalen, die vi die Winkelgeschwindigkeiten. li sind die relativen Pendellängen, mi ihre relativen Massen, g die (Erd)beschleunigung (um die 9.8 m/s2); Die Summe der Pendellängen ist 1 und die der Massen ebenfalls. Die Dauer der Animation dur, Anfangswerte und k sind als GET-Parameter verfügbar.

Sofern GZIP nicht vom server alleine durchgeführt wird, kann die Dateigröße natürlich noch ordentlich reduziert werden, indem man das mit PHP macht:
Doppelpendel mit unterschiedlichen Massen und Pendellängen, Alternative 'GZIP't,
(Quelltext zum Doppelpendel mit unterschiedlichen Massen und Pendellängen, Alternative 'GZIP't).

Schwärme

Ein weitere kompliziertere Anwendung ist die Simulation von Schwärmen, wie sie bei einigen Vogelarten, Fledermäusen, Fischen, Herdentieren etc vorkommen, verwandte Probleme kann es aber auch in der Physik geben, wenn etwa Atome oder Ionen in einer Teilchenfalle ein (Bose-Einstein-)-Kondensat oder auch ein klassisches Kondensat bilden - oder gar bei niedrigen kinetischen Energien (Temperaturen) gar zwei- oder dreidimensionale Kristallstrukturen ausbilden. Ferner kann man auch vermuten, daß die Entstehung von Planetensystemen, Galaxien, Superhaufen etc auch teilweise mit solchen Phänomenen zusammenhängt.

Dabei geht es immer um viele (als mehr als zwei) Objekte, die miteinander wechselwirken. Zumeist wechselwirkt jedes Objekt mit jedem anderen, was die Rechnung aufwendig macht, der Aufwand steigt zum einen ungefähr quadratisch mit der Anzahl der Objekte und linear mit der Anzahl der gewünschten Rechenschritte. Die Art der Wechselwirkungen ist zumindest bei der Simulation von Tierschwärmen eher heuristisch angenommen, von daher ist es hier nicht notwendig, die Numerik besonders genau zu rechnen, ein einfaches Einschritt-Verfahren (nach Euler, also jenes, was einem als erstes auch selber ohne weitere Fehlerrechnung einfallen wird) wird hier zumeist reichen, da es nur darum geht, ein typisches Verhalten zu simulieren und nicht genaue Lösungen eines Anfangswertproblems zu finden. Man sollte nur darauf achten, daß die Schritte nicht so groß werden, daß durch die Ungenauigkeit das Verhalten komplett verfälscht wird.

Der allgemeine Befund der Schwarmforschung führt jedenfalls dazu, daß folgende Komponenten erforderlich sind, damit sich Schwärme bilden oder sich ein spezifisches Schwarmverhalten ausbildet:

  1. Kleine Abstände zwischen zwei Objekten des Schwarms haben einen repulsiven Effekt, um das Risiko von Zusammenstößen zu minimieren
  2. Bei größeren Abständen gibt es einen attraktiven Effekt, welcher es begünstigt, daß sich Gruppen zusammenschließen
  3. Ein Objekt versucht seine Geschwindigkeit in Richtung und Betrag an die von nahen Nachbarn anzupassen

Bei unbelebten Objekten wie Molekülen oder Atomen kann man die ersten beiden Regeln auf das Molekülpotential zurückführen, welches meist bereits eine solche Form hat. Es tragen also alle Objekte zum Effekt bei, weit entfernte allerdings nur sehr wenig. Die dritte Regel entspricht einer Reibung, in dem Falle etwa hervorgerufen durch andere Teilchen oder Felder (Licht).

Bei Tieren kann es hingegen plausibel sein anzunehmen, daß nur die nächsten Nachbarn berücksichtigt werden und vielleicht auch nur solche, die grob in eine ähnliche Richtung fliegen wie das gerade betrachtete Tier selbst, während Zusammenstöße nach Möglichkeit immer vermieden werden. Ein Tier hat zudem eine maximale Geschwindigkeit und kann sich meist besser nach vorne als nach hinten orientieren. Bei fliegenden Schwärmen muß man zudem berücksichtigen, daß die Tiere eine minimale Geschwindigkeit nicht unterschreiten dürfen, um nicht abzustürzen. Tiere können zudem eine eigene Motivation oder Zielvorgabe haben, etwa bei Zugvögeln, die ein bestimmtes Ziel erreichen wollen oder auch wenn Freßfeinde ins Sichtfeld geraten, die dann zumindest einen stark repulsiven Effekt haben werden, bei intelligenteren Tieren aber natürlich auch dazu führen können, gezielt andere Mitglieder des Schwarms zwischen sich und den Freßfeind zu bringen, um das Risiko, selbst erwischt zu werden, zu minimieren.

Je nach Problem kann es also notwendig sein, mehr oder weniger spezifische Modifikationen am allgemeinen Modell vorzunehmen. Aber auch relativ simple Annahmen führen bereits zu Schwärmen, deren Verhalten man mit diversen Parametern leicht beeinflussen kann, welche dann die Effekte relativ zueinander gewichten. Die Annahme etwa einer etwas anderen Wechselwirkung führt dabei nicht notwendig zu einem komplett anderen Verhalten. Zusätzliche Regeln oder aber Verfeinerungen wie die genannte, daß die Geschwindigkeit nur an Objekte angepaßt wird, die sich in eine ähnliche Richtung fortbewegen wie das betrachtete Objekt selbst oder auch die Beschränkung auf eine maximale und minimale Geschwindigkeit können natürlich zu einem deutlich anderen Ergebnis führen.

Bei den im Folgenden betrachteten Beispielen wird aus praktischen Gründen meist ein harmonisches Zentralpotential verwendet, um die Objekte in einem endlichen Bereich zu halten. Solche Umweltbedingungen haben natürlich auch Einfluß auf das Schwarmverhalten, zwangsläufig bleiben die simulierten Individuen dabei enger zusammen als ohne ein solches Potential, was grob simuliert, daß die Tiere in einem endlichen Gebiet eingesperrt sind, etwa bei Fischen ein See oder bei Herden eine Insel, bei Vögeln eine Voliere.

Um die Bildung eines Schwarm zu beschleunigen und damit die Rechenzeit so zu begrenzen, daß mit den üblichen Einstellungen von PHP ein Effekt zu sehen ist, wird eine Wechselwirkung angenommen, bei der der attraktive Teil umgekehrt proportional zum Abstand ist und der repulsive umgekehrt proportional zum Quadrat des Abstandes ist. Das ist auch bei Tieren insofern nicht unplausibel, weil diese auch weit entfernte Objekte sehen können.

Schwarm (1) (Quelltext zum Schwarm (1))
Ein Schwarm aus bunten Teilchen bildet sich in einem harmonischen Potential, wobei sich die Teilchen anfangs ungefähr auf einer Kreisbahn um das Zentrum des Potentials bewegen. Die Wechselwirkung zwischen zwei Teilchen ist bei kleinem Abstand repulsiv, bei großem Abstand etwas attraktiv. Zusätzlich gibt es eine kleine Geschwindigkeitskorrektur für jedes Teilchen bei jedem Rechenschritt, die aus dem gewichteten Mittel der Geschwindigkeiten der anderen Teilchen besteht. Bei großem Abstand ist dabei der Korrekturbeitrag klein. Insgesamt wirkt diese Korrektur ähnlich wie Reibung, entzieht dem System also kinetische Energie. Zudem kann eine Zufallsgeschwindigkeit bei jedem Schritt hinzuaddiert werden, dies entspricht einem Wärmebad, welches ungerichtet kinetische Energie hinzufügt.
Je nach Gewichtung kommt es meist früher oder später zu einer Schwarmbildung oder sogar zu einer Kondensation zu einer flüssigkeitsähnlichen Anordnung oder sogar zu einer kristallinen Anordnung.
Geometrische Parameter:
x: horizontale Ausdehnung des Bildes
y: vertikale Ausdehnung des Bildes
m: Maßeinheit dafür, %, px, cm, mm
Schwarmparameter:
n: Anzahl der Durchläufe und Animationsdauer in zehntel Sekunden, n maximal 10000
anz: Anzahl der Teilchen; um die Rechenzeit und die Dateigröße zu begrenzen, ist festgelegt, daß n * anz² maximal 1e12
Gewichte (typisch um 1):
pot: Zentralpotential
ww: Teilchen-Teilchen-Potential
av: Geschwindigkeitsanpassung an Nachbarn
zv: Zufällige Geschwindigkeitsänderung
sz: Szenario: 1: Rotation; 2: Beschuß; sonst: Zufall

Schwarm (2) (Quelltext zum Schwarm (2))
Ähnlich wie der erste Schwarm, aber es gibt nur das erste Szenario. Ansonsten ist der Parameter für das Zentralpotential kleiner als 0.1, gibt es eine Bewegung in eine zufällige Richtung, die Anzeige folgt der mittleren Geschwindigkeit der Teilchen) Zudem ist die Begrenzung des repulsiven Teils zum Abstand 0 etwas eleganter gelöst, das Zentralpotential wird dadurch etwas verändert.

Schwarm (3) (Quelltext zum Schwarm (3))
Ähnlich wie der zweite Schwarm. Hinsichtlich der Gewichtung der Geschwindigkeitsanpassung ist allerdings sowohl der Abstand als auch die Richtung entscheidend, es gibt keine Korrektur bei gegenläufigen Geschwindigkeiten, je ähnlicher die Richtungen, desto stärker die Angleichung.
Weitere Parameter:
sv: Geschwindigkeitsverteilung anfangs stark zufällig? 0 ohne Zufall
vmin: minmale Geschwindigkeit ( größer oder gleich 0)
vmax: maximale Geschwindigkeit (größer oder gleich vmin)

Schwarm (4) (Quelltext zum Schwarm (4))
Ähnlich wie der dritte Schwarm, allerdings wirkt hier die richtungsabhängige Anpassung nicht nur auf die Geschwindigkeit, sondern auch auf das Teilchen-Teilchen-Potential. Dies ergibt etwa Sinn bei Tieren, die sich nur nach vorne oder zu Seite orientieren, nicht aber nach hinten.

Derartige richtungsabhängige Korrekturen führen eher dazu, daß sich ein Schwarm langsamer ausbildet, ist aber wohl realistischer für Lebewesen, bei denen die Wechselwirkung primär darauf basiert, daß diese einander wahrnehmen und diese Wahrnehmung richtungsabhängig ist, wie das etwa für das Sehen der Fall ist, weniger etwa fürs Hören.

Bei folgendem Beispiel wird statt auf eine maximale und eine minimale Geschwindigkeit zu setzen eine ungefähre Sollgeschwindigkeit vorgegeben. Zudem bestimmt ein Sichtwinkel, welche anderen Teilchen Einfluß auf die Flugrichtung und die Anziehung haben. Auch hier gibt es aber immer eine stärkere Gewichtung für Teilchen ungefähr in Vorwärtsrichtung.

Schwarm (5) (Quelltext zum Schwarm (5))
Weitere Parameter:
vsoll: ungefähre Sollgeschwindigkeit (~200)
vv: Variation der Sollgeschwindigkeit (0 bis 1)
sw: Sichtwinkel in Grad (sofern nicht angegeben nach vorne, 180)

Als Variation des Themas sei nun noch eine Simulation eines Schwarms von Leuchtkäfern betrachtet. Die bewegen sich nun in der Simulation nicht, sondern sitzen in regelmäßiger Struktur nebeneinander und reagieren darauf, wenn ihre Nachbarn leuchten. Auch hier führen einfache Regeln dazu, daß es zu Schwarmverhalten kommen kann, wobei die Käfer nur endlich lange und mit Pause leuchten können und nur auf die nächsten Nachbar reagieren. Statt Differentialgleichungen zu lösen, wird hier im Wesentlichen nur noch geguckt, wieviele Nachbarn leuchten. Je nach Parameterkombination kann ein ähnlicher Eindruck entstehen wie bei Conways Spiel des Lebens.

Leuchtschwarm (Quelltext zum Leuchtschwarm)
In hexagonaler Anordnung sind Punkte angeordnet, die entweder hell oder dunkel sind. Die Anordnung erfolgt jeweils mit periodischen Randbedingungen. Dies simuliert grob eine Kolonie von Leuchtkäfern, von denen jeder nur endlich lang leuchten kann und danach eine Pause braucht. Ein Käfer entscheidet entweder spontan, daß er leuchten will, oder er richtet sich nach seinen nächsten Nachbarn. Beides nur, wenn er gerade nicht leuchtet oder keine Pause nach dem Leuchten einlegen muß.
Durch die Pausen und die Verzögerung hängt es nun stark von diesen Parametern ab und davon, wieviele Nachbarn leuchten müssen, damit ein Käfer selbst anfängt zu leuchten, wie sich die Situation entwickelt, die mit einer zufälligen Verteilung von leuchtenden oder nicht leuchtenden Käfern beginnt.
Geometrische Parameter:
x: horizontale Ausdehnung des Bildes
y: vertikale Ausdehnung des Bildes
m: Maßeinheit dafür, %, px, cm, mm
Schwarmparameter:
n: Anzahl der Durchläufe und Animationsdauer in Sekunden, n maximal 10000
abs: Abstand zu den 6 nächsten Nachbarn, 200-2000
ld: Leuchtdauer in Sekunden, 0-100
pd: Minimale Pause nach dem Leuchten, 0-100
wl: Wahrscheinlichkeit spontanen Leuchtens, falls gerade nicht geleuchtet wird, in Promille-Werten, maximal zwei Nachkommastellen, 0-1000
nb: Wieviele Nachbarn müssen mindestens leuchten, um eigenes Leuchten zu induzieren, 1-6
wa: Welcher Anteil ist anfangs an, in Promille-Werten, maximal zwei Nachkommastellen, 0-1000
kombi: Interessante Kombinationen, überschreibt Schwarmparameter nb, ld, pd, wl und wa, 1-14

Ähnlich wie zuvor, nur die Abhängigkeit von den Nachbarn wird mit dem reziproken des Abstandsquadrates gewichtet, für die 6 nächsten Nachbarn gibt es jeweils 12 Gewichtungspunkte, für die nächsten 6 gibt es 4 Gewichtungspunkte und für die nächsten 6 noch 3 Gewichtungspunkte.
Entsprechend ist die Zuordnung des Parameters kombi etwas anders: nb.?, wobei ? für eine Ziffer steht, es gibt für nb 3,4, 6-60 mindestens einen Eintrag für ? = 1:
Leuchtschwarm (2) (Quelltext zum Leuchtschwarm (2))

Die Gewichtung ändert nichts Wesentliches am prinzipiellen Verhalten. Bei der Notwendigkeit von nur wenigen Nachbarn ergibt sich meist ein starkes Wachstum der leuchtenden Bereiche, welches aber auch bedingt durch die Pause nach dem Leuchten im Anschluß zu dunkel bleibenden Bereichen führt. Bei der Notwendigkeit von vielen Nachbarn gibt es allenfalls eine interessante Dynamik für die dunklen Bereiche, bei nicht hinreichend guten Anfangsbedingungen aber ein mehr oder weniger schnelles Einstellen der Aktivitäten des Schwarm. Spontanes Leuchten kann insbesondere hier dem dauerhaften Dunkel entgegenwirken. Lange Pausen führen meist auch zu einem dauerhaften Dunkel. Nur bei bestimmten Parameterkombinationen ergibt sich dabei eine Expansion, die auch im Ursprungsbereich der Dynamik wieder ein Leuchten anregt.

Wie nicht anders zu erwarten tritt die spannende Dynamik irgendwo im mittleren Bereich auf. Dort kann es Überstrukturen geben, die sich periodisch ändern, sich dabei fortbewegen oder aber auch Strukturen, die immer wieder neue Wellen von Aktivitäten erzeugen, die expandieren.

Als Beispiel für umgedrehte Gewichtung der 24 nächsten Nachbarn:
Leuchtschwarm (3) (Quelltext zum Leuchtschwarm (3))
Die 6 nächsten Nachbarn bekommen hier je nur einen Gewichtungspunkt, die nächsten 6 je 3 und die nächsten 6 dann 4 Gewichtungspunkte. Ein paar interessante Kombinationen gibt es wieder mit dem Parameter kombi, hier von 1 bis 5, gibt natürlich viel mehr interessante Strukturen.

Weil das schon erwähnt wurde, noch eine hexagonale Variation von Conways Spiel des Lebens:
Spiel des Lebens (Quelltext zum Spiel des Lebens)
In hexagonaler Anordnung sind Punkte angeordnet, die entweder hell oder dunkel sind. Die Anordnung erfolgt jeweils mit periodischen Randbedingungen.
Analog zu Conways Spiel des Lebens richtet sich hell oder dunkel eines Durchlaufes danach, wieviele und welche Nachbarn im vorherigen Durchlauf hell waren. Zu viele oder zu wenige bewirken dunkel, bei einer Zahl dazwischen ist es hell.
Geometrische Parameter:
x: horizontale Ausdehnung des Bildes
y: vertikale Ausdehnung des Bildes
m: Maßeinheit dafür, %, px, cm, mm
Schwarmparameter: n: Anzahl der Durchläufe und Animationsdauer in drittel Sekunden, n maximal 10000
abs: Abstand zu den 6 nächsten Nachbarn, 200-2000
Kriterium für hell oder dunkel bestimmt sich mit den Parametern nt, ne und typ.
typ bestimmt, welche Nachbarn wie berücksichtigt werden:
1 - die sechs nächsten Nachbarn, jeweils Gewicht 1
2 - die achtzehn nächsten Nachbarn mit gleichem Gewicht 1
3 - die achtzehn nächsten Nachbarn, Gewicht 12 jeweils für die 6 nächsten Nachbarn, 4 jeweils für die 6 übernächsten Nachbarn, 3 jeweils für die 6 Nachbarn, die darauf folgen.
4 - Gewicht 2 für die nächsten 6 Nachbarn, Gewicht 1 für die übernächsten 12
Ist die Summe der Gewichte der Nachbarn gleich nt oder größer, so wird es dunkel.
Ist die Summe der Gewichte der Nachbarn kleiner als ne, so wird es ebenfalls dunkel. Dazwischen wird es hell.
wl: Wahrscheinlichkeit spontanen Leuchtens, falls gerade nicht geleuchtet wird, in Promille-Werten, maximal zwei Nachkommastellen, 0-1000
wa: Welcher Anteil ist anfangs an, in Promille-Werten, maximal zwei Nachkommastellen, 0-1000
kombi: Interessante Kombinationen, überschreibt Schwarmparameter ne, nt, typ, wl und wa, 1-21