Quadratwurzeln mit dem Halley-Verfahren berechnen

Das Halley-Verfahren (nach dem Astronomen Edmond Halley) kann, wie das Newton-Verfahren, verwendet werden, um eine Nullstelle einer Funktion zu bestimmen. Es hat die Iterationsvorschrift

x+ = x - (2 f(x) f'(x)) / (2 (f'(x))² - f(x) f''(x))

um eine Näherung x zu einer Näherung x+ zu verbessern.

Für die Anwendung zum Bestimmen der Quadratwurzel verwendet man als Funktion f(x) = x² - a mit den beiden Ableitungen f'(x) = 2x und f''(x) = 2. Somit ist die Iterationsvorschrift

x+ = x - (2 (x² - a) 2x) / (2 (2x)² - (x² - a) 2)
   = x - 2 (x³ - ax) / (3 x² + a)

in diesem Spezialfall.

Die Realisierung in Lisp drückt diese Iterationsvorschrift durch die Hilfsfunktion halley-square-root-next aus.

(defproc halley-square-root-next (a x)
  (- x
    (/ (* 2 (- (* x x x) (* a x)))
       (+ (* 3 x x) a))))

Die Iteration läuft solange, bis die absolute Differenz des Quadrats der Wurzel und der ursprünglichen Zahl eine Grenze unterschreitet.

(defproc halley-square-root-iteration (a x delta)
  (let
    ((x-next (approximate (halley-square-root-next a x) delta)))
    (if
      (< (abs (- (* x-next x-next) a)) delta)
      x-next
      (halley-square-root-iteration a x-next delta))))

Die Funktion square-root prüft ihr erstes Argument und steigt dann in die oben beschriebene Iteration ein. Das Ergebnis der Iteration wird auf die gewünschte Genauigkeit gerundet.

(defproc square-root (a delta)
  (if
    (< a 0)
    (throw (quote error) "negative argument for square-root")
    (approximate
      (halley-square-root-iteration a 1 (* delta delta))
      delta)))

Quellen
https://de.wikipedia.org/wiki/Halley-Verfahren
http://mathworld.wolfram.com/HalleysMethod.html


square-root.sheet.xml