Die Pendelsimulation

Die folgende Animation zeigt eine Echtzeit-Simulation in Aktion. Mit den Schiebereglern können während die Animation läuft diverse Parameter verändert werden und deren Wirkung am Modell beobachtet werden. Eine Beschreibung aller Parameter ist weiter unten.

Die Animation ist physikalisch korrekt. Es wird kein vereinfachtes Modell oder eine vorberechnete Bildsequenz verwendet. Es wird ein physikalisch-mathematisches Modell erstellt, davon die Bewegungsgleichungen abgeleitet und diese werden numerisch gelöst und in der Animation angezeigt.

Mathematisches Modell des Pendels

Kräftediagramm des Pendels

Für die Simulation des Pendels sind nur die tangentialen Kräfte $\vec D$ und $\vec G_\mathrm{T}$ relevant, da sich das Pendel in der Achse der Pendelschnur immer auf einer festen Kreisbahn bewegt, in dieser Achse also nicht frei ist. Das bedeutet, dass die Kraft der Pendelschnur immer gerade so gross wird, dass die Komponente $\vec G_\mathrm{S}$ der Gravitation plus die Zentrifugalkraft $\vec Z$ gerade aufgehoben werden: $\vec S = - (\vec G_\mathrm{S} + \vec Z)$. Aber wie gesagt, auf den Auslenkwinkel $\varphi$ haben die Komponenten parallel zur Pendelschnur keinen Einfluss.

Die tangentiale Komponente der Gewichtskraft ist:

Die Reibungskraft $\vec D$ wirkt immer tangential entgegen der Bewegungsrichtung der Pendelmasse und ist abhängig von der tangentialen Geschwindigkeit $\vec v$ der Pendelmasse:

Die Summe aller tangentialen Kräfte ist:

Die tangentiale Beschleunigung $a$ bzw. Winkelbeschleunigunug $\ddot \varphi$ der Pendelmasse ist nach Newton:

Setzen wir die einzelnen Kräfte in obige Gleichung ein erhalten wir die Differetialgleichung:

und zusammengefasst:

Bewegungsgleichung des Pendels

Diese Differentialgleichung ist die sog. Bewegungsgleichung des Pendels. Aus der aktuellen Auslenkung $\varphi$, der aktuellen Winkelgeschwindigkeit des Pendels $\dot \varphi$ und den diversen anderen Parametern wie Pendellänge $l$, Pendelmasse $m$, Gravitationsbeschleunigung $g$ und Dämpfungsfaktor $d_2$ lässt sich also die Winkelbeschleunigung $\ddot \varphi$ berechnen.

Die Winkelgeschwindigkeit $\dot \varphi$ ist das Integral der Winkelbeschleunigung $\ddot \varphi$ und der Winkel $\varphi$ ist das Integral der Winkelgeschwindigkeit $\dot \varphi$:

Da die Differentialgleichung (6) nicht linear ist, können diese Integrale nur numerisch gelöst werden. Dazu wird für den aktuellen Zustand des Pendels ($\varphi$ = a und $\dot \varphi$ = da in der Animation) die aktuelle Beschleunigung nach (6) berechnet. Durch numerisches integrieren wird für den nächsten Zeitpunkt $t + \Delta t$ die neue Winkelposition $\varphi$ (a) und Winkelgeschwindigkeit $\dot \varphi$ (da) berechnet. Dies ergibt dann die neuen Ausgangsdaten für den nächsten Iterationsschritt usw.

Numerisches Lösen der Bewegungsgleichung

Beim numerischen Lösen einer Bewegungsgleichung werden Schritt für Schritt aus den momentanen Randbedingungen (z.B. aktuelle Position s(t) und Geschwindigkeit v(t) eines Teilchens) die aktuellen Kräfte berechnet. Aus der Summe ∑F = f(v,s,t) (d.h. die Summe der Kräfte ist eine Funktion aller Parameter) ergibt sich dann die Beschleunigung a_i zu einem gewissen Zeitpunkt t_i. Daraus lässt sich für den nächsten Schritt t_i+1 die neue Geschwindigkeit v_i+1 und die neue Position s_i+1 durch numerisches Integrieren berechnen.

