Heron-Verfahren

Grundidee und Bedeutung

Das Heron-Verfahren, auch Heronsches Näherungsverfahren oder Babylonisches Wurzelziehen genannt, berechnet näherungsweise die Quadratwurzel einer reellen Zahl a > 0. Es benötigt keinen bereits korrekten Ausgangswert, ist relativ robust gegen Rundungsfehler und konvergiert meist schneller als das schriftliche Wurzelziehen. Allerdings liefert es grundsätzlich Näherungswerte.

Die geometrische Idee geht von einem Rechteck mit Flächeninhalt A aus. Ein Quadrat mit demselben Flächeninhalt hätte die Seitenlänge √A. Deshalb verändert man das Seitenverhältnis des Rechtecks schrittweise, bis seine Seiten nahezu gleich lang sind, während sein Flächeninhalt unverändert bleibt.

Zunächst wählt man eine beliebige Seitenlänge x₀. Die zweite Seitenlänge ist dann y₀ = A/x₀. Anschließend ersetzt man die längere Seite durch den Mittelwert der bisherigen Seiten:

x₁ = (x₀ + y₀)/2.

Die andere Seite wird so angepasst, dass der Flächeninhalt A erhalten bleibt:

y₁ = A/x₁.

Dieser Ablauf wird wiederholt. Für A = 9 kann man mit den Seitenlängen 9 und 1 beginnen. Der erste Mittelwert ist 5; die andere Seite beträgt dann 9/5 = 1,8. Nach weiteren Schritten wird das Rechteck immer quadratischer. Eine entstehende Seitenlänge von 3,024 liegt bereits nahe am exakten Wert √9 = 3.

Das Verfahren wurde bereits zur Zeit des babylonischen Königs Hammurapi I. um 1750 v. Chr. in Mesopotamien verwendet. Heron von Alexandria beschrieb es um 100 n. Chr. im ersten Buch seines Werkes Metrika.

Iterationsvorschrift und Beispiel

Aus der geometrischen Idee entsteht ein Iterationsverfahren, also eine wiederholte Anwendung derselben Rechenvorschrift. Für a > 0 wählt man einen Startwert x₀ ≠ 0, möglichst in der Nähe von √a, und setzt y₀ = a/x₀. Danach gilt in jedem Schritt:

xₙ₊₁ = (xₙ + yₙ)/2 und yₙ₊₁ = a/xₙ₊₁.

Setzt man yₙ = a/xₙ direkt ein, erhält man die übliche Rekursionsgleichung

xₙ₊₁ = 1/2 · (xₙ + a/xₙ).

Die dadurch entstehende Folge x₀, x₁, x₂, … heißt Heron-Folge. Ihre Glieder sind zunehmend genaue Näherungen für √a. Dieselbe Vorschrift ergibt sich auch aus dem Newton-Verfahren, wenn man die Nullstelle der Funktion f(x) = x² − a sucht.

Beispiel für √2 mit x₀ = 2:

• x₁ = 1/2 · (2 + 2/2) = 3/2 = 1,5 • x₂ = 1/2 · (3/2 + 2/(3/2)) = 17/12 = 1,41666666… • x₃ = 1/2 · (17/12 + 2/(17/12)) = 577/408 = 1,41421568…

Der wahre Wert ist √2 = 1,41421356…. Bereits nach drei Iterationen ist die Näherung auf fünf Nachkommastellen genau; die Abweichung beträgt weniger als 0,0001 %.

Konvergenz und Fehler

Für jeden positiven Startwert x₀ konvergiert die Heron-Folge gegen die positive Quadratwurzel √a. Alle Folgenglieder sind positiv. Ab dem ersten Iterationsschritt gilt außerdem xₙ ≥ √a, denn

xₙ² − a = 1/4 · (xₙ₋₁ − a/xₙ₋₁)² ≥ 0.

Zugleich ist die Folge monoton fallend, weil

xₙ₊₁ − xₙ = (a − xₙ²)/(2xₙ) ≤ 0.

Sie ist somit nach unten beschränkt und monoton und besitzt deshalb einen Grenzwert x. Aus der Iterationsgleichung folgt für diesen Grenzwert

x = 1/2 · (x + a/x),

also x² = a und wegen x > 0 schließlich x = √a. Auch jeder negative Startwert konvergiert, sofern er nicht null ist; in diesem Fall geht die Folge gegen die negative Quadratwurzel.

Das Verfahren besitzt Konvergenzordnung 2, also quadratische Konvergenz. In der Nähe des gesuchten Wertes verdoppelt sich mit jedem Schritt ungefähr die Zahl der richtigen Stellen. Ein schlechter Startwert kann die Berechnung jedoch zunächst stark verlängern. Soll beispielsweise die Wurzel einer ganzen Zahl a mit 200 Binärstellen berechnet werden und wird x₀ = a gewählt, benötigt die Näherung ungefähr 100 Schritte, bis sie die richtige Länge von 100 Stellen erreicht. Danach genügen etwa sechs bis sieben weitere Schritte, entsprechend log₂(100), um alle 100 Stellen vor dem Komma richtig zu bestimmen. Zweckmäßig ist daher ein Startwert mit ungefähr der halben Bitlänge von a. Eine Schätzung aus den ersten beiden Gliedern der Taylorentwicklung der binomischen Reihe um 1 ist x₀ = (a + 1)/2.

