beltoforion.de
 Impressum |  CSS |  XHTML

Das magnetische Pendel

Demonstration des Schmetterlingseffektes mit einem Magnetpendel

Einleitung


Vor einigen Jahren gab mir mein Physikprofessor an der UNI einen Artikel aus der Zeitschrift Spektrum der Wissenschaft, der "Experimente zum Chaos" hieß. Der Artikel stammt aus dem Jahr 1994, und stellt unter anderem eine Simulation vor, welche die chaotische Bewegung eines Pendels unter dem Einfluß der Schwerkraft und dreier Magnete simuliert. Diese Simulation fand ich faszinierend und implementierte sie daher selbst in einem kleinen Programm. Seitdem geriet es bei mir fast in Vergessenheit.
Kürzlich bin ich wieder darüber gestolpert und beschloß das Programm neu zu schreiben und hier zu veröffentlichen. Hier ist es nun also...

Es sollte vielleicht noch erwähnt werden, daß eine Berechnung, je nach vorhandener Computerhardware und Bildgröße einige Zeit dauern kann. Zeiten von 4-5 Stunden sind nicht ungewöhnlich für ein Bild mit einer Auflösung von 2000 x 2000 Pixeln. Zunächst ein kleiner Überblick, was euch im weiteren Verlauf dieses Artikels erwartet.

Was ihr mitbringen solltet: Was ihr bekommt:

Der Schmetterlingseffekt


Das hier vorgestellte Programm demonstriert den Schmetterlingseffekt. Für diejenigen, die mit dem Begriff nicht vertraut empfielt sich der entsprechende Eintrag in der Wikipedia. Auszug:

Als Schmetterlingseffekt (engl. butterfly effect) bezeichnet man den Effekt, dass in komplexen, dynamischen Systemen eine große Empfindlichkeit auf kleine Abweichungen in den Anfangsbedingungen besteht. Geringfügig veränderte Anfangsbedingungen können im langfristigen Verlauf zu einer völlig anderen Entwicklung führen. Die Bezeichnung Schmetterlingseffekt stammt von einer bildhaften Veranschaulichung dieses Effekts von Edward N. Lorenz am Beispiel des Wetters, in der er fragt ob der Flügelschlag eines Schmetterlings in Brasilien einen Tornado in Texas auslösen kann.

Die heutige Antwort hierauf müßte wohl ein "Unter Umständen ja" sein. Im Übrigen: Bevor jetzt Schmetterlinge gajagt werden um Stürme zu verhindern muß man sich klar machen, das der Flügelschlag den Sturm auch verhindern könnte. Es gibt Situationen in denen Wetterlagen nicht stabil sind und sich in völlig verschiedene Richtungen entwickeln können. Vergleichbar in etwa mit einer Kugel, die einen Bergrücken entlangrollt. Je nachdem ob sie ein wenig mehr auf der einen oder anderen Seite gestartet wird, wird sie letztendlich eine andere Hangseite hinab rollen. Die Endpositionen unterscheiden sich dramatisch, obwohl die Startpositionen vielleicht nur Millimeter voneinander entfernt waren. Bezogen auf das Wetter bedeuted das nun, daß eine kleine Veränderung in der Ausgangssituation, so zum Beispiel ein Hochdruckgebiet, daß sich 1 mm weiter nördlich befindet das Wetter in einem anderen Teil des Kontinents einige Zeit später völlig verändern kann. Aus genau diesem Grund kann der Wetterbericht auch nicht beliebig lange Zeiträume vorhersagen, denn umso weiter in die Zukunft gerechnet wird, umso empfindlicher wird das Endergebnis in Bezug auf minimale Variationen der Ausgangsbedingungen. Die Ausgangsbedingungen sind jedoch nicht beliebig genau bestimmbar.