Beschleunigung a(t), Geschwindigkeit v(t) und Weg s(t) sind mathematisch miteinander wiefolgt verknüpft:

Aus der Bewegunugsgleichung kennen wir die Beschleunigung a(t) (die Winkelbeschleunigung im Falle des Pendels) und können durch Integrieren daraus Geschwindigkeit v(t) (Winkel-Geschwindigkeit des Pendels) und Position s(t) (Auslenkwinkel des Pendels) berechnen:

Numerische Integration

Wenn man den Verlauf der Beschleunigung a(t) grafisch aufzeichnet, erhält man eine bestimmte Kurve a = f(t). An jeder Stelle t_i kann die zugehörige Beschleunigung a(t_i) abgelesen werden 1. Die Fläche unter der Kurve a(t) vom Zeitpunkt t = 0 bis t_i entspricht gerade der Geschwindigkeit v(t_i) zu diesem Zeitpunkt, abzüglich der Startgeschwindigkeit (Integrationskonstante) v_o 2.

Wenn man die Beschleunigung a(t_i) zum Zeitpunkt t_i kennt, kann man daraus die ungefähre Geschwindigkeitsänderung Δv bis zum nächsten Zeitpunkt t_i+1 berechnen. Wenn wir annehmen, die Beschleunigung ändere sich im Zeitraum Δt nicht, so ist die Geschwindigkeitsänderung Δv = a(t_i) · Δt 3.

Solange man Δt 4 klein genug wählt, ist dies eine hinreichend genaue Annäherung. Man kann jedoch die Berechnung noch verbessern, wenn man die Steigung der Kurve a(t) berücksichtigt. Die Änderung der Beschleunigung Δa im Zeitraum Δt 5 kann man abschätzen, indem man den vorletzten Beschleunigungswert a(t_i−1) 6 mit berücksichtigt:

Damit ergibt sich ein genauerer Wert für die Geschwindigkeitsänderung und Geschwindigkeit:

Analog kann aus der Geschwindigkeit der Weg berechnet werden:

Die obigen Formeln sind exakt, wenn die Beschleunigung konstant oder linear ist und die Geschwindigkeit linear zu- oder abnimmt.

Algorithmus

Mit obigen Formeln lässt sich nun ein Programm schreiben, welches Beschleunigung, Geschwindigkeit und Weg Schritt für Schritt berechnet:

var dt = 0.1;   // Delta t
var da, dv, dt; // Zwischenwerte

// Startbedingungen
var t = 0;      // Zeit
var v = Vstart; // Anfangsgeschwindigkeit
var s = 0;      // Anfangsposition

// Initiale Berechnung
var a = Beschleunigung( v, s, t );
var aOld = a;   // Beschleunigung der letzten Iteration
Zeichne( a, v, s, t );  // Resultate anzeigen

while (true) {
  // Änderungen für nächsten Schritt berechnen
  da = a - aOld;
  dv = (a + da/2) * dt;
  ds = (v + dv/2) * dt;
  // neue Werte berechnen
  v = v + dv;
  s = s + ds;
  t = t + td;
  aOld = a;
  // neue Beschleunigung berechnen
  a = Beschleunigung( v, s, t );
  Zeichne( a, v, s, t );
}

In der Funktion Beschleunigung( v, s, t ) wird die neue Beschleunigung aufgrund der aktuellen Situation über die Berechnung aller Kräfte berechnet. Der obige Algorithmus ist im Sim-Modul implementiert.

Echtzeit-Simulationen

Die Simulation eines Modells geschieht in kleinen Zeiteinheiten Δt. Je kleiner diese Einheit gewählt wird, umso genauer wird die Simulation in der Regel. In jeder Zeiteinheit wird der neue Zustand des Modells aufgrund des aktuellen Zustandes anhand von physikalischen Formeln berechnet. Beim Pendel werden zum Beispiel anhand des aktuellen Zustandes alle Kräfte am Pendel berechnet. Daraus resultiert die aktuelle Beschleunigung des Pendels. Durch Numerische Integration kann daraus die neue Geschwindigkeit und Position des Pendels nach der Zeit Δt berechnet werden.

