Intervall-Newton-Verfahren

Einige Sternfreunde befassen sich mit der Himmelsmechanik. Ein zentrales Problem dabei ist die Lösung der Keplergleichung

0 = E - 180°/pi*e*sin(E) - M.

Das Problem mit ihr ist, daß man sie nicht algebraisch nach E auflösen kann - E kann nur numerisch bestimmt werden, z.B. mit dem Newton-Verfahren:

En+1 = En - f(En)/f '(En),

wobei f(E) die rechte Seite der Keplergleichung ist und f '(E) ihre erste Ableitung.

Ausgehend von einem Startwert für E, den man in die Vorschrift einsetzt, erhält man einen Näherungswert zum erneuten Einsetzen in die Vorschrift, usw. Dies wird fortgesetzt, bis man einen Wert bekommen hat, der genau genug ist, denn mit jedem Iterationsschritt nähert sich E der realen Lösung an. Das Problem mit dem Newton-Verfahren ist, daß es keine Gewähr für diese Annäherung (Konvergenz) gibt. Ob dieses Verfahren funktioniert, hängt von der Wahl der Startwerte ab.

Doch welchen Startwert soll man nehmen? Häufig hat man für ihn keinen Anhaltspunkt. Gibt es bessere Verfahren, mit denen sich das Dilemma vermeiden läßt?

Einen Ausweg bietet das Intervall-Newton-Verfahren auf Basis der Intervallarithmetik an. Im Gegensatz zur konventionellen Arithmetik, bei der mit diskreten Zahlenwerten gerechnet wird, betrachtet die Intervallarithmetik das Rechnen mit kompletten Intervallen. Die Grundrechenarten, wie wir sie kennen, gibt es auch bei der Intervallarithmetik (die Addition zweier Intervalle wird so definiert, daß die Summe zweier Zahlen aus den beiden Intervallen im Ergebnisintervall garantiert enthalten ist: [a,b]+[c,d] = [a+c,b+d]). Die übrigen Grundrechenarten werden mit der selben Anforderung definiert. Eine Besonderheit stellt die Division dar: ihr Ergebnis können zwei disjunkte Intervalle sein.

Mit Hilfe der Intervallarithmetik kann man das Newton-Verfahren neu implementieren. Das Intervall-Newton-Verfahren hat nämlich zwei sehr schöne Nebeneffekte: nicht nur, daß die Konvergenz gegen die Lösungen unabhängig von den Startwerten sichergestellt ist (vorausgesetzt, es gibt auf dem Startintervall mindestens eine Lösung), sondern man bekommt zudem alle Lösungen, wenn die Gleichung mehrere besitzt! In diesem Falle zerfällt das Intervall bei irgendeinem Iterationsschritt in mehrere Teilintervalle (durch die Division), von denen jedes wiederum eine Lösung enthält (das Abbruchkriterium ist die abnehmende Intervallbreite - der Abstand der Nullstellen muß darüber liegen). Gibt es hingegen keine Lösung, so bekommt man ein Leerintervall und kann das Iterieren abbrechen.

Mir ist keine Programmiersprache bekannt, die von Haus aus Intervallarithmetik unterstützt. Es gibt aber mit Boost eine C++-Bibliothek, die mit verhältnismäßig geringem Aufwand für eigene Anwendungen benutzt werden kann und sogar ein passendes Beispiel enthält.