Also was hat das alles mit dem magnetischen Pendel zu tun? Nun gut ich werde es euch zeigen, ich gebe euch die Gleichungen, die Erklärungen, die Algorithmen und die Software aber bevor ich das tue gebe ich euch etwas besseres: Ein Video, welches die Simulation ausführlich erklärt. Folgendes Video wurde von Aldo Cavini Benedetti erzeugt, der mir freundlicherweise die Genehmigung erteilt hat es hier zu zeigen. Er hat mit viel Aufwand den Versuch aufgebaut um zu zeigen, worum es bei der Simulation eigentlich geht und wenn ihr das Video bis zum Ende anschaut werdet ihr sehen, das selbst Katzen das verstehen.


Video: Chaos Theorie erklärt am Beispiel des Magnetischen Pendels ; © 2007 Aldo Cavini Benedetti

Übersicht über die Simulation

Im nächsten Abscnhitt werde ich das zugrundeliegende numerische Modell vorstellen und zeigen, wie kleine Änderungen in den Anfangsbedingungen der Simulation die Endposition des Pendels stark beeinflußen können. Dies führt in letzter Konsequenz zu einer Unvorhersehbarkeit der Pendelendposition, da bereits kleinste Änderungen am Anfang das Ergebnis dramatisch beeinflußen.

Kommen wir jetzt zu den Details. Das klassische Modell geht von einem Pendel aus, das von drei Magneten angezogen wird. Jedem Magneten ist eine bestimmte Farbe zugeordnet. Die Magnete sind unterhalb des Pendels auf einer Ebene in einem gleichseitigen Dreieck angeordnet. Das Pendel ist mittig darüber aufgehängt. Die Magnete sind stark genug, um das Pendel niemals in der Mitte zur Ruhe kommen zu lassen. Im folgenden Bild ist das Modell in der Draufsicht dargestellt. Die farbigen Kreise symbolisieren die Magnete, das Kreuz in der Mitte ist der Aufhängepunkt des Pendels. Die Simulation berechnet die Bewegung des Pendels unter dem Einfluß der drei Magnete und der auf das Pendel wirkenden Schwerkraft und Reibungskraft. Infolge von Energieverlusten durch die Reibungskraft wird das Pendel früher oder später über einem der Magnete zur Ruhe kommen. Der Startpunkt des Pendels wird dann mit der Farbe des Magneten eingefärbt, über dem es zur Ruhe kam. Führt man dies für alle Punkte im Simulationsbereich durch, so erhält man eine farbiges Muster aus roten, grünen und blauen Punkten.

Pendulum traces above the magnets
Bild 1 (links): Grafische Darstellung des Experimentes (Bild mit freundlicher Genehmigung von Paul Nylander).
Bild 1 (rechts): Kleine Veränderungen in den Anfangsbedingungen führen zu großen Abweichungen in der Pendelbewegung.

Das Ergebnis einer Berechnung ist in den Bildern 1 und 3 dargestellt. Der weiße Fleck in der Mitte von Bild 3 wird dadurch verursacht, daß das Pendel infolge der Parameterwahl in der Mitte zur Ruhe gekommen ist. Mit der bislang beschriebenen Methode sind allerdings nur dreifarbige Bilder möglich. Es wäre wünschenswert in den Bildern mehr Farben oder zumindest Schattierungen zu verwenden. Die einzige Information, die dafür noch verfügbar ist, ist die Spurlänge. Es scheint also angebracht diese in die Berechnung der Farbe einzubeziehnen. Dies kann mit einer Farbfunktion geschehen, welche die Spurlänge berücksichtigt. Zur Erleichterung der Skalierung wird auch noch die maximale Spurlänge in die Berechnung mit einbezogen. Prinzipiell kann man beliebige Funktionen verwenden, aber erfahrungsgemaß bringen die folgenden gute Ergebnisse:

Color functions
Bild 2: Gegenüberstellung von drei unterschiedlichen Farbfunktionen.