Animationen sind eine Sequenz von Bildern, Frames genannt. Zwischen Frames vergeht eine bestimmte Zeit. Zum Beispiel versuchen Browser 60 Frames pro Sekunde anzuzeigen. Das heisst, zwischen den Frames vergehen ca. 16,7 ms. Damit das Modell am Ende dieser Zeit den Zustand anzeigt, wie er nach Ablauf dieser 16,7 ms in der Realität aussehen würde, müssen entsprechend viele Berechnungsschritte am Modell in dieser Zeit durchgeführt werden, damit die Darstellung synchron verläuft. In der Regel benötigt der Computer keine 16,7 ms um 16,7 ms Realzeit des Pendels zu simulieren. Es braucht aber auch noch Rechenzeit für das Zeichnen der neuen Grafik. Das Zeichnen ist in dieser Animation wesentlich rechenintensiver als die Berechnung des neuen Zustandes des Pendels nach 16,7 ms.

Das Problem bei Echtzeit-Animationen ist, dass die Zeit zwischen zwei Frames variabel ist. Der Computer erledigt neben der Simulation noch viele andere Aufgaben, die auch unvorhersehbar viel Rechenzeit benötigen. Ebenso ist die Rechenzeit der Simulation variabel. Man muss nun dafür sorgen, dass die Berechnung der Simulation mit der zur Verfügung stehenden Zeit zwischen zwei Frames auskommt. Dies kann durch die Wahl von Δt der Simulation beeinflusst werden. Bei zu grossen Werten wird die Simulation ungenau und kann im Extremfall komplett aus dem Ruder laufen. Bei zu kleinen Werten wird mehr Rechenzeit zwischen zwei Frames benötigt, als zur Verfügung steht. In diesem Fall kann die Simulation nicht in Echtzeit ablaufen. Damit der Browser in diesem Fall nicht blockiert, wird die Simulation entsprechend verlangsamt. Wenn man die Simulation im Zeitlupen-Tempo ablaufen lässt, muss weniger berechnet werden, sodass die Zeitspanne zwischen zwei Frames ausreichen kann.

Die Anzahl Berechnungsschritte zwischen zwei Frames ist in der Regel nicht konstant. Um eine flüssige Animation zu erhalten, müssen immer genau soviele Schritte berechnet werden, dass die Simulationszeit mit der Realzeit synchron bleibt.

Das Sim-Modul

Ich habe ein Modul Sim programmiert, das sich um all diese Details kümmert. Man kann sogar einstellen, ob die Simulation schneller oder langsamer als real ablaufen soll. Es können statistische Berechnungen aktiviert werden, damit man überprüfen kann, wie stark der Rechner durch die Simulation belastet ist und ob Bilder ausfallen. Das Modul nimmt auf Wunsch auch die Berechnung der numerischen Integration vor, mit der man die Bewegungsgleichung einer Simulation löst.

Das Prinzip von Sim ist folgendes:

Initialisieren der Animation, Berechnung des Startzustand und zeichnen diesen Zustandes
Start der Simulation
Bei jedem Frame
  Berechne die Zeit seit dem letzten Frame
  Berechne wieviele Schritte n simuliert werden müssen, 
    um die gewünschte Simulationszeit zwischen zwei Frames zu erreichen
  Für jeden dieser n Schritte
    Rufe die Berechnungsfunktion des Modells auf
  Rufe die Zeichenfunktion des Modells auf

Das Modell muss also drei Funktionen zur Verfügung stellen:

• Reset: Setzt das Modell auf einen Anfangszustand zurück
• NextStep: Berechnet den nächsten Simulations-Schritt (numerisches Integrieren)
• Draw: Zeichnet den aktuellen Simulations-Zustand

