Einführung in die Numerik

Die vielleicht wichtigste Gleichung in der Physik ist das Hook'sche Gesetz: \[ \overrightarrow{F} = - k \cdot \overrightarrow{x} \] .
Wir nehmen diese Gleichung (die Gleichung eines sogenannten harmonischen Oszillators) zum Ausgang, um zwei wichtige "numerische" Verfahren kennenzulernen.  Mit Hilfe dieser Verfahren (Euler-Algorithmus, Euler-Cromer-Algorithmus) werden wir dann in der Lage sein, sogenannte Differenzialgleichungen numerisch lösen zu können.  Das wiederum bedeutet, dass wir fortan Java nutzen können, um viele Beispiele aus der Physik genauer zu untersuchen. Es ist sinnvoll, diese Differentialgleichungen zu verstehen, da sie in der Physik so häufig wie 'Sand am Meer' vorkommen. Übrigens betrachten wir hier nur den Idealfall (keine Reibung), d.h. wenn das Federpendel einmal ausgelenkt und losgelassen wird, dann wird es nie wieder zur Ruhe kommen.

Das Hook'sche Gesetz

Ut tensio, sic vis (1678)

Das Hook'sche Gesetz wird in der  Physik behandelt. Wir wollen es hier nur so weit in Erinnerung rufen (bzw. in eine Form bringen), wie es für die Numerik nützlich ist. 

Das Hook'sche Gesetz lautet: \[ \overrightarrow{F} = - k \cdot \overrightarrow{x} \] 

Hierbei ist:

  • x die Strecke (in m), um die die Feder (von ihrer Ruhelage aus gemessen) ausgedehnt wird.

  • F die rücktreibende Kraft (in Newton N). Die rücktreibende Kraft ist der Kraft, die zur Auslenkung der Feder benötigt wird, entgegengesetzt, deshalb ist ein Minuszeichen nötig.

  • k die Federkonstante (Einheit: N/m). Sie hängt von der Federsorte ab. Ist also experimentell zu bestimmen. Eine Federkonstante von z.B. k = 3 N/m bedeutet, dass sich die Feder um einen halben Meter dehnt, wenn sie mit einer Kraft von 1.5 N belastet würde.

 

Die Kraft, die nötig ist, um die Feder zu dehnen, wird umso größer je weiter die Feder gedehnt werden soll. Wir gehen davon aus, dass wir die Feder nur so weit dehnen wollen, dass sie, ohne beschädigt zu werden, wieder in den Ausgangszustand gelangt, ansonsten gilt das Hook'sche Gesetz natürlich nicht. Diesen Sachverhalt schreiben wir mathematisch so:

  • F~ x  →  F/x = const =k →  F = k x

Entsprechend gilt für die rücktreibende Kraft:

  • F~ x  →  F/x = const =-k →  F = - k x

Wir schreiben diese Gleichung jetzt in eine Form, wie sie für die Numerik benötigt wird. Dazu muss man allerdings das 2.Newtonsche Axion (F=m a) kennen, außerdem sollte man wissen, wie die Geschwindigkeit und Beschleunigung definiert ist.

 
Die Orginalversuchsanordnung
von Robert Hook (1635-1704).

Die Ableitung nach der Zeit kennzeichnen wir mit einem Punkt über der Variablen.


Damit schreiben wir für die Geschwindigkeit:

\[ v= \frac{dx}{dt}= \dot x \]

Für die Beschleunigung erhalten wir jetzt:

\[a= \frac{dv}{dt}= \frac{d \dot x}{dt}= \ddot x \]

Damit lautet nun das Hook'sche Gesetz:

\[ F= m a = m \ddot x = - k x \rightarrow  \ddot x = - \frac{k}{m} x \]



Aus der ursprünglichen einen Gleichung sind jetzt zwei geworden, dafür liegen sie jetzt in einer Form vor, wie es für die Numerik benötigt wird.