Wie man sieht ist das Ergebnis ziemlich interessant, wenigstens wenn man sich ein wenig fü Chaos-Theorie interessiert. (Wenn nicht frage ich mich wieso du noch liest...) Anfangspositionen, die in längeren Spuren Enden werden in dunkleren Schattierungen dargestellt. Die Software erlaubt die verwendung beliebiger Farbfunktionen. Die Formeln werden mit einem Mathematikparser interpretiert.

Simulation setup overview
Bild 3 (links nach rechts): Simulationsschema; Farbzuweisung nach Magnetindex; Farbzuweisung mit Farbfunktion

Idealisiertes Modell

Wem sich nun die Frage stellt, wie das Pendel mit seinen drei Raumdimensionen und der Bewegung auf einer Kugeloberfläche in ein mathematisches Modell umgesetzt wird, dem sei gesagt: Das ist einfach! Wie heißt es doch so schön? "Was nicht paßt wird passend gemacht" Getreu diesem Grundsatz wird einfach angenommen, das Pendel sei sehr lang. Unter dieser Vorraussetzung kann man die Bewegung der Pendelmasse als Bewegung in einer Ebene über den Magneten idealisieren. Für ein langes Pendel mit geringer Auslenkung kann man weiterhin annehmen, die Rücktriebskraft, welche es zum Zentrum zurück zieht, ist proportional zur Auslenkung. Sie verhält sich also dem Hookesches Gesetz entsprechend. Das alles ist äußerst praktisch, denn nun kann man in zwei Dimensionen rechnen, ohne sich um Drehmomente, Rotationswinkel, Kreuzprodukte und ähnliches Gedanken machen zu müssen.

So viel mehr Aufwand, wäre es zugegebenermaßen auch nicht gewesen die komplette Physik zu implementieren, aber die Realität hat mal wieder einen entscheidenden Nachteil: Das Ergebnis würde auf eine Kugel gemappt werden. Das wäre ein wenig unpraktisch, denn im Gegensatz zu Bildern kann man Kugeln weder an die Wand hängen noch als Desktophintergrund verwenden.

Force distance relationships

Bild 4: Vergleich von Pendelrücktriebskraft mit der von den Magneten verursachten Anziehungskraft.

Für die Berechnungen wird angenommen, die Magnete üben eine Kraft auf das Pendel aus, die proportional zum inversen Quadrat der Pendelentfernung vom Magneten ist. Dies ist ein typischer Zusammenhang, wie er auch vom Gravitationsgesetz oder bei Coloumbschen Gesetz bekannt ist. Das ist allerdings nur die halbe Wahrheit, denn dieses Kraftgesetz würde nur gelten, wenn die Kraft von magnetischen Monopolen verursacht würde und die gibt es nach bisherigem Kenntnisstand nicht. Magnete sind immer Dipole und ein Dipol verursacht eigentlich eine Kraft, die propertional zu 1/r³ ist. Dies wird hier ignoriert, kann aber vom interessierten Leser im Konfigurationsfile selbst eingestellt werden. Weiterhin werden elektrische Effekte, wie zum Beispiel Induktion und dadurch im Pedel selbst verursachte Magnetfelder vernachlässigt.

Die Bewegungsgleichungen

Die Positionen des Pendels ergeben sich durch zweifache Integration der auf das Pendel wirkenden Gesamtbeschleunigung. Die Beschleunigung folgt unmittelbar aus Newtons erstem Bewegungsgesetz. Die Kraft, die notwendig ist, um einen Körper zu bewegen, oder dessen Bewegung zu ändern entspricht dem Produkt aus Masse und Beschleunigung. Damit kann die Beschleunigung direkt berechnet werden.

Newtons first law

Da die Ausgangsposition und die Initialgeschwindigkeit sich aus den Anfangsbedingungen ergeben, muß nur noch die Gesamtbeschleunigungen berechnet werden. Diese ergibt sich aus der Summer der durch die Magnete verursachten Beschleunigungen plus der durch die Pendelrücktriebskraft verursachten Beschleunigung minus der durch die Reibung verursachte Bremswirkung. Aus praktischen Erwägungen wird für alle weiteren Betrachtungen eine Pendelmasse von 1 Masseneinheit angenommen. Es ergeben sich folgende Gleichungen:

