Halley-Verfahren

Grundidee und Standarditeration

Das Halley-Verfahren, auch Verfahren der berührenden Hyperbeln, ist ein Verfahren der numerischen Mathematik zur Bestimmung von Nullstellen reeller Funktionen f: ℝ → ℝ, also von Lösungen der Gleichung f(x) = 0. Es ähnelt dem Newton-Verfahren, verwendet aber zusätzlich zur ersten Ableitung f′(x) auch die zweite Ableitung f″(x). Dadurch besitzt es unter geeigneten Voraussetzungen die Konvergenzordnung 3: Die Zahl der korrekten Stellen wird in der Nähe einer einfachen Nullstelle ungefähr verdreifacht.

Ist f ∈ C²(ℝ, ℝ) und a eine einfache Nullstelle mit f′(a) ≠ 0, so konvergiert die durch

xₖ₊₁ = xₖ − 2 f(xₖ) f′(xₖ) / (2 f′(xₖ)² − f(xₖ) f″(xₖ)), k = 0, 1, 2, …

erzeugte Folge für Startwerte x₀ nahe a gegen a. Der Nenner muss dabei von null verschieden sein. Pro Schritt sind mehr Rechenoperationen nötig als beim Newton-Verfahren, dafür ist die lokale Konvergenz deutlich schneller.

Motivation und Herleitung

In der Nähe einer einfachen Nullstelle wird der Funktionsverlauf von f in zweiter Ordnung „gerade gebogen“. Dazu betrachtet man statt f die Funktion

g(x) = f(x) / √|f′(x)|.

Diese Konstruktion hängt nicht von der gesuchten Nullstelle ab. Wendet man das Newton-Verfahren auf g an, ergibt sich mit

g′(x) = (2 f′(x)² − f(x) f″(x)) / (2 f′(x) √|f′(x)|)

die Iteration

H_f(x) = N_g(x) = x − 2 f(x) f′(x) / (2 f′(x)² − f(x) f″(x)).

Dieselbe Vorschrift folgt aus dem Householder-Verfahren zweiter Ordnung:

H_f(x) = x + 2 (1/f)′(x) / (1/f)″(x).

Das Verfahren kann daher als Newton-Verfahren auf eine geeignet veränderte Funktion verstanden werden, deren Verlauf in der Umgebung der Nullstelle genauer angenähert wird.

Parabolische und hyperbolische Varianten

Das ursprünglich von E. Halley vorgeschlagene irrationale oder parabolische Halley-Verfahren ersetzt den Graphen von f in der Umgebung des aktuellen Werts x durch die Taylor-Parabel zweiten Grades

f(x + h) = f(x) + f′(x)h + ½ f″(x)h² + O(h³).

Die daraus bestimmte Nullstelle liefert die Iteration

xₖ₊₁ = xₖ − 2 f(xₖ) / [f′(xₖ) + sign(f′(xₖ)) √(f′(xₖ)² − 2 f(xₖ) f″(xₖ))].

Durch die Entwicklung √(1 + y) = 1 + ½y + O(y²) erhält man daraus das heute übliche rationale oder hyperbolische Halley-Verfahren, also die Standardformel. Beim hyperbolischen Verfahren wird f durch eine Hyperbel angenähert, die den Graphen in x ebenfalls bis zur zweiten Ordnung berührt. Der Nullpunkt des Zählers dieser Hyperbel führt zu

h = −2 f(x) f′(x) / (2 f′(x)² − f(x) f″(x)).

Eine weitere Verallgemeinerung ist das Laguerre-Verfahren:

xₖ₊₁ = xₖ − n f(xₖ) / f′(xₖ) · [1 + (n − 1) √(1 − n/(n − 1) · f(xₖ) f″(xₖ) / f′(xₖ)²)]⁻¹.

Für Polynome wird n gleich dem Grad des Polynoms gesetzt. Da der Ausdruck unter der Quadratwurzel negativ werden kann, können das irrationale Halley-Verfahren und das Laguerre-Verfahren auch bei reellen Polynomen und reellen Startwerten zu komplexen Nullstellen konvergieren. Bei Quadratwurzeln aus komplexen Zahlen ist die Lösung mit positivem Realteil zu wählen, damit der Nenner den größtmöglichen Betrag besitzt.

Warum die Konvergenz kubisch ist

Sei f dreimal stetig differenzierbar und a eine Nullstelle. Auf einem Intervall I, das a enthält, gilt nach dem Mittelwertsatz

|f(x)| / maxᵧ∈I |f′(y)| ≤ |x − a| ≤ |f(x)| / minᵧ∈I |f′(y)|.

Daher sind x − a = O(f(x)) und f(x) = O(x − a) äquivalent. Beim parabolischen Verfahren ist der Schritt h = O(f(x)); durch die Wahl von h verschwinden die ersten drei Glieder der Taylorentwicklung, sodass f(xₖ₊₁) = O(h³) = O(f(xₖ)³) gilt.

Beim hyperbolischen Verfahren wird die Taylorentwicklung mithilfe von (a + bh)(a − bh) = a² − b²h² = a² + O(h²) in einen Bruch linearer Funktionen umgeformt. Dadurch entsteht eine Hyperbelapproximation zweiter Ordnung. Für den Halley-Schritt gilt ebenfalls h = O(f(x)) und

f(H_f(x)) = O(h³) = O(f(x)³).

Zusammen mit den Abschätzungen zwischen Funktionswert und Abstand zur Nullstelle folgt

xₖ₊₁ − a = O((xₖ − a)³).

Das ist die kubische Konvergenzordnung.

Beispiel und mehrdimensionale Erweiterung

Für die Quadratwurzel von a = 5 verwendet man f(x) = x² − a. Die Halley-Iteration lautet

H_f(x) = x (x² + 3a) / (3x² + a) = x (x² + 15) / (3x² + 5).

Mit x₀ = 3 ergeben sich unter anderem folgende Werte:

k = 0: x₀ = 3, f(x₀) = 4 k = 1: x₁ = 2,25, f(x₁) = 0,0625 k = 2: x₂ = 2,23606811145510835913312693498452012383900928792569659442724, f(x₂) = 5,99066414899E−7 k = 3: x₃ = 2,23606797749978969640929385361588622700967141237081284965284, f(x₃) = 5,37483143712E−22 k = 4: x₄ = 2,23606797749978969640917366873127623544061835961152572427090, f(x₄) = 0,000000000000.

Die Folge besitzt 0, 1, 5, 21, >60 gültige Stellen; dies entspricht ungefähr einer Verdreifachung in jedem Schritt. Das Newton-Verfahren verwendet hier G_f(x) = (x² + a)/(2x) = (x² + 5)/(2x) und konvergiert im direkten Vergleich langsamer.

Für F: ℝⁿ → ℝⁿ lässt sich das Verfahren erweitern. Dabei ist F′(x) eine als invertierbar vorausgesetzte Matrix, während F″(x) ein Tensor dritter Stufe beziehungsweise eine vektorwertige symmetrische Bilinearform ist. Außerdem kommutiert F″(x)(v, ·) im Allgemeinen nicht mit F′(x).

Der mehrdimensionale Schritt kann in drei Teilaufgaben berechnet werden:

• Newton-Schritt: Löse F′(xₖ)sₖ = −F(xₖ).
• Korrektur: Löse [F′(xₖ) + ½F″(xₖ)(sₖ, ·)]tₖ = −F(xₖ).
• Aktualisierung: Setze xₖ₊₁ = xₖ + tₖ.

Ist die zweite Ableitung Lipschitz-stetig, konvergiert dieses Verfahren lokal kubisch. Eine Näherung ohne die große Inverse führt zum Euler-Tschebyschow-Verfahren:

t = −F′(x)⁻¹F(x) − ½F′(x)⁻¹F″(x)(F′(x)⁻¹F(x), F′(x)⁻¹F(x)).