Das Hook'sche Gesetz (anders geschrieben):

\[  \frac{dv}{dt}= - \frac{k}{m} x \]

mit

\[  \frac{dx}{dt}=  v \]


      


Federpendel: Eine Animation der Gleichung  F  = - k x
von Oleg Alexandrow

 

Was ist ein Phasendiagramm?

 
Phasendiagramm - ein Diagramm mit der Geschwindigkeit (Impuls) als Ordinate und dem Ort als Abszisse.

Phasenpunkt - ein Punkt in der x-v-Ebene zu einem bestimmten Zeitpunkt t.

Trajektorie - eine Kurve, die sich ergibt, wenn der Phasenpunkt sich in der Zeit bewegt.		 Wenn sich ein Teilchen entlang der x-Achse bewegt, wird  dessen Bewegung häufig als ein Graph der Koordinate (wir nennen sie hier x) oder aber der Geschwindigkeit (v=dx/dt)  als Funktion der Zeit dargestellt.
Eine andere aber überaus nützliche Darstellungsweise erhält man, wenn man die  Geschwindigkeit v des Teilchens in Abhängigkeit von der Position x zu jedem Zeitpunkt t graphisch darstellt. Für einen bestimmten Zeitpunkt t erhält man
so einen Punkt (x(t),v(t)), den man als Phasenpunkt bezeichnet. Dieser Punkt wird sich im Laufe der Zeit weiter bewegen. Die Kurve die sich dabei ergibt, bezeichnet man als Trajektorie und die Fläche (b.z.w. Raum), in der sich das Teilchen bewegt als Phasenfläche (Phasenraum). Das Diagramm als ganzes nennt sich dann Phasendiagramm. Solche Diagramme spielen eine wichtige Rolle in allen Forschungsrichtungen der nichtlinearen Dynamik (z.B. Chaosphysik).		  

 

 

 

 

 

 

Wie sieht das Phasendiagramm für den harmonischen Oszillator (Federpendel) aus?

  
Die potentielle Energie des Federpendels ist:
Epot = 1/2 kx²

Die kinetische Energie des Federpendels ist:
Ekin = 1/2 m v²

Für den Fall des ungedämpften und sich ohne Reibung bewegten Federpendels hört die Bewegung niemals auf. Permanent wandelt sich potentielle Energie in kinetische Energie um (und umgekehrt). Das bedeutet aber, dass die Gesamtenergie des Systems eine Konstante sein muss. Die kinetische Energie ist  Ekin = 1/2 m v². Die potentielle Energie des Federpendels ist Epot = 1/2 kx². Demnach gilt:


 \[  \frac{1}{2} m v^2 + \frac{1}{2} k x^2 = E = const.  \]

 

 





Das ist aber die Gleichung für eine Ellipse in x und v, d.h. das Phasendiagramm des Federpendels ist damit eine Ellipse. Wenn sich das Teilchen mit der Zeit bewegt, bewegt sich der Phasenpunkt entlang einer ellipsenförmigen Trajektorie. Wenn das Teilchen sich einmal herumbewegt hat und wieder am Ausgangspunkt angekommen ist, dann ist eine Periode vorbei und es geht von vorne los.



Für jede periodische Bewegung gilt, dass das Phasendiagramm geschlossene Trajektorien besitzt.

 

Das Euler-Verfahren zur numerischen Lösung von Differentialgleichungen

Diskretisieren die unendlich kleine
Größe, z.B. dt wird durch eine kleine
aber endliche Größe (hier Δt)  ersetzt.

Anfangsbedinungen:
Jede numerische
Lösung setzt voraus,
dass man die Anfangsbedingungen kennt. Im vorliegenden Fall bedeutet das, dass der Ort und die Geschwindigkeit am Anfang bekannt sein muss (x(0), v(0)).