Equations of motion

Integration der Pendelspur

Gibt man eine Startposition und eine Startgeschwindigkeit für das Pendel vor, so kann man mit einem passenden Integrationsverfahren die komplette Spur berechnen. Für diese Simulation wird das Beeman-Integrationsschema verwendet. Dieses Schema wird häufig auf Probleme der Molekulardynamik angedwendet und ist der ein zur Verlet-Integration verwandtes Verfahren. Es ist wirklich einfach zu implementieren und sieht als Pseudocode so aus:
  
while tracing pendulum
    position += Velocity
    acceleration = 0

    for all force sources
    acceleration += acceleration_caused_by_source

    if (pendulum is close to source and velocity is small) then
      stop_magnet = source_index
      break
    end if
    end for

    acceleration -= acceleration_caused_by_friction
    velocity     += acceleration
    trace_length += length(velocity)

    store stop_magnet
    store trace_length
  end while

In C++ sieht das ganze dann so aus

for (int ct=0; ct<m_nMaxSteps && bRunning; ++ct)
{
  // compute new position
  pos[0] += vel[0]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[0] - 1.0 / 6.0 * (*acc_p)[0]);
  pos[1] += vel[1]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[1] - 1.0 / 6.0 * (*acc_p)[1]);

  (*acc_n) = 0.0;    // reset accelleration

  // Calculate Force, we deal with Forces proportional
  // to the distance or the inverse square of the distance
  for (std::size_t i=0; i<src_num; ++i)
  {
    const Source &src( m_vSources[i] );
    r = pos - src.pos;
    if (src.type==Source::EType::tpLIN)
    {
      //---------------------------------------
      // Hooke's law:          _
      //         _             r         _   
      //  m * a = - k * |r| * --- = -k * r
      //                      |r|
      //
      (*acc_n)[0] -= src.mult * r[0];
      (*acc_n)[1] -= src.mult * r[1];
    } 
    else
    {     
      //---------------------------------------
      // Magnet Forces: _
      //      _         r
      //  m * a = k * -----
      //               |r³|
      //
      double dist( sqrt( sqr(src.pos[0] - pos[0]) +  
                         sqr(src.pos[1] - pos[1]) + 
                         sqr(m_fHeight) ) );
      (*acc_n)[0] -= (src.mult / (dist*dist*dist)) * r[0]; 
      (*acc_n)[1] -= (src.mult / (dist*dist*dist)) * r[1]; 
    }

    // Check abort condition
    if (ct>m_nMinSteps && abs(r)<src.size && abs(vel)<m_fAbortVel)
    {
      bRunning = false;
      stop_mag = (int)i;
      break;
    }
  }   // for all sources

  //--------------------------------------------------------------
  // 3.) Friction proportional to velocity
  (*acc_n)[0] -= vel[0] * m_fFriction;
  (*acc_n)[1] -= vel[1] * m_fFriction;    

  //--------------------------------------------------------------
  // 4.) Beeman integration for velocities
  vel[0] += dt * ( 1.0/3.0 * (*acc_n)[0] + 5.0/6.0 * 
                 (*acc)[0] - 1.0/6.0 * (*acc_p)[0] );
  vel[1] += dt * ( 1.0/3.0 * (*acc_n)[1] + 5.0/6.0 * 
                 (*acc)[1] - 1.0/6.0 * (*acc_p)[1] );

  //--------------------------------------------------------------
  // 5.) flip the acc buffer
  tmp = acc_p;
  acc_p = acc;
  acc   = acc_n;
  acc_n = tmp;
}

Bedienung der Software