Diese drei Funktionen werden im Sim installiert und automatisch in der nötigen Reihenfolge und Anzahl aufgerufen. In der Zeichenfunktion kann falls gewünscht auch die Statistik des Sim ausgegeben werden.

Rundungsfehler

Numerische Simulation sind nur eine Annäherung an die Realität und müssen je nach Modell mit mehr oder weniger starken Rundungsfehlern kämpfen. Je länger eine Simulation läuft, umso grösser wird die Abweichung zur Realität. Allerdings kann die Realität sowieso auch nur angenähert werden, da man kaum alle erdenklichen Einflüsse in die Simulation einbeziehen kann oder möchte. Im Modell des Pendels wirken sich Rundungsfehler meist als schwache Dämpfung aus. Eventuell könnte sich das Modell aber auch aufschaukeln, wenn die Dämpfung d2 = 0 gesetzt wird.

Rundungsfehler können verkleinert werden, indem die Simulation mit kleineren zeitlichen Abständen Δt rechnet. Die obige Simulation rechnet mit Δt = 1 ms. Es werden also pro Sekunde 1000 Simulationsschritte berechnet! Das sind ca. 16..17 Schritte pro Frame. Verkleinert man Δt um die Genauigkeit zu steigern, wird entsprechend mehr Rechenzeit pro Frame benötigt. Reicht die Rechenzeit zwischen zwei Frames nicht mehr aus, wird die Simulation verlangsamt, läuft also nicht mehr in Echtzeit ab.

Code der Pendel-Simulation

Im folgenden Code habe ich nur die für die Simulation relevanten Teile aufgelistet. Das Zeichnen der Kurven habe ich weggelassen.

#INCLUDE JsGraph.inc
#INCLUDE Sim.inc
#INCLUDE ControlPanel.inc

<jscript>
// Pendulum Model

var Pendulum = {
  v0: 0,          // initial angular speed
  state: {        // current angular acceleration, angular speed and angle
    OldAccel: 0,
    Accel: 0,     // dda = angular acceleration
    Speed: 0,     // da  = angular speed
    Pos: 0,       // a   = angle
  },
  d1: 0,          // linear damping factor
  d2: 0,          // quadratic damping factor
  m: 1,           // mass of pendulum (does not matter)
  l: 1,           // length of pendulum
  g: 9.81,        // gravitational acceleration
  Fg: 0,          // gravitational force
  Ft: 0,          // sum of all tangential forces
  Fd: 0,          // drag force
  
  graph: null,    // link to graphic model JsGraph
  sim: null,      // link to Sim

  Create: function( sim ) {
    // set link to Sim and create JsGraph object on the webpage
    this.sim = sim;
    var me = this;
    this.graph = NewGraph2D( {
      Id: 'Pendulum-Graph',
      Width: '100%',
      Height: '60%',
      DrawFunc: function(g)  { me.Draw(g,sim);  },
      OnClick: function(e,g) { sim.Pause(true); },
      AutoReset: false,
      AutoClear: false,
    } );
    this.d1 = 0.0;
    this.d2 = 0.5;
  },

  Reset: function( sim ) {
    // this function is called from Sim on the start or on each reset of the animation
    // it computes the start conditions of the animation
    var state = this.state;
    state.Pos = 90*Math.PI/180;
    state.Speed = this.v0;
    state.Accel = 0;
    this.sim.InitStates( state, 0, 0 );
  },

  Update: function( sim ) {
    // this function is called from Sim for each simulation step to compute the next delta t step
    var state = this.state;
  
    // compute all forces:
    // * gravity
    this.Fg = this.m * this.g;
    // * drag
    this.Fd = this.l * (this.d1 * Math.abs( state.Speed ) + this.d2 * state.Speed * state.Speed);
    // * tangential component of these forces
    var vdir = state.Speed >= 0 ? -1 : 1;
    this.Ft = -this.Fg * Math.sin( state.Pos ) + vdir * this.Fd;

    // acceleration of pendulum mass from sum of forces: a = F / m
    // * resulting angular acceleration depends on pendulum length: angular_accel = a / length
    state.Accel = this.Ft / this.m / this.l;

    // numerical integration of acceleration using Sim gives new speed and position of pendulum
    this.sim.CompNewStates( state, 0, 0 );

    // restrict angle to +/- PI
    state.Pos = (state.Pos + 3*Math.PI) % (2*Math.PI) - Math.PI;
  },

  Draw: function( g, sim ) {
    // this function is called from Sim or from browser to redraw the new pendulum state 

    g   = g   || this.graph;
    sim = sim || this.sim;
    var state = this.state;

    g.Reset();
    g.MapWindow( 0, 0.5, 2.2, 0, 0 );

    // compute pendulum position
    var x = this.l * Math.sin( state.Pos );
    var y = -this.l * Math.cos( state.Pos ) + 1;
  
    // draw zero line
    g.SetLineAttr( '#ddd', 1 );
    g.Line( 0, 1, 0, -0.5 );

    // draw pendulum
    g.SetLineAttr( 'black', 2 );
    g.Line( 0, 1, x, y );
    g.SetMarkerAttr( 'Circle', 20, 'black', 'yellow', 2 );
    g.Marker( x, y );

    // draw force Ft
    var ax = 0.025 * this.Ft * Math.cos( state.Pos );
    var ay = 0.025 * this.Ft * Math.sin( state.Pos );
    g.SetMarkerAttr( 'Arrow1', 8, 'black', 'red', 1 )
    g.Arrow( x, y, x+ax, y+ay, 1, 3 );

    // display Sim statistic values
    g.SelectTrans( 'viewport' );
    g.SetTextAttr( 'Arial', 16, 'black', 'normal', 'normal', 'left', 'top', 4 );
    g.Text( 'FPS = ' + sim.Fps.toFixed(0), 0, 0 );
    g.Text( 'CPU = ' + sim.CpuLoadFrameAvg.toFixed(3), 0, 20 );
    g.Text( 'SimTime = ' + sim.SimulTime.toFixed(1), 0, 40 );

  },

};