Wir demonstrieren das Euler-Verfahren am Beispiel des Federpendels (Hook'sche Gesetz).




Das Differential dt stellt eine unendlich kleine Größe dar. So etwas können wir mithilfe des Computers nicht simulieren. Deshalb ist es notwendig, solche unendlich kleinen Größen durch zwar kleine aber endliche Größen zu ersetzen. Für dt schreiben wir dann  Δt und aus dx wird: dx ≈ x(t+Δt)-x(t),
sowie für dv ≈ v(t+Δt)-v(t), d.h. wir bilden jeweils Differenzen. Hiermit ergibt sich dann für das Hook'sche
Gesetz (bitte unbedingt nachrechnen):



und


Wir sind hier nur daran interessiert, das Verfahren kennenzulernen, deshalb setzen wir im Weiteren k=m=1
Hieraus resultieren dann die von uns verwendeten Gleichungen zur numerischen Lösung:


Mithilfe der ersten Gleichung kann man die Position (Ort) des Teilchen im nächsten Zeitschritt Δt bestimmen. Und mithilfe der zweiten Gleichung ist es möglich, die Geschwindigkeit im nächsten Zeitschritt zu bestimmen. Grundsätzlich
können wir also, wenn wir den Ort und die Geschwindigkeit kennen auch Ort und Geschwindigkeit im nächsten Schritt bestimmen. Trotzdem können wir diese Gleichung momentan noch nicht benutzen, denn erst wenn wir den Ort x(0) und die Geschwindigkeit v(0) am Anfang kennen, können wir diese Gleichung verwenden, um für jeden beliebigen Zeitpunkt Ort und Geschwindigkeit iterativ zu bestimmen.

 


Die Methode init muss zuerst aufgerufen werden, um die Anfangswerte für den Ort, Geschwindigkeit sowie die Schrittweite zu setzen.

Mit der Methode nextEuler() wird jeweils eine Iteration durchgeführt.

Mit getX() und getV() kann nach jeder Iteration der neuberechnete Ort (Geschwindigkeit) ausgegeben
werden.
 
public class Euler {

 private double x;
 private double v;
 private double t;

 public Euler() { 	}

 public void init(double x0, double v0, double timeStep) {
 t = timeStep;
 x = x0 + v0 * t;
 v = v0 - x0 * t;
 }

 public void nextEuler() {
 double x_old = x;
 x = x + v * t;
 v = v - x_old * t;
 }
 public double getX() {
 return x;
 }
 public double getV() {
 return v;
 }
 }

 

 
 
 

Das Euler-Cromer-Verfahren zur numerischen Lösung von Differentialgleichungen

Die Gleichungen für das Euler-Cromer- Verfahren sind fast identisch mt denen des Euler-Verfahrens.
Lediglich bei der Berechnung des Ortes wird die gerade erst berechnete Geschwindigkeit eingesetzt, wohingegen das Euler-Verfahren dort mit dem "alten" Geschwindigkeitswert weiter rechnet. 
 
 
Dennoch im Ergebnis unterscheiden sich die Verfahren stark. Will man z.B. Planetenbahnen berechnen, ist das Euler-Verfahren im allgemeinen unbrauchbar.
 
 

public class EulerCromer {
   private double x;
   private double v;  
  private double t; 

 public EulerCromer() {}

   public void init(double x0, double v0, double timeStep) {
   t = timeStep;
   v = v0 - x0 * t;
   x = x0 + v * t;
 }

 public void nextEulerCromer() { 
                  v = v + x * t;
   x = x - v * t;
 }

 public double getX() {
 return x;
 }

 public double getV() {
 return v; }
 }

 

Vergleich der beiden Verfahren

 
Numerische Lösung für das Federpendel. Das Eulerverfahren ist nicht geeignet  die Trajektorien des Federpendels auch nur eine Periode lang richtig zu simulieren. (Schrittweiter 0.1, k=1, m=1, x0 =0 , v0 = 1.5) - Die Grafik wurde mit GnuPlot von Finn Ole Schröder und Vanessa Sobotta erstellt.