Die Anwendung benötige die Laufzeitumgebung von Visual Studio 2005 um zu starten. Die Benutzeroberfläche zu erklären ist einfach, denn es gibt (fast) keine. Alle Simulationsparameter werden aus einem Konfigurationsfile eingelesen, das als Parameter beim Start übergeben werden muß. Das Programm besteht aus einem einzigen Fenster, die Berechnung startet sofort. Da das Warten auf ein Ergebnis recht langweilig sein kann habe ich ein wenig Interaktivität eingebaut. Wannimmer die Maus bewegt wird berechnet das Programm die zur Mausposition gehörende Pendelspur. Das ist zwar nur spielerei, kann aber auch ganz hilfreich sein, um manuell bestimmte Gebiete der Simulation abzutesten. Hier kommen wir auch gleich wieder auf den Schmetterlingseffekt zurück, denn es ist schon interessant zu sehen, wie stark sich die Spuren von benachbarten Startpositionen unterscheiden können.

Durch drücken der rechten Maustaste kann angezeigt werden, welche Pendelspur gerade berechnet wird. Nochmaliges drücken deaktiviert diese Anzeige wieder.

Application window

Bild 5: Die Simulation in Aktion.

Die Anwendung nutzt die Vorteile von Mehrkernsystemen durch verteilen der Rechenaufgaben auf die einzelnen Kerne der CPU.

Format der Konfigurationsdatei

Das Format der Konfigurationsdatei entspricht dem einer einer Windows INI-Datei. Es können folgende Parameter verwendet werden:

