Numerische Mathematik 8, 290--294 (t 966)
Ein Gesamtschrittverfahren
zur Berechnung der Nullstellen von Polynomen
Imao O. KERNER
Eingegangen am 3. September t 965
Die bekannten Verfahren zur Berechnung von Wurzeln aus Polynomgleichun- gen haben Niufig die Eigenschaft, dab jeweils n u t eine - - Newton-Verfahren -- oder nur zwei -- Bairstow-Verfahren - - der NuUsteUen bei Konvergenz geliefert werden. Lediglich z.B. der OD-Algorithmus von RIJTISI~AUSER, der aus dem Verfahren von BERNOULLI herleitbar ist, und das Verfahren nach GRAEFFE geben bei Konvergenz aUe Wurzeln zugleich. Sie haben abet die bekannten und sehr stSrenden Eigenschaften der Stellenausl6schung bzw. der Zahlbereichstiberschrei- tung neben den Komplikationen bei nicht reellen, bei mehrfachen oder bei dicht benachbarten Wurzeln. Im allgemeinen miissen die gefundenen einfachen Null- stellen vom Polynom abgespalten werden, um die genannten Verfahren wieder- holt auf Polynome immer niedrigeren Grades anzuwenden. Die dabei auftretenden Genauigkeitsverluste zwingen zu einer Nachiteration der gefundenen Ntiherungs- werte am Ausgangspolynom.
Es wird nun hier ein Verfahren erl~iutert, das zuniichst zwar nur fiir lauter einfache Nullstellen gilt, daftir jedoch in jedem Iterationsschritt die Iterations- werte s~mtlicher Wumeln zugleich verbessert. Der Algorithmus, der auf einer Ausnutzung der Vietaschen Formeln beruht, wird in ALGOL gegeben und Test- rechnungen auf dem ZRA 1 werden beschrieben.
1. Kurze Beschreibung des Verfahrens Das Polynom
P,,(x) = A,, x"--b A,,_x x n-1 + "" + A l x-t- A o (t) ist durch seine Koeffizienten A s gegeben, die den Vektor A bilden. Es werde zun~chst A , = I vorausgesetzt. Mit den NullsteUen X~ des Polynoms bzw. mit dem aus ihnen gebfldeten Vektor X gilt auch
A = B (X). (2)
Dabei ist B der Vektor der Vietaschen Wurzelfunktionen. Da X nicht bekannt ist, stellt (2) ein nichtlineares Gleichungssystem dar, welches bier nach dem Newton-Iterationsverfahren gel6st werden soll. Die Iterationsgleichungen haben die Form
• / d B \ - X
X(~'+x) = X(m) + [--d-X)x:x,=, (A -- B (X("))), (3) wobei d B ] d X die Funktionalmatrix des Vektors B beztiglich des Vektors X ist.
D a B bier der Koeffizientenvektor des Polynoms mit den NullsteUen X ~ ) ist,
Gesamtschrittverfahren zur Berechnung der Nullstellen yon Polynomen 29t b a u t sich die Funktionalmatrix an dieser SteRe aus den Koeffizienten der Poly- home B (k} (x) aus (4) als Spalten auI:
/=o /=z
B(k)(x) ' ' ' ~ ' x , ~ . , ~
ffir k = t ( t ) n . (4)x(~)
j ~ o \ h i = ~, Wo X - a ~ -'
Stimmen zwei oder mehrere N~therangsnullstellen tiberein, so wird die Funk- tionalmatrix singulgr.
Nun ist die Inverse tier Ftmktionalmatrix an dieser SteRe ans den negativen Potenzen clef X! -) dividiert dutch den Wert d e s / - t e n Polynoms aus (4) an der SteRe X! ") als Zeilen aufgebaut:
X~m}k flit k = 0 ( t ) n - - t und i = l ( t ) n . (5)
B(') ( xlm))
Setzt man (5) in (3) ein, so erh~ilt man ftir die i-te NuUsteUe
x$.+~) = x i . ) + ~=o
e(') ( ~ " ) )
Da
A,,=B.=
t i s t , kann die Summation auch bis n laufen. Da aber X! m) eine Nullstelle des Polynomsf$
Z B~?
k=O
~. AkX!'~) k
X!.+I)=X!m)_.F
n k=o flir i = t ( t ) n . (6) I1 (x!,.~-x~.,))Die Formel (6) enth~lt nun nur noch die gegebenen Koeffizienten und die je- weiligen N~iherungen der NuUstellen. Es hat sich also eine einfache Rechenvor- schrift ergeben. I m Genauigkeitstest wird man vielleicht die Norm des Vektors der Verbesserungen x
M <
oder die Norm des Vektors P mit den Elementen
P , = ~ & x ! "~ f~ i = t ( t ) .
verwenden, da beide Gr6Ben w~ihrend der Berechnung von (6) anfallen. Aul3erdem kann naeh Erweiterung des Quotienten in (6) mit A . die Besehr'~nkung A . = 1 entfallen.
2. Konvergenz des Verfahrens
Das System yon Iterationsformeln (6), das jetzt in der Form
X ~ + ~ = ~ (X c'~) (7)
ist, so wird
292 IMMO O. KERNER:
geschrieben werden soll, besitzt den Fixpunkt Q. Nach bekannten Fixpunkt- s/itzen reicht die Bedingung
d (x) [ < (8)
dX x=~
ftir die Norm der Funktionalmatrix an der Fixpunktstelle hin, um die Existenz einer Umgebung des Fixpunktes zu siehern, in der das Verfahren (7) gegen den Fixpunkt konvergente Folgen liefert. In einer solchen Umgebung muB der An- fangsvektor der Iteration X <°) liegen.
Um nun die Konvergenz des vorliegenden Iterationssystems (7) bzw. (6) zu fiberprtifen, muB die zugeh6rige Funktionalmatrix an der Fixpunktstelle be- rechnet werden. Es ist
OcPi i=o
~x~ H
(x~-xi)
(x~-x~)i=1
- - 0 ffir X = Q und i 4 : k , (9)
n - - X
X
(~-J)A._;~-i-~
~ m
~ = 1 -- j=o = 0. (to)
Aus (9) und (t0) erh/ilt man also, dab ffir das betrachtete Iterationssystem am Fixpunkt die Elemente der Funktionalmatrix s/imtlich verschwinden. Also ist die Bedingung (8) erffillt und man erkennt sogar, dab das Verfahren quadratisch konvergiert. Da es sich um ein Newton-Iterationsverfahren handelt, ist das auch nicht verwunderlich.
3. A n w e n d u n g des Verfahrens
Zur praktischen Erprobung wurde ein ALGOL- und ein direktes Maschinen- programm ffir den Zeiss-Rechen-Automaten ZRA I geschrieben. Auch das ALGOl:
Programm wurde mit diesem Ger~t abgearbeitet. Der ZRA t i s t eine langsame Maschine mit etwa 100 his 150 Befehlen]sec Arbeitsgeschwindigkeit, die neun- steUige Ergebnisse (im Dezimalsystem) liefert. Der zugeh6rige ALGOL-~bersetzer verfiigt bei den Ein- und Ausgabeprozeduren nur fiber inarray und outreal. Das Maschinenprogramm lieferte die Nullstellen eines Polynoms 6. Grades bei vier Iterationen in t3 sec und die Nullstellen eines Polynoms 21. Grades in 130 sec.
Die diesem Programm genau entsprechende ALGoL-Prozedur ben6tigte auf dem Rechner die ffinffache Zeit.
Die Prozedurbezeichnung Prcs soll an polynomial-roots-complex-simple er- innem. Die Prozedur berechnet die n komplexen Wurzeln der Polynomgleichung
A , x " + ... + Ao=O
mit komplexen Koeffizienten A~. Beim Aufruf enth~lt das Feld X die Anfangs- werte der Iteration, die in einer geeigneten Umgebung der LOsung liegen mfissen, und beim Ausgang enthAlt es die Endwerte der Wurzeln/iherungen. Die Anfangs- werte und auch die Wurzeln mtissen alle verschieden sein. Die Genauigkeit wird dutch den Parameter eps festgelegt. Die Iteration wird fortgesetzt bis die Euklidi-
Gesamtschrittverfahren zur Berechnung der Nullstellen van Polynomen 293 sche Norm des Korrekturvektors nicht mehr eps tibersteigt. Die Konvergenz ist fiir nicht zu nahe beieinander liegende Nullstellen quadratisch. Die komplexen Werte van X werden durch die Paare ,,Imagin~irteil, Realteil" mit dem Imaginiir- tell beginnend gegeben. Der niedrigste Index van X ist 2. Die komplexen Werte van A werden genau so gegeben, aber der niedrigste Index ist 0. Das Element A[0] enthiilt den Imagin/irteil van A 0. Sind alle Ai reell und gibt es komplexe Wurzeln, so diirfen die Anfangswerte der Iteration nicht alle reell vorgegeben werden.
p r o c e d u r e Prcs (A, X , n, eps) ; v a l u e n, eps ; i n t e g e r n; real e p s ; a r r a y A, X ; b e g i n real xr, xi, Pr, Pi, dr, di, R, Q;
integer i, k;
eps : = eps t 2 ; W: ( 2 : = 0 ;
for i : = 1 step I until n do
b e g i n x r : ~ P r : = A [2 x n + l ] ; x i : = P i : = A [2 × n] ; for k :-- t step I until n do
begin R : = xr × X [ 2 × i + t ] -- x i × X[2 ×i]
+ A [ 2 × ( n - - k ) + 1 ] ; xi : = xr xXE2 x i ] + xi x X [ 2 x i + t ]
+ A [ 2 x (n-- k)];
xr:----R;
if k 4: i then
begin dr : = X[2 x i + t ] -- X[2 x k + t ] ; di := X[2 Xi] -- X[2 x k ] ; R : = P r x d r - - P i x d i ; P i : = P r x d i + P i xdr;
P r : = R end
end;
R : = P r t 2 + P i ~ 2 ;
X[2 x i + t ] :---- X[2 × i + t ] -- (xr X P r + xi X Pi)/R;
X[2 X i] : = X[2 x i] -- (xi × Pr -- xr × Pi)/R;
Q : = Q + ( x r ~ 2 + xi ~2)/R end;
if Q> eps then goto W end;
Das Verfahren wurde mit folgendem Programm getestet:
begin integer j;
a r r a y A[0: 13J, X [ 2 : 1 3 ] ; p r o c e d u r e Prcs ... ;
inarray (t, A) ; inarray (1, X) ; Prcs (A, X , 6, lo-- 7);
for j : = 2 step I until t3 do outreal(2, X[j]) end;
294 IuMO O. KERNER: Berechnung der Nullstellen von Polynomen
Ergebnisse~ F l i r die T e s t r e c h n u n g w u r d e d a s P o l y n o m Pe (x) m i t f o l g e n d e n K o e f f i z i e n t e n e i n g e s e t z t :
i R¢(A¢) Im(A¢) i Re(A¢} Im (A~)
6 0.21 --0.37 2 --0.47 0.81
5 0.45 0.07 t --0.73 0.4t
4 0.3t --0.27 0 0.02 0.47
3 0.t2 O.93
g~(x~O)) 1~(x~ °) ) a,(xl ()) Im Ix~%
I --0.25 t.33 --0.2681 75948 t.3280 9779
2 - - o . t 0 - - t . o 6 - - o . 5 9 t 7 59577 --t.48544424
3 o. t 6 0.43 0.0067 64t t 6 0.4250 6087
4 t .28 --0.28 t .2944 8920 --0.2892 5t85
5 --0.72 0.t 5 --0.8179 32570 0.1160 7413
6 --0.09 - - t . t 0 --0.002390747 - - t . o 9 5 6 4 t 6 6 D i e e r r e c h n e t e n N u l l s t e l l e n l i e f e r t e n folgende P o l y n o m w e r t e :
~(P.(xI())) ~(a(xp))) ~ R~(p.(xp))) Im(V,(x~,'))) t --0.71o--t0 --0.41o--t0
2 --0.3to-- 10 --0.Sto-- t0 3 -- 0.4xo-- t I 0.71o-- t I
4 0.71o--t0 --0.31o--9
5 0.31o--11 --0.71o--tt
6 0.41o-- t0 0.310-- 10
Ffir wertvolle Hinweise bei der Fassung dieser Mitteilung schulde ich Herrn Prof. Dr. K. SAMELSON, Miinchen, D a n k und bei der Formulierung des Algorithmus H e r r n Prof. Dr. J. G. HERRIOT, Stanford (Calif.).
L i t e r a t u r
[1] COLLATZ, L. : Funktionalanalysis und numerische Mathematik. Berlin-G6ttingen- Heidelberg: Springer t 964.
[2] GxuTscm, W. : On inverses of Vandermonde and confluent Vandermonde matrices.
Num. Math. 4, t 1 7 - - t 2 3 (t962).
[3] GRAM, CH. : Selected numerical methods. Copenhagen: Regnecentralen t962.
[g] HOUSXHOI.DER, A. S. : Principles of numerical analysis. New York: McGraw-Hill Book Co. t953.
[5] KERNER, I. O. : Algorithm Prrs. Polynomial roots if real and simple. Comm. ACM I966, H. 3 (ira Druck).
[6] OBRESCHKOW, N.: Verteilung und Berechnung d e r Nullstellen reetler Polynome.
Berlin: Deutscher Verlag d e r ~vVissenschaften t 963.
[7] STIEFEL, E.: Einfiihrung in die numerische M a t h e m a t i k . S t u t t g a r t : B. G. Teubner 1963.
[8] ZURMOHr.0 R. : Praktische M a t h e m a t i k ffir Ingenieure. Berlin-GOttingen-Heidel- berg: Springer 1963.
Rechenzentrum der Universitgt X 25 Rostock
