Das Babylonische Wurzelziehen (oft auch Heron-Verfahren) ist ein alter iterativer Algorithmus zur Bestimmung einer rationalen Näherung der Quadratwurzel einer Zahl. Es ist ein Spezialfall des Newton-Verfahrens.

Die Iterationsvorschrift lautet:

$x\_\{n+1\}=\backslash frac\{x\_n\; +\; \backslash frac\{a\}\{x\_n\}\}\{2\}$.

Hierbei steht $a$ für die Zahl, deren Quadratwurzel bestimmt werden soll. Der Startwert $x\_0$ der Iteration kann, solange er nicht gleich Null ist, beliebig festgesetzt werden, wobei zu beachten ist, dass negative Werte gegen die negative Quadratwurzel konvergieren.

Triviales Beispiel

---

Im Folgenden ein triviales Beispiel für die Wurzel aus 9 und die Annäherung nach vier Berechnungsschritten an den wahren Wert $\backslash sqrt\{9\}\; =\; 3$:

$a=9$ und $x\_0=1$

$x\_1=\backslash frac\{1\; +\; \backslash frac\{9\}\{1\}\}\{2\}\; =\backslash frac\{10\}\{2\}=5$

$x\_2=\backslash frac\{5\; +\; \backslash frac\{9\}\{5\}\}\{2\}\; =\backslash frac\; \{\backslash frac\; \{34\}\{5\}\}\; \{2\}\; =\backslash frac\{34\}\{10\}\; =\; 3\{,\}4$

$x\_3=\backslash frac\{\backslash frac\{34\}\{10\}\; +\; \backslash frac\{9\}\{\backslash frac\{34\}\{10\}\}\}\{2\}\; =\backslash frac\{\backslash frac\{34\}\{10\}\; +\; \backslash frac\{90\}\{34\}\}\{2\}\; =\backslash frac\{257\}\{85\}\; \backslash approx3\{,\}0235$

$x\_4=\backslash frac\{\backslash frac\{257\}\{85\}\; +\; \backslash frac\{9\}\{\backslash frac\{257\}\{85\}\}\}\{2\}\; =\backslash frac\{\backslash frac\{257\}\{85\}\; +\; \backslash frac\{765\}\{257\}\}\{2\}\; =\backslash frac\{65537\}\{21845\}\; \backslash approx3\{,\}00009$

Geometrische Veranschaulichung des Heron-Verfahrens

---

Der Flächeninhalt eines Quadrates kann über das Quadrat der Länge seiner Seiten berechnet werden. Die Bestimmung der Quadratwurzel einer gegebenen Zahl $a$ kann also geometrisch gedeutet werden als Bestimmung der Seitenlänge (genauer: als rationale Näherung zu $\backslash sqrt\{a\}$) eines Quadrates mit dem Flächeninhalt $a$.

Die Idee ist nun, von einem Rechteck des Flächeninhaltes $a$ auszugehen und die Seitenlängen einem Quadrat immer weiter anzunähern.

Dazu wird ein Startwert gewählt, im obigen Fall gilt $a=9$ und als Startwert wurde $x\_0=1$ gewählt. Geometrisch bedeutet dieses, dass von einem Rechteck der Seitenlänge $x\_0=1$ ausgegangen wird. Die andere Seitenlänge dieses Rechtecks ergibt sich aus dem vorgegebenen Flächeninhalt: $y\_0=\backslash frac\{9\}\{1\}=\; 9$

Bei der Betrachtung ist unmittelbar ersichtlich, dass es eine geeignetere Näherung an ein Quadrat gibt, denn die eine Seitenlänge $x\_0=1$ ist zu klein, die andere mit $y\_0=9$ zu groß.

Um eine verbesserte Annäherung an die Länge einer Quadratseite zu erhalten, kann das arithmetische Mittel der Seitenlängen $x\_0$ und $y\_0$ dienen (Hier gibt es eine ganze Reihe von Möglichkeiten, das Verfahren zu verfeinern.):

$x\_1=\backslash frac\{x\_0\; +\; y\_0\}\{2\}\; =\backslash frac\{9\; +\; 1\}\{2\}=5$

Die Länge der zweiten Seite dieses neuen Näherungs-Rechtecks ergibt sich wieder durch den vorgegebenen Flächeninhalt $a$:

$y\_1=\backslash frac\{a\}\{x\_1\}=\backslash frac\{9\}\{5\}=\; 1\{,\}8$

Die Werte $x\_1=5$ und $y\_1=1\{,\}8$ sind geometrisch gedeutet die Seitenlängen eines zweiten Näherungs-Rechtecks. Dieses und die folgenden Rechtecke lassen sich nun weiter verbessern durch erneute Bildung des arithmetischen Mittelwertes als verbesserte Näherung an die Seitenlänge eines Quadrates mit der Seitenlänge $\backslash sqrt\{a\}$.

