Differentialgleichungen sind Funktionen, in denen neben den Unbekannten auch deren Ableitungen vorkommen. Ich habe dazu eine Einführung geschrieben, bitte diese zuerst lesen. Die formelmäßige Lösung, das Auffinden einer sog. Stammfunktion, gelingt in vielen Fällen nicht. Es gibt jedoch numerische Verfahren, welche eine Näherungslösung berechnen. Das einfachste und damit am leichtesten verständliche ist das Polygonzugverfahren nach Leonhard Euler und Augustin-Louis Cauchy.
Das Verfahren eignet sich nicht für alle Arten von Differentialgleichungen. Insbesondere können sog. partielle Differentialgleichungen damit nicht gelöst werden – in diesen kommen Ableitungen nach mehreren Größen vor. Es ist aber gut geeignet für die in der Astronomie wichtigen Bewegungsanalyse und Bahnberechnung. Zunächst nehmen wir an:
All diese drei Einschränkungen lassen sich umgehen. Ich schreibe dazu weiter unten noch etwas.
Das Polygonzugverfahren ist simpel: Ausgehend vom Startwert berechnen wir die Ableitung, das ist der Anstieg der Kurve. Diese ist ja formelmäßig gegeben. In einer kleinen Umgebung des Startwertes wird sich diese Ableitung auch nur wenig ändern, zumindest bei gutmütigen Funktionen ohne Sprungstellen. Wir sehen daher diese Ableitung als gültig für einen kleinen Schritt in x-Richtung an, dieser kleine Schritt wird oft mit h bezeichnet. Der Anstieg führt vom Startpunkt (x0, y0) zum nächsten Punkt (x1, y1). Für diesen kann man wieder den Anstieg berechnen und zum Punkt (x2, y2) kommen usw.
Die Abbildung macht ein Problem all dieser Näherungsverfahren deutlich: Jeder Rechenschritt erzeugt einen kleinen Fehler, der in den nächsten Schritt wieder eingeht: Die Fehler summieren sich. Man kann den Fehler verringern, in dem man die Schrittweite h verkleinert. Dies erhöht aber den Rechenaufwand und erzeugt auch irgendwann numerische Probleme. Andererseits besticht das Verfahren, weil es anschaulich ist und einfach programmierbar.
Es kann vorkommen, dass sich die Ableitung nicht isolieren lässt, weil sie an mehreren Stellen der Formel vorkommt. In solchen Fäller erfordert die Berechnung der Ableitung selbst ein Näherungsverfahren, z.B. das allgemeine Iterationsverfahren. Es ist jedoch unwichtig, wie aufwendige eine zumindest näherungsweise Berechnung der Ableitung ist: Wenn man sie irgendwie berechnen kann, dann lässt sich dieses Verfahren auch einsetzen.
In den Bewegungsgleichungen der Astronomie trten immer höhere Ableitungen auf. Neben dem Ort eines Körpers ist dessen Geschwindigkeit (die erste Ableitung des Ortes nach der Zeit) und die Beschleunigung in Flger der auf ihn wirkenden Kräfte (die zweite Ableitung des Ortes nach der Zeit) von Bedeutung. Solche Differentialgleichungen können gelöst werden, wenn man die eine Gleichung höherer Ordnung ein System von Gleichungen einfacher Ordnung umformt. Ich erläutere das an dem bereit bekannten Beispiel des Fallschirmspringers:
Das einzige, was wir zunächst wissen,
ist die Kräftebilanz – die Gesamtkraft ist die Summe aus Erdanziehung und
Luftwiderstand
F = FG
+ FW
Die Gravitationskraft kann als
konstant angenommen werden, FG = m · g . Der
Luftwiderstand hängt ab von der Querschnittsfläche, der Geschwindigkeit und
der Luftdichte. Die Querschnittsfläche und
die Luftdichte werden in unserem Fall als konstant angesehen und durch eine
Konstante β ausgedrückt. Die Widerstandkraft ist dann FW
= - β · v.
Wenn Ausgangsbedingungen gegeben sind:
dann kann man aus der Kräftebilanz die Beschleunigung berechnen, also die zweite Ableitung des Ortes (bei uns: der Höhe h) nach der Zeit:
h'' = [ FG(h) + FW(h') ] /m
Die Gewichtskraft ist hier eine Konstante, der Luftwiderstand hängt von der
Geschwindigkeit ab. Gegeben und berechenbar ist also die Beschleunigung, unbekannt sind der
Verlauf der Geschwindigkeit und der Höhe des Fallschirmspringers über die Zeit. Um nun das
Euler-Cauchy-Verfahren (oder ein anderes, genaueres verfahren) einsetzen zu
können, wird eine Hilfsvariable eingeführt, welche die erste
Ableitung der Höhe enthält. In unserem Beispiel hat diese eine physikalische
Bedeutung, es ist nämlich die Geschwindigkeit v.
Damit können wir aus der
einen Gleichung mit einer Ableitung zweiten Grades zwei gekoppelte Gleichungen
machen, welche nur noch Ableitungen ersten Grades enthalten. Des sind
Anstiege, mit denen Polygone gezogen werden können:
v' = FG/m + FW/m
h' = v
Für eine Differentialgleichung zweiter Ordnung benötigen wir auch zwei Anfangsbedingungen, nämlich eine Höhe h0 und eine Geschwindigkeit v0. Beide können wir zum Beispiel auf Null setzen. Das Polygonzugverfahren lässt sich jetzt wieder anwenden. Es wird ein Zeitschritt Δt festgelegt, dann berechnet man für jeden Schritt
Im Python-Programm habe ich zuerst die neue Höhe und dann die Geschwindigkeit berechnet, weil ich dadurch eine Hilfsvariable für die neue Geschwindigkeit eingespart habe. Wichtig für das Euler-Cauchy-Verfahren ist nur das Hauptprogramm. Es hat als Voraussetzung, dass sich die Beschleunigung prinzipiell berechnen lässt.
# Euler-Cauchy-Verfahren
from math import *
# Konstanten
g=9.81
b=15
# Anfangsbedingungen
v0=0
h0=0
m=70
def beschleunigung(v):
a=g
a=a-(b*v)/m
return a
# Hauptprogramm
v=v0
h=h0
dt=1 # Zeitschritt
t=0
while t<20:
print(t, v)
a=beschleunigung(v)
h=h+v*dt
v=v+a*dt
t=t+dt
Die Zeitschrittweite wurde mit 10 ms sehr kurz gewählt. Hiermit ergibt sich aber eine gute Genauigkeit und die vorhergesagte Endgeschwindigkeit von 45 m/s wird auch errechnet. Auch die Zeitkonstante stimmt:
Ihr könnte mit dem Programm etwas experimentieren. Wenn man die Schrittweite zu kurz wählt, dann schwingen die Ergebnisse oder werden gar ganz unsinnig.Für wirkliche Anwendungen sind die Verfahren nach Karl Heun oder nach Carl Runge und Martin Wilhelm Kutta besser geeignet. Ich habe dazu einen eigenen Aufsatz verfasst.
Für eine Differentialgleichung zweiter Ordnung sind zwei
Randbedingungen nötig. Im eben besprochenen Fall waren diese am Anfang des
Rechenintervalls als sogenannte Anfangsbedingungen gegeben. Es ist aber auch
möglich, dass die Bedingungen für das Ende des Rechenintervalls bekannt sind.
Auch das ist kein Problem, dann rechnet man eben in der Zeit rückwärts.
Schwierig wird es, wenn ein Teil der Bedingungen für den Anfang und ein Teil
für das Ende gegeben wurden. Ein Beispiel: Ein Pilot verlässt das Flugzeug mit
einem Schleudersitz, welcher ihn senkrecht nach oben schießt. Nach 5 Sekunden
öffnet der Fallschirm. Zu diesem Zeitpunkt hat der Pilot eine Geschwindigkeit
von 10 m/s. Wie hoch war die Anfangsgeschwindigkeit, die ihm der Schleudersitz
verlieh?
Solche Fälle sind generell schwierig zu berechnen. Geeignet ist das
sogenannte
Schießverfahren. Der Begriff stammt aus dem sog. "Einschießen" der
Artillerie, wo zunächst mit geschätzten Werten geschossen wird, das Ergebnis
beobachtet und danach korrigiert. Man nimmt also irgendwelche Startwerte an,
rechnet 5 Sekunden vorwärts und schaut auf das Ergebnis. Dann korrigiert man
die Startwerte, bis man nahe genug am gewünschten Resultat ist.
Mit den
heutigen schnellen Rechner kann man auch erwägen, für einige Dutzend
Startgeschwindigkeiten ein Ergebnis zu berechnen, und aus der Tabelle der
Ergebnisse das best geeignete herauszusuchen. Das geht aber nur, wenn die
Rechnung in sich nicht allzu aufwendig ist.
Schießverfahren sind meist schlecht zu automatisieren.
(Uwe Pilz, Dezember 2018)