Bairstowverfahren

Zweck und Grundidee

Das Bairstow-Verfahren ist ein Iterationsverfahren der numerischen Mathematik zur Bestimmung der Nullstellen eines Polynoms mit reellen Koeffizienten. Es sucht nicht direkt nach einzelnen Nullstellen, sondern nähert quadratische Faktoren mit reellen Koeffizienten an. Dadurch können auch komplexe Nullstellen berechnet werden, ohne während der Iteration komplexe Arithmetik zu verwenden.

Nach dem Fundamentalsatz der Algebra kann jedes reelle Polynom in lineare und quadratische irreduzible Faktoren mit reellen Koeffizienten zerlegt werden. Reelle Nullstellen ergeben lineare Faktoren. Nichtreelle komplexe Nullstellen treten bei reellen Polynomen stets als konjugiertes Paar z = u + iv und z̄ = u − iv auf. Zu diesem Paar gehört der reelle quadratische Faktor

a(x) = (x − z)(x − z̄) = (x − u)² + v² = x² − 2ux + (u² + v²).

Das Verfahren bestimmt einen solchen Faktor a(x) = x² + a₁x + a₀. Seine beiden Nullstellen erhält man anschließend mit der Lösungsformel für quadratische Gleichungen. Der gefundene Faktor wird vom Ausgangspolynom abgespalten; danach wird das Verfahren auf das verbleibende Polynom kleineren Grades angewendet.

Das Bairstow-Verfahren ist aus dem Newton-Verfahren abgeleitet und konvergiert lokal mit 2. Ordnung. Eine Iteration mit reellen Zahlen ist erheblich schneller als eine Iteration mit komplexer Arithmetik. Die Konvergenz hängt jedoch stark von einem geeigneten Startfaktor ab. Bei mehrfachen oder dicht beieinanderliegenden Nullstellen kann das Verfahren scheitern.

Polynomdivision und mathematische Grundlage

Gegeben ist das Polynom n-ten Grades

f(x) = fₙxⁿ + fₙ₋₁xⁿ⁻¹ + … + f₁x + f₀

mit ausschließlich reellen Koeffizienten. Für einen versuchsweise gewählten quadratischen Faktor a(x) = x² + a₁x + a₀ wird eine Polynomdivision ausgeführt:

f(x) = a(x)b(x) + r(x).

Dabei hat b(x) den Grad n − 2. Weil a(x) quadratisch ist, besitzt der Rest höchstens den Grad 1 und hat somit die Form r(x) = r₁x + r₀. Genau dann, wenn r₀ = r₁ = 0 gilt, ist a(x) ein exakter Faktor von f(x).

Um einen nur angenäherten Faktor zu verbessern, wird ein lineares Korrekturpolynom Δa(x) gesucht. Unter Vernachlässigung von Termen höherer Ordnung ergibt die Newton-Linearisierung

r(x) = a(x)·Δb(x) + Δa(x)·b(x).

Eine weitere Division von b(x) durch a(x) liefert

b(x) = q(x)a(x) + p(x),

wobei p(x) = p₁x + p₀ wieder ein linearer Rest ist. Durch Rechnung modulo a(x) reduziert sich das Problem auf ein lineares Gleichungssystem mit den beiden Unbekannten Δa₀ und Δa₁:

(r₀, r₁)ᵀ = ((p₀, −p₁a₀), (p₁, p₀ − p₁a₁))·(Δa₀, Δa₁)ᵀ.

Die Determinante der Systemmatrix lautet

D = p₀² + p₁(a₀p₁ − a₁p₀).

Ist D ≠ 0, kann das System gelöst und a(x) korrigiert werden. Danach beginnen die beiden Polynomdivisionen erneut. Abgebrochen wird, sobald beide Restkoeffizienten eine vorher festgelegte Schranke unterschreiten.

Iterationsschema und Startwerte

Ein Iterationsschritt besteht aus folgenden Teilen:

• Division des Polynoms f(x) durch a(x) = x² + a₁x + a₀; sie liefert b(x) und den Rest r(x).

• Erneute Division von b(x) durch a(x); sie liefert q(x) und den Rest p(x).

• Berechnung der Korrekturen für a₀ und a₁ aus den beiden Resten.

Für die erste Division setzt man bₙ = bₙ₋₁ = 0 und berechnet rekursiv

bⱼ = fⱼ₊₂ − a₁bⱼ₊₁ − a₀bⱼ₊₂ für j = n − 2, n − 3, …, 0, −1, −2.

Dann gelten r₁ = b₋₁ und r₀ = b₋₂ + a₁b₋₁. Bei der zweiten Division wird dieselbe Art von Rekursion verwendet; für ihren Rest gelten p₁ = q₋₁ und p₀ = q₋₂ + a₁q₋₁.