Konvergenz

---

Das Verfahren konvergiert relativ rasch innerhalb weniger Schritte. Da es sich aus dem Newtonschen Näherungsverfahren ableiten lässt, ist die Konvergenzordnung 2.
Es gelten:
$\backslash sqrt\{a\}\; \backslash le\; x\_\{n+1\}\; \backslash le\; x\_n\; (n\; =\; 1,\; 2...)$
und
$\backslash lim\_\{n\; \backslash to\; \backslash infty\}x\_n\; =\; \backslash sqrt\{a\}$

Fehler

---

Für den Fehler der Heron-Folge $\{x\_n\}\; (n\; \backslash ge\; 1)$ gilt:
$\backslash frac\{a\}\{x\_n\}\; \backslash le\; \backslash sqrt\{a\}\; \backslash le\; x\_n$ (Einschließung), sowie
$x\_n-\backslash sqrt\{a\}\; \backslash le\; \backslash frac\; \{1\}\{2\backslash sqrt\{a\}\}\; \backslash left(\; x\_\{n-1\}-\backslash sqrt\{a\}\; \backslash right)^2$ (quadratische Konvergenz)

Implementierung in Software

---

Das Verfahren eignet sich besonders gut zur Implementierung in Software, wird heute angesichts der breiten Verfügbarkeit numerischer Prozessorhardware aber nur noch selten benötigt. Man kommt dabei nur mit Grundrechenarten aus, s. o.

Wenn dazu noch eine Gleitkommadarstellung mit einem Zweier-Exponenten benutzt wird, wird der Ansatz relativ einfach:

• Zunächst spaltet man vom Exponenten eine gerade Anzahl ab, so dass als Exponent entweder eine 0 oder 1 übrigbleibt, die Zahl also auf das Intervall { ½ , 2 } normalisiert wird. In diesem Intervall ist die Wurzelfunktion eine nur schwach gekrümmte Kurve, lässt sich also numerisch gut behandeln.
• Als Startwert für die eigentliche Iteration gibt es mehrere Möglichkeiten, diese Kurve durch eine noch einfachere zu approximieren. Mit den komplizierteren Ansätzen lässt sich ggf. ein Iterationsschritt einsparen:
  • eine einfache Konstante (z. B. 1);
  • eine Gerade mit Steigung 1/2 und einer additiven Konstante von 1 (als Vereinfachung des nachfolgenden Falls);
  • eine Gerade mit Steigung 1/2 und einer additiven, optimierten Konstante von $\backslash left(\; 2\backslash cdot\; \backslash sqrt*\{2\}\; -\; \backslash sqrt\{2\}\; \backslash right)^2\; \backslash approx\; 0\{,\}929683$ ;
  • eine Gerade mit optimierter Steigung und einer additiven Konstante.
• Ausgehend von dem so ermittelten Startwert führt man eine feste Anzahl von Iterationsschritten durch. Die nötige Anzahl, um die gewünschte Genauigkeit zu erreichen, lässt sich dank der obigen Fehlerabschätzung als Worst Case innerhalb des Startintervalls direkt ausrechnen. Bei 32 Bits Mantisse und dem mittleren Startansatz braucht man z. B. nur drei Schritte. Diese fest gewählte Anzahl erspart wesentlich aufwendigere Abfragen auf Erreichung der Genauigkeit. Der Ersatz der genannten Konstanten durch die Zahl 1,0 ändert daran nichts. Auch der noch kompliziertere Ansatz brächte zumindest bei dieser Genauigkeit keine Einsparung eines weiteren Iterationsschritts. Bei höheren Genauigkeitsanforderungen kann das anders aussehen.
• Zum Schluss muss man den Exponenten restaurieren, indem man die Hälfte des im ersten Schritt abgespalteten Werts wieder hinzufügt.

Verallgemeinerung des Verfahrens

---

Dieses Verfahren kann man leicht verallgemeinern, so dass man die n-te Wurzel berechnen kann. Je größer n aber ist, umso mehr Schritte werden benötigt, um die Wurzel genau zu berechnen.

Die Iterationsvorschrift lautet hier:

$x\_\{m+1\}=\backslash frac\{x\_m\; +\; \backslash frac\{a\}\{x\_m^\{n-1\}\}\}\{2\}$

Geschichte

---

Dieses Verfahren ist auch unter dem Namen Heron-Verfahren bekannt. Es wurde nach Heron von Alexandria benannt und entstammt seiner Formelsammlung.

Numerische Mathematik