// Create Sim object and install pendulum functions

var PendulumSim = new Sim( {
  EnableStatistics: true,
  TimeStep:      0.001,
  TimeSpeed:     1,
  SimObj:        Pendulum,
  ResetFuncs:    function(sim) { sim.SimObj.Reset(sim); },
  TimeStepFuncs: function(sim) { sim.SimObj.Update(sim); },
  FrameFuncs:    function(sim) { sim.SimObj.Draw(null,sim); },
} );

// create pendulum and pendulum animation canvas

Pendulum.Create( PendulumSim );

// run animation

PendulumSim.Run( true );

// UpdateAll is called to reflect changes of the model to the sliders

function UpdateAll() {
  ControlPanels.Update();
}

// install UpdateAll() as a reset function, called by the Sim to reset all sliders

PendulumSim.AddResetFunc( function(sim){ UpdateAll(sim); } );

// create sliders that directly control animation parameters

ControlPanels.NewSliderPanel( { 
  Name: 'Pandulum-Sliders',
  ModelRef: 'Pendulum',
  NCols: 1, 
  ValuePos: 'left',
  OnModelChange: UpdateAll, 
  Format: 'fix0', 
  Digits: 2,
  ReadOnly: true, 
  PanelFormat: 'InputShortWidth'

} ).AddValueSliderField( {
  Name: 'd2',
  Color: 'red',
  Min: 0,
  Max: 2,

} ).AddValueSliderField( {
  Name: 'v0',
  Color: 'green',
  Min: -20,
  Max: 20,

} ).AddValueSliderField( {
  Name: 'g',
  Color: 'orange',
  Min: 0,
  Max: 20,

} ).AddValueSliderField( {
  Name: 'l',
  Color: 'blue',
  Min: 0.2,
  Max: 1.5,

} ).AddValueSliderField( {
  Name: 'PendulumSim.TimeSpeed',
  Label: 'Speed',
  Color: 'black',
  Min: 0.25,
  Max: 4,

} ).Render();

UpdateAll();
</jscript>