Mit

M = −a₀q₋₁ − a₁q₋₂ und D = q₋₂² − Mq₋₁

werden die Korrekturen berechnet:

Δa₁ = (q₋₁b₋₂ − q₋₂b₋₁)/D,

Δa₀ = (Mb₋₁ − q₋₂b₋₂)/D.

Danach setzt man

a₁(neu) = a₁(alt) − Δa₁ und a₀(neu) = a₀(alt) − Δa₀.

Eine mögliche Wahl der Startwerte verwendet die drei führenden Koeffizienten von f:

a₁ = fₙ₋₁/fₙ und a₀ = fₙ₋₂/fₙ.

Sind bereits Näherungen x̃₁ und x̃₂ für zwei Nullstellen bekannt, können nach den Vietaschen Formeln auch

a₁ = −x̃₁ − x̃₂ und a₀ = x̃₁x̃₂

gewählt werden. Das Verfahren endet, wenn r₀ = r₁ = 0 beziehungsweise b₋₁ = b₋₂ = 0 innerhalb der gewünschten Genauigkeit gilt. Dann ist a(x) ein Teiler von f(x). Nach dem Abspalten setzt man f(x) = b(x) und sucht die übrigen Nullstellen. Die Rekursionen lassen sich schnell und speichergünstig verschachteln, sodass nicht alle Koeffizienten von b(x) gespeichert werden müssen.

Rechenbeispiel

Gesucht werden die Nullstellen von

f(x) = 6x⁵ + 11x⁴ − 33x³ − 33x² + 11x + 6.

Aus den führenden Koeffizienten folgen die Startwerte

a₁ = 11/6 und a₀ = −33/6.

Die Näherungen verändern sich während der Iteration zunächst deutlich. Nach der 5. Iteration liegen sie bei a₁ = 3,326244386565 und a₀ = 0,978742927192; die angegebene Schrittweite beträgt 0,022431883898. Nach der 7. Iteration sind a₁ = 3,333333333340 und a₀ = 1,000000000020 bei einer Schrittweite von 0,000000000021 erreicht. Nach der 8. Iteration sind die Koeffizienten im Rahmen der Rechengenauigkeit exakt:

a(x) = x² + (10/3)x + 1.

Die Lösungen dieser quadratischen Gleichung sind

x₁ = −1/3 und x₂ = −3.

Beide sind daher auch Nullstellen von f(x). Nach dem Abdividieren dieser Faktoren kann dasselbe quadratische Polynom als Startwert für die Suche nach den nächsten Nullstellen verwendet werden.

Konvergenz und Grenzen

Die Abhängigkeit vom Startwert lässt sich ähnlich wie bei einem Newton-Fraktal darstellen. Jedem Punkt (u,v) wird dabei der Startfaktor

a(x) = (x − u)² + |v|v = x² − 2ux + (u² + |v|v)

zugeordnet. Für v > 0 beschreibt ein gefundener Faktor das komplex konjugierte Nullstellenpaar u ± iv; für v < 0 beschreibt er das reelle Paar u ± v. Besitzt f(x) insgesamt k komplexe Nullstellenpaare und s = n − 2k reelle Nullstellen, gibt es k entsprechende quadratische Faktoren oberhalb der x-Achse und s(s − 1)/2 mögliche Faktoren unterhalb der x-Achse, weil je zwei reelle Nullstellen kombiniert werden können.

In den dargestellten Konvergenzbildern für (u,v) ∈ [−3,3]² kennzeichnen weiße Punkte Startwerte, die schon nach einem Iterationsschritt einen guten Faktor liefern. Farbige Bereiche zeigen, zu welchem Faktor die Iteration konvergiert; dunklere Farben bedeuten mehr notwendige Schritte. Bei schwarzen Punkten wurde innerhalb der ersten 100 Schritte kein hinreichend genauer Faktor gefunden oder die Iteration divergierte.

Die Konvergenz ist im Allgemeinen nicht garantiert und kann langsam sein. Dies zeigt sich besonders beim Polynom f(x) = x⁵ − 1. Teilweise lässt sich das Verhalten verbessern, indem man durch Hinzufügen der Nullstelle x = 0 zum geraden Grad übergeht und stattdessen f(x) = x(x⁵ − 1) = x⁶ − x betrachtet. Praktisch sinnvoller ist es bei Polynomen ungeraden Grades häufig, zunächst mit Newton-Verfahren, Sekantenverfahren oder Regula falsi eine reelle Nullstelle zu bestimmen und das Bairstow-Verfahren anschließend auf das reduzierte Polynom geraden Grades anzuwenden.