Mitteilungen aus dem Institut angewandte MathematlCS for AN DER EIDGENÖSSISCHEN TECHNISCHEN HOCHSCHULE IN ZüRICH HERAUSGEGEBEN VON PROF. DR. E. STIEFEL Nr·7 Der Quotienten-Differenzen-Algorithmus von Heinz Rutishauser Professor an der Eidgenössischen Technischen Hochschule in Zürich BIRKHAuSER VERLAG BASEL/STUTTGART 19,,7 ISBN 978-3-7643-0323-5 ISBN 978-3-0348-7175-4 (eBook) DOI 10.1007/978-3-0348-7175-4 Nachdruck verboten. Alle Rechte vorbehalten, insbesondere das der Übersetzung in fremde Sprachen und der Reproduktion auf photostatischem Wege oder durch Mikrofilm Birkhäuser Verlag Basel, 1957 © INHALTSVERZEICHNIS 1) Einleitung. . . . . . . . . . . . 5 1. Kapitel. Theoretische Grundlagen § 1. Problemstellung . . . . . . . 7 § 2. Der Quotienten-Differenzen-Algorithmus . 8 § 3. Die Rhombenregeln . . . . . . . . 10 § 4. Die zugeordneten Polynome p~)(z) 11 § 5. Beziehungen zur Kettenbruchtheorie 13 § 6. Schwierigkeiten bei der Bildung des QD-Schemas 16 § 7. Grundlegende Eigenschaften des QD-Algorithmus 17 § 8. Beziehungen zum BO-Algorithmus von C. LANczos 20 § 9. Beziehungen zum cg-Algorithmus . . . 22 § 10. Ein Additionstheorem für Kettenbrüche . . . . . 23 II. Kapitel. Anwendungen des QD-Algorithmus § 1. Umwandlung einer Potenzreihe in einen Kettenbruch 26 § 2. Summation schlecht konvergenter Reihen 27 § 3. Auflösung von algebraischen Gleichungen 29 § 4. Die progressive Form des QD-Algorithmus 30 § 5. Auflösung algebraischer Gleichung mit Hilfe des progressiven QD-Algorithmus . . . . . . . . 31 § 6. Die Wronskische Formel. . . . . 33 § 7. Bestimmung komplexer Nullstellen 34 § 8. Quadratische Konvergenz des QD-Algorithmus 36 § 9. Massnahmen bei Division durch Null. . . . . 41 § 10. Interpolation durch Exponentialsummen. . . 43 III. Kapitel. Bestimmung der Eigenwerte und Eigenvektoren einer Matrix mit H ilje des Quotienten-Dijjerenzen-Algorithmus § 1. Die Bestimmung der Eigenwerte . . . . . . 49 § 2. Das Problem der Eigenvektorberechnung . . . 52 x1 § 3. Rekursive Berechnung der Vektoren 2/-1) , y~2/-1) 53 § 4. Ein quadratisch konvergentes Verfahren zur Eigenvektorbestim- mung . . . . . . . . . . . . . . . . . . . . . . . . 59 § 5. Eigenwerte und Eigenvektoren unendlicher symmetrischer Ma- trizen . . . . . . . . . . . . . . . . . . . . . . . . 62 IV. Kapitel. Anhang § 1. Die LR-Transformation . . . . . . . . . . . . . 65 § 2. Ein kontinuierliches Analogon zum QD-Algorithmus. 69 § 3. QD-Relaxation 73 Literaturverzeichnis 74 1) Bei Hinweisen im Text, die sich auf eine Formel oder einen Paragraphen eines andern Kapitels beziehen, wird die Nummer des betreffenden Kapitels beigefügt; zum Beispiel (H, 7) = Formel 7 aus Kapitel II. Einleitung Im Anschluss an eine praktische Anwendung des BO-Algorithmus (Biortho gonalisierungs-Algorithmus von C. LANCZOS [4], [5]1) machte mich Herr Prof. E. STIEFEL, ETH, auf das Problem aufmerksam, die höheren Eigenwerte direkt aus den sogenannten Schwarzsehen Konstanten zu bestimmen, das heisst ohne den Umweg über die Orthogonalisierung. Auf diese Anregung hin entwickelte der Verfasser einen Algorithmus, der die gestellte Aufgabe löst. Allerdings gab bereits A. C. AITKEN [1] eine Methode an, welche haupt sächlich zur Auflösung algebraischer Gleichungen gedacht war, aber auch die Bestimmung höherer Eigenwerte aus Schwarzsehen Konstanten gestattet. Ferner stammt von C. LANCZOS ein Algorithmus2) zur Bestimmung des charak teristischen Polynoms einer Matrix aus Schwarzsehen Konstanten. Überdies entwickelte J. HADAMARD in seiner Dissertation [2] eine Methode zur Bestim mung der Pole einer durch ihre Potenzreihe gegebenen Funktion. Er hat damit, wie § 1 zeigen wird, auch das eingangs erwähnte Eigenwertproblem gelöst. Wenn hier das schon gelöste Problem nochmals aufgegriffen wird, so geschieht dies deshalb, weil der entwickelte Algorithmus eine Reihe von weiteren An wendungen gestattet und insbesondere auch wertvolle Beziehungen zur Ketten bruchtheorie vermittelt3). Die Arbeit gliedert sich in drei Kapitel, von denen sich die Kapitel I und n mit Theorie und Anwendungen befassen, während III eine Ausdehnung des QD-Algorithmus auf Vektoren behandelt. Schliesslich folgt ein Anhang über verwandte Methoden (insbesondere die LR-Transformation). Die Kapitel I, n, In sind einzeln bereits in der ZAMP erschienen'), doch ist zu beachten, dass I und n zum Teil erhebliche Veränderungen erfahren haben. 1) Die Ziffern in eckigen Klammern verweisen auf das Literaturverzeichnis, Seite 74. 2) Es handelt sich nicht um den BO-Algorithmus, vgl. vielmehr Kapitel VI bei [4] oder S. 173-179 bei [5]. 3) Herrn Prof. STIEFEL verdanke ich auch die Anregung zur Vereinfachung einiger Beweise mit Hilfe der Kettenbruchtheorie. ') Kapitel I in ZAMP 5, 233-251 (1954); II in ZAMP 5, 496-508 (1954); III in ZAMP 6, 387 bis 401 (1955). I. KAPITEL Theoretische Grundlagen § 1. Problemstellung Sei A eine n-reihige Matrix, AT ihre Transponierte, ferner seien xo, Yo zwei Vektoren, von denen Xo in bezug auf die Matrix A und Yo in bezug auf die Matrix AT in allgemeiner Lage sei, das heisst, im Koordinatensystem der Eigen und Hauptvektoren von A (bzw. AT) verschwinde keine Komponente von Xo (bzw. Yo). Ich bilde nun aus Xo und Yo die unendlich vielen Vektoren Xv = A' Xo und Y. = (A 1'). Yo sowie die innern Produkte SMv = (xl" Y.), die bekanntlich nur von der Summe der Indizes p" '11 abhängen und als Schwarzsehe Konstanten der Matrix A in bezug auf die Anfangsvektoren Xo und Yo bezeichnet werden. Die mit diesen Schwarzsehen Konstanten s,. gebildete Funktion 00 fez) = }; Z:~1 (1) o ist nun rational, und ihre Pole sind die Eigenwerte von A, denn es ist fez) = (B xo, Yo) mit B = (z E - A)-I. Damit ist das gestellte Problem, nämlich die Bestimmung der Eigenwerte einer Matrix A aus Schwarzsehen Konstanten, zurückgeführt auf die Bestim mung der Pole einer rationalen Funktion f( S ) = ZSo + ZS2I- + .... Ist jedoch A eine unendliche Matrix oder ein Integraloperator, so ist die Funktion fez) zwar nicht mehr rational, aber ihre Pole sind immer noch Eigen werte von A. In gewissen Fällen kann man beweisen, dass fez) für z =1= 0 mero morph ist, beispielsweise dann, wenn die Schwarzsehen Konstanten einem selbstadjungierten und volldefiniten Eigenwertproblem entstammen, für wel ches der Entwicklungssatz gilt. Die dargelegten Gründe rechtfertigen es. wohl ohne weiteres, nicht das Pro blem der Eigenwertbestimmung, sondern wie J. HADAMARD [2] das etwas all gemeinere Problem der Bestimmung der Pole einer durch eine Reihe (1) gege benen Funktion fez) zu behandeln, wobei die s" beliebig komplex sein dürfen. Darüber hinaus wird in gewissen Fällen auch noch die Partialbruchzerlegung der Funktion fez) gesucht sein, falls eine solche existiert. 8 I. Kapitel. Theoretische Grundlagen § 2. Der Quotienten-Differenzen-Algorithmus (QD-Algorithmus) Für den Fall, dass die erzeugende Funktion .E00 fez) = z:~T o I I Al ausserhalb eines gewissen Kreises z = R genau einen einfachen Pol besitzt, gab bereits D. BERNOULLI die Formel /,L I = l1' m -S-.+1- (2a) 1'-+00 Sv an. Zur Bestimmung weiterer Pole dienen die Formeln (21) bei AITKEN [lJ; beispielsweise erhält man unter der Voraussetzung IA ll> IA 21> IA 31f ür den Pol A die Beziehung 2 1S .+1 s.+21' , S.+2 S.+3 (2b) und entsprechende Formeln können für die weiteren Pole aufgestellt werden. Nun gilt aber die Identität I s. S.+1 I , S.+1 Sv+2 = S (sV+2 _ ~V+1 ) Sv v + 1 SV+1 Sv' der wir entnehmen, dass man A auch durch folgenden Prozess aus den Schwarz 2 sehen Konstanten s. = S\,,)6) ermitteln kann: Man bildet die Quotienten q\,,) = s'!+ 1) / s'i), dann die Differenzen dieser q'i) : dr) = q'i + 1) - q'i), und schliess lieh die Produkte Dann ist A2 nach (2b) der Grenzwert der Quotienten q~) = S~+l)/S~) ('1'+00), Al genau so wie nach (2a) der Grenzwert von qr) ist. Für die praktische Rechnung wird man die Zahlen sr), qr), ... unter einanderschreiben (vgl. das Beispiel auf Seite 9); man spricht dann von einer ({ Quotientenkolonne » q"!) bzw. von einer ({ Differenzenkolonne » usw. Man wird nun vermuten, dass sich durch konsequente Fortsetzung dieses Prozesses auch die weiteren Pole A3, A4, ••• als Grenzwerte entsprechender Quo IAll>IA21>IA31>IA41>"')' tienten q~), q~), ... ergeben (vorausgesetzt, dass Dies ist in der Tat der Fall, nur muss man die weiteren Differenzenkolonnen d);') (0' = 2, 3, ... ) vor der Multiplikation mit S);'+l) modifizieren; es ergibt sich 6) Wir bezeichnen die Schwarzsehen Konstanten s. fortan auch gelegentlich mit s~·J. § 2. Der Quotienten-Differenzen-Algorithmus 9 dann die folgende allgemeine Vorschrift zur Aufstellung des QD-Schemas zu den gegebenen Reihenkoeffizienten (bzw- Schwarzsehen Konstanten) Sv = st): Ausgehend von den st) berechne man nacheinander ( (J = 1,2,3, ... ) q~) = s(v+1) /s(v) Quotientenreihe der s~) , (] (] d(v) = q~'+I) _ q~) Differenzenreihe der q~) , (] (3) e(v) = d(v) I e(v+ 1) modifizierte Differenzen 7) , (] a T G-l s(v) = S(v+l) . e(v) neue s-Reihe . a+l a a Die so erhaltenen Grössen werden wie folgt in ein Schema eingeordnet (QD Schema der sv) : s(O) 1 qiO) S(l) d~O) = eiO) S(O) 1 2 q~l) q~O) S(2) d?) = e~l) s(1) d(O) e(O) s(O) 1 2 2 2 3 qi2) q~]) q(O) 3 S(3) di2i = e~2) S(2) d(l) e(1) s(l) 1 2 2 2 3 Beispiel zum QD-Algorithmus. Die s(v) selen die sogenannten Fibonaccischen 1 Zahlen; es ist dann j(z) = -z 2 - zz - 1: s(v) qiv) div) = e~) s(V) q~v) d(v) e(v) s(V) 1 2 2 2 3 1 1 1 1 1 2 -1 2 -0,5 -1 +0,5 ° ° 1,5 -0,5 3 0,166667 0,5 -0,166667 ° ° 1,666667 - 0,666667 5 - 0,066667 - 0,333333 + 0,066667 ° ° 1,6 -0,6 8 +0,025 +0,2 - 0,025 ° ° 1,625 - 0,625 13 -0,009615 -0,125 1,615384 I 21 t t Ä1 s=::; 1,618034 Ä2 s=::; - 0,618034 7) Für (1 = 1 stimmen die e<il mit den d~) überein, also ist immer e~) = 0 zu setzen. 10 I. Kapitel. Theoretische Grundlagen An diesem einfachen Beispiel tritt bereits eine Erscheinung zutage, die sich als allgemeingültig erweisen wird: Wenn f(z) eine rationale Funktion ist, deren Nenner den Grad n besitzt, so bricht das QD-Schema nach der n-ten modifizierten Ditterenzenkolonne ab, ganz ähnlich, wie das gewöhnliche Ditterenzenschema eines Polynoms nach der n-ten Ditterenz abbricht. Diese Eigenschaft erlaubt uns, das QD-Schema in einem solchen Fall auch von oben nach unten fortschreitend aufzubauen, wenn die in der obersten Schrägzeile stehenden Werte so' qiO), eiO), q~O), ••• , e~o~ l' q~O) bekannt sind. Dieses Vorgehen, welches einem wohlbekannten Prozess beim gewöhnlichen Differenzenschema völlig analog ist, wird vor allem für die An wendungen von grosser Bedeutung sein. § 3. Die Rhombenregeln Falls man die s~)-Werte (er> 1) nicht benötigt, kann man auch mit der folgenden komprimierten Vorschrift arbeiten (ausgehend von e~) = 0 und q~) = Sv + 1/SV) : + e(V) = e(v+1) q(v+1j _ q(v) I o 0-1 a 0 ' e(v+l) (er = 1,2, ... ) (4) q(v) = q(V+l) • ~ __ a+1 a (,,) ea und erhält so das QD-Schema in komprimierter Form siO) s(1) e(O) 1 1 q~1) q~O) S(2) e(1) e(O) 1 qi 1 2 2) q~1) q~O) S(3) e(2) e(l) 1 1 2 Die obigen Formeln (4), nach denen das komprimierte QD-Schema auf gebaut ist, können nach einem Vorschlag von Herrn Prof. E. STIEFEL als soge nannte Rhombenregeln8) formuliert werden: 8) Siehe E. STIEFEL [13], Seite 42. § 4. Die zugeordneten Polynome p~\z) 11 Jede der beiden Formeln verknüpft 4 Elemente, die im obigen Schema wie die Ecken eines Rhombus angeordnet sind: Q-Rhombus E-Rhombus q(v) a+l Dies gibt offenbar Anlass zu den folgenden Regeln: Q-Regel: Die Summe der im Q-Rhombus über der punktierten Linie stehen den Grössen ist gleich der Summe der unter dieser Linie stehenden GrÖssen. E-Regel: Das Produkt der im E-Rhombus über der punktierten Linie stehen den Grössen ist gleich dem Produkt der unter dieser Linie stehenden GrÖssen. § 4. Die zugeordneten Polynome p~) (z) Gewissermassen als Vorläufer des BO-Algorithmus hat C. LANczos das bereits erwähnte Verfahren zur Bestimmung des charakteristischen Polynoms einer Matrix aus Schwarzsehen Konstanten angegeben. Im Verlaufe dieses Prozesses, der übrigens mit dem QD-Algorithmus nahe verwandt ist9), wird eine Kette von Polynomen Po == 1, PI(A), P2(A), ... aufgebaut, deren letztes Glied Pm(A) das Minimalpolynom der Matrix A ist. Entsprechend bilden wir hier das zweidimensionale Schema der Polynome P);,) (z) vom Grade (P-Schema): (j p~o) piO) p~l) p~O) pil) P~O) p~2) p~l) p(l) pi2) 3 das wir nach der folgenden Regel aufbauen: p~) == 1 für alle v , ~ "~ (5) p~)(z) = Z P~~II)(Z) - q~) P~~l(Z) (v 0, I, ... ; I, 2,. .. ). ) Die Regel (5) kann auch schematisch dargestellt werden, so dass man sofort 9) V gl. hierzu § 8. 12 1. Kapitel. Theoretische Grundlagen ersieht, in welcher Weise ein Element im P-Schema aus den links davon ste henden Elementen p(v+!) und p(v) gebildet wird: <1-1 <1-1 p(V) q~) <1-1 ~ p~,). p(V+!) Z <1-1 Wir wollen nun für diese Polynome p~)(z) die beiden folgenden Relationen be weisen: (6) (7) a> für O. ) Beweis von (6) Wir verwenden als Induktionsvoraussetzung die Formel (6) mit 1 an (J - Stelle von also die Beziehung (J, welche für (J = 1 wegen P~) == 1 und e~) = 0 sicher für alle verfüllt ist. Wir bilden nun die Grösse die wegen (5) zu + + z [p(V+2) _ plv+1)] _ q(v+1) p(v+1) q(v) pr,,) e(v) p(v+1) a-l a-l a G-l a a-l (J G-l wird. Unter Verwendung von q~'+1) _ e~) = q~') - e~:11) [siehe (4)J ergibt sich 15 = z [P~:12) - p~:l)] - q~) [P~:11) - P~~1] + e~:11) p~:l) , oder unter Verwendung der Induktionsvoraussetzung: 15 = -z e(v+1) p(v+2) + q(v) e(v) p(V+1) + e(v+1) p(v+1) <1-1 <1-2 <1 <1-1 <1-2 <1-1 <1-1' was sich schliesslich unter Verwendung von q(v) e(v) = q(V+ 1) e(v+ 1) a 0'-1 0'-1 a-l [siehe (4)J zu + + 15 = e(V+1){_zp(V+2) q(v+1) p(v+1) p(V+1)} <1-1 <1-2 <1-1 <1-2 <1-1 reduziert. Dieser Ausdruck verschwindet aber wegen (5); was zu beweisen war.