Sektion Schlüssel Beschreibung
[FIELD] COLS Anzahl der Spalten des Simulationsgebietes. Dieser Parameter definiert die räumliche Diskretisierung in X-Richtung.
HEIGHT Anzahl der Reihen des Simulationsgebietes. Dieser Parameter definiert die räumliche Diskretisierung Y-Richtung.
SIM_WIDTH
(optional)
Breite der Simulationsdomain in metern. Wird auf COLS gesetzt, wenn unspezifiziert.
SIM_HEIGHT Länge der Simulationsdomain in Metern. Wird auf ROWS gesetzt, wenn unspezifiziert.
WIN_WIDTH
(optional)
Breite des Ausgabefensters in Pixel. Benutze diesen Parameter um das Ausgabefenster zu verkleinern, wenn die Simulationsgröße über die Bildschirmgröße heraus geht. Ist der Wert ungesetzt, so wird er mit dem Wert von COLS belegt.
WIN_HEIGHT
(optional)
Die Höhe des Ausgabefensters in Pixeln. Dieser Parameter kann verwendet werden, um die Höhe des Ausgabefensters anzupassen, wenn Bilder berechnet werden, die Größer als der Bildschirm sind. Wird der Parameter weggelassen, do wird er mit dem Wert von ROWS belegt.
[SIMULATION] THREADS
(optional)
Legt die Anzahl der Threads fest, die von der Simulation erzeugt werden fest. Wenn der Parameter ausgelassen wird, dann wird automatisch die Anzahl der Kerne verwendet. Jeder Thread läuft auf einem eigenen Kern. Dadurch kann die Geschwindigkeit auf Mehrkernsystemem deutlich erhöht werden. (ca. 40% mit einem Dual-Core Rechner).
FRICTION Der Reibungskoeffizient. Die Reibungskraft ist direkt proportional zur Geschwindigkeit des Pendels wobei dieser Wert als Proportionalitäskonstante dient.
PEND_HEIGHT Die Höhe des Pendels über der Ebene in der die Magnete angeordnet sind.
DELTA_T Die Integrationsschrittweite. Eine kleinere Schrittweite bedeuted eine genauere Simulation. Dieser Parameter entspricht dem Parameter h in der Beeman-Integrationsformel.
MIN_STEPS Minimale Anzahl an Schritten, die eine Pendelspur haben muß bevor die Abbruchbedingungen überprüft werden.
MAX_STEPS Maximale Anzahl an Schritten, die eine Pendelspur haben darf. Wird diese Anzahl erreicht bricht die Berechnung ab, auch wenn das Pendel noch nicht über einem Magneten zur Ruhe gekommen ist.
ABORT_VEL Wenn die Pendelgeschwindigkeit unter diesen Wert fällt wird die Verfolgung der Spur beendet und die Berechnung des nächsten Pixels wird gestartet.
COLOR_SCHEME Eine Gleichung die verwendet werden soll um die Farbe eines Pixels in Abhängigkeit der Gesamtspurlänlge zu berechnen. Die Farbe eines Pixels wird durch die Länge der gerade berechneten Spur und die Länge der Längsten Spur bestimmt. Es muß sich bei diesem Parameter um eine Funktion handeln, der die Variablen len und max_len (siehe Bild 2) enthält.
zum Beispiel: 1/(exp(0.000001*(len*len)))
BATCH_MODE
(optional)
Dieses Flag aktiviert den Batch-Modus. Im Batch-Modus beendet sich die Simulation selbstständig nach einer Berechnung. Dieses Flag kann benutzt werden um mittels Shell-Skripten Animationen zu erstellen.
[SOURCE ...]
... steht für einen Index
TYPE Gibt die Art der Quelle an. Erlaubt sind folgende Schlüsselwörter: INV_SQR oder LINEAR. Die von der Quelle ausgehende Kraft ist entweder proportional zur Entfernung des Pendels oder proportional zum inversen Quadrat der Pendelentfernung. Der erste Fall beschreibt ein Verhalten ähnlich dem Hookesches Gesetz, der zweite Fall beschreibt eine Quelle ähnlich dem Coulombschen Gesetz (Rein physikalisch gesehen haben wir es hier allerdings mit Magneten zu tun, nicht mit elektrischen Ladungen also ist das Coulombsche Gesetz hier nur Beispielhaft zu verstehen).
COLOR Die Farbe der Quelle. Kommt das Pendel über der Quelle zu stehen, dann wird seine Startposition mit dieser Farbe eingefärbt.
RAD Die Positionen der Quellen werden in Polarkoordinaten angegeben. Dies ist der Abstand der Quelle vom Simulationszentrum aus gesehen.
THETA Die Positionen der Quellen werden in Polarkoordinaten angegeben. Dies ist der Winkel, den die Quelle vom Simulationszentrum aus gesehen einnimmt.
MULT Multiplikator für die Stärke der Kraft. Je größer dieser Wert ist, umso stärker wirkt die Kraft auf das Pendel.
SIZE Größe der Quelle. Die Abbruchbedingungen werden nur überprüft, wenn das Pendel sich näher als so viele Pixel an der Quelle befindet.

Beispielekonfigurationsdateien befinden sich im data Unterverzeichnis. Wenn du eine eigene Konfiguration erstellen möchtest, dann fange am besten mit einer bestehenden Datei an und modifiziere diese.

Verwandte Links


Die folgenden links führen zu Webseiten, die ergänzende Informationen zur Magnetpendelsimulation oder zu Fraktalen im Allgemeinen bereit stellen.

Gallerie


Abschließend möchte ich einige Bilder präsentieren, die mit dem hier vorgestellten Programm berechnet wurden. Die Simulationen unterscheiden sich durch die Anzahl der Magnete, deren Stärke, Position und des Punktes der Pendelaufhängung.

Gallery
Bild 6: Eine Auswahl von Bildern berechnet mit unterschiedlichen Parametern.

Versionen


Version 1.00: 30/10/2006 - Erstveröffentlichung

Eigenschaften:

Version 1.01: 19/11/2006

Änderungen:


Creative Commons License
Autor: Ingo Berg; Dieser Text und die dazugehörigen Mediendateien stehen, sofern nicht anderweitig angegeben, unter der Creative Commons Namensnennung 3.0 Deutschland Lizenz.