Dr. O. Hoffmann
Das übt ja auch...
Olaf Hoffmann
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.
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.
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.
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.
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.
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.
Eine Verallgemeinerung des Lorenz-Attraktors ist der Chen-Attraktor:
Chen-Attraktor,
(Quelltext zum Chen-Attraktor).
Differentialgleichung:
dx/dt = -k (y - x)
dy/dt = n x + l y - x z
dz/dt = x y - m z
Um die Richtung (x, y, z) = (p, q, r) wird das Problem um den Winkel o (in Grad) gedreht.
s ist ein Skalierungsparameter.
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
Translation 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
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).
Als anderes Beispiel soll gezeigt werden, wie sich zwei Teilchen in einem harmonischen Potential bewegen. Ohne Wechselwirkung untereinander handelt es sich bei den Trajektorien offenkundig einfach um Ellipsen, also eher langweilig. Ist die Anfangsgeschwindigkeit anfangs geeignet gewählt, erhält man sogar bloß jeweils die Schwingung in einer Ebene. Durch eine Wechselwirkung der beiden Teilchen kommt es allerdings zu komplizierteren Trajektorien. Je nach Anfangsbedingungen können diese entstehenden Bewegungsmuster allerdings dennoch Symmetrien aufweisen.
Zwei Objekte in Potential (1),
(Zwei Objekte in Potential (1)).
Die Wechselwirkung zwischen den Teilchen ist bei diesem Beispiel gleichfalls harmonisch.
k steht für das Zentralpotential (positiv für ein attraktives Potential), l steht entsprechend für die Stärke der Wechselwirkung zwischen den beiden Teilchen.
Ferner können die Anfangsbedingungen angegeben werden, also (x1,y1) als Koordinate des ersten Teilchens, (v1,w1) dessen Geschwindigkeit, entsprechend für das zweite Teilen (x2,y2) und (v2,w2).
Die Dauer der Animation dur und diese Parameter sind als GET-Parameter verfügbar.
Zwei Objekte in Potential (2),
(Zwei Objekte in Potential (2)).
In diesem Falle ist die Wechselwirkung zwischen den Teilchen lokal begrenzt.
Dabei steht der Parameter w für den ungefähren Radius der Wechselwirkung.
Die Abweichung von der harmonischen Wechselwirkung wird hier erreicht durch einen zusätzlichen Faktor, welcher der sinc-Funktion bis zur ersten Nullstelle entspricht.
Die Parameter sind die gleichen wie beim vorherigen Beispiel.
Zwei Objekte in Potential (3),
(Zwei Objekte in Potential (3)).
Abweichend vom vorherigen Beispiel ist der Faktor hier beim Abstand der Teilchen r voneinander: 1-r/w.
Die Parameter sind die gleichen wie beim vorherigen Beispiel.
Zwei Objekte in Potential (4),
(Zwei Objekte in Potential (4)).
Abweichend vom vorherigen Beispiel wirken hier die Kräfte zwischen den Teilchen nicht mehr in Richtung des jeweils anderen Teilchens, stattdessen senkrecht dazu, was eine ziemlich exotische Wechselwirkung darstellt.
Die Parameter sind die gleichen wie beim vorherigen Beispiel.
Zwei Objekte in Potential (5),
(Zwei Objekte in Potential (5)).
Harmonisches Zentralpotential und Wechselwirkung zwischen den Objekten.
Die Wechselwirkung ähnelt dem Gravitationspotential, ist allerdings so abgewandelt, daß der Parameter w einem weiteren Abstand in einer zusätzlichen Dimension entspricht, also effektiv den Abstand mindestens auf w vergrößert, so daß insbesondere keine Singularität entsteht.
Parameter k steht für die Stärke des Zentralpotentials, l für die Wechselwirkung, also negatives l ergibt eine attraktive Wechselwirkung.
Ansonsten sind die Parameter die gleichen wie beim vorherigen Beispiel.
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:
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.
Anmerkung: Ohne Graphikkarte oder leistungsfähigen Prozessor werden folgende Beispiele allenfalls sehr zäh laufen. Entsprechend kann es etwas dauern, bis der verwendete Brauser die Graphik darstellt; die Animation in sinnvollen Schritten wird eventuell sodann nicht möglich sein. Gleichfalls kann es etwas dauern, bis die Animation oder Darstellung abgebrochen werden kann.
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