Das Heron-Verfahren ist außerdem ein Fixpunktverfahren: Für φ(x) = 1/2 · (x + a/x) erfüllt ein Fixpunkt φ(x) = x und damit x² = a.

Für n ≥ 1 schließen die beiden Rechteckseiten die gesuchte Wurzel ein:

a/xₙ ≤ √a ≤ xₙ.

Der absolute Fehler lässt sich abschätzen durch

xₙ − √a = 1/(2xₙ₋₁) · (xₙ₋₁ − √a)² ≤ 1/(2xₙ₋₁xₙ²) · (xₙ₋₁xₙ − a)².

Im Beispiel für √2 ergibt dies bei x₃ eine obere Fehlerschranke von 0,00000212…. Für den relativen Fehler εₙ = (xₙ − √a)/√a gilt

εₙ₊₁ = εₙ²/[2(1 + εₙ)].

Bei einem festgelegten relativen Anfangsfehler ε₀ ist die Folge der relativen Fehler daher unabhängig von a.

Umsetzung in Software

Das Verfahren lässt sich leicht programmieren, weil es nur Grundrechenarten benötigt. Wegen moderner numerischer Prozessorhardware wird es heute jedoch nur noch selten ausdrücklich benötigt.

Bei einer Gleitkommadarstellung mit Zweier-Exponent kann die Zahl zunächst auf das Intervall [1/2, 2] normalisiert werden. Dazu spaltet man vom Exponenten eine gerade Anzahl ab, sodass als Rest 0 oder 1 bleibt. In diesem Intervall ist die Wurzelfunktion nur schwach gekrümmt und numerisch gut zu behandeln. Für √5 gilt beispielsweise

√5 = √(4 · 1,25) = 2 · √1,25 ≈ 2 · 1,118034 = 2,236068.

Die Iteration muss daher zunächst nur √1,25 bestimmen. Als Startwert kann man etwa die Konstante 1 verwenden; ihr relativer Fehler beträgt hier 1,1 · 10⁻¹. Eine genauere einfache Näherung ist die Gerade x₀ = 1/2 + a/2. Für a = 1,25 liefert sie x₀ = 1,125 mit einem relativen Fehler von 6,2 · 10⁻³. Möglich ist auch eine Gerade mit Steigung 1/2 und der optimierten additiven Konstante (2 · ⁴√2 − √2)²/2 ≈ 0,4648415.

Bei einer Mantisse von 32 Bits reichen mit dem mittleren Startansatz beispielsweise drei Iterationsschritte. Eine vorher festgelegte Schrittzahl vermeidet zusätzliche Prüfungen, ob die Zielgenauigkeit schon erreicht wurde. Ausgehend von x₀ = 1 erhält man für a = 1,25:

• x₁ = 1,125, relativer Fehler 6,2 · 10⁻³ • x₂ ≈ 1,118056, relativer Fehler 2,0 · 10⁻⁵ • x₃ ≈ 1,118034, relativer Fehler kleiner als 10⁻⁶.

Hier zeigt sich die quadratische Konvergenz: Der relative Fehler wird von Schritt zu Schritt ungefähr quadriert. Abschließend wird der zuvor abgespaltene Exponent zur Hälfte wieder eingesetzt; im Beispiel ergibt sich 2 · x₃ ≈ 2,236068.

Höhere Wurzeln und Kehrwerte

Das Verfahren kann mithilfe des Newton-Verfahrens auf die k-te Wurzel von a > 0 erweitert werden. Für f(x) = xᵏ − a und f′(x) = kxᵏ⁻¹ lautet die Iterationsvorschrift

xₙ₊₁ = 1/k · ((k − 1)xₙ + a/xₙᵏ⁻¹).

Für die Kubikwurzel ergibt sich beispielsweise

xₙ₊₁ = 1/3 · (2xₙ + a/xₙ²).

Die Folge muss mit einem geeigneten Startwert für ⁿ√a begonnen werden. Mit wachsendem k werden mehr Schritte benötigt. Für jedes positive ganzzahlige k gelten dieselben Konvergenzaussagen wie für k = 2.

Für k = −1 entsteht ein quadratisch konvergierendes Verfahren zur Näherung des Kehrwerts 1/a, das keine Division benötigt:

xₙ₊₁ = 2xₙ − axₙ² = (2 − axₙ)xₙ.

Es konvergiert für alle Startwerte x₀ im offenen Intervall (0, 2/a) quadratisch gegen 1/a. Bei frühen Computern ohne eingebaute Division konnte man so eine Division auf Multiplikation und Subtraktion zurückführen: Zuerst wurde der Kehrwert des Nenners angenähert und anschließend mit dem Zähler multipliziert.

Für 1/3 und x₀ = 1/2 erhält man:

• x₁ = (2 − 3 · 1/2) · 1/2 = 1/4 = 0,25 • x₂ = (2 − 3 · 1/4) · 1/4 = 5/16 = 0,3125 • x₃ = (2 − 3 · 5/16) · 5/16 = 85/256 = 0,33203125.

Die Bedingung an den Startwert ist wesentlich: Für x₀ = 2/3, also den nicht zum offenen Intervall gehörenden Randwert, ist x₁ = 0. Danach bleiben alle Folgenglieder null und konvergieren nicht gegen 1/3.