Numerik Algorithms
Numerik Algorithms
press
Die Reihe Xpert.press vermittelt Professionals
in den Bereichen Softwareentwicklung,
Internettechnologie und IT-Management aktuell
und kompetent relevantes Fachwissen über
Technologien und Produkte zur Entwicklung
und Anwendung moderner Informationstechnologien.
Gisela Engeln-Müllges · Klaus Niederdrenk
Reinhard Wodicka
Numerik-Algorithmen
Verfahren, Beispiele, Anwendungen
123
Prof. Dr. rer. nat. Gisela Engeln-Müllges Prof. Dr. rer. nat. Klaus Niederdrenk
Fachbereich Maschinenbau und Mechatronik Fachbereich Chemieingenieurwesen
Fachhochschule Aachen Fachhochschule Münster
Goethestraße 1 Stegerwaldstraße 39
52064 Aachen 48565 Steinfurt
[email protected] [email protected]
Die Neuauflage fasst die 5. Auflage von „Numerische Mathematik für Ingenieure“ von Gisela Engeln-
Müllges und Fritz Reutter(ursprünglich Bibliographisches Institut/Brockhaus AG, Mannheim) und die
8. Auflage von „Numerik-Algorithmen“ von Gisela Engeln-Müllges und Fritz Reutter (ursprünglich
VDI-Verlag, Düsseldorf) zusammen.
ISSN 1439-5428
ISBN 3-540-62669-7 Springer Berlin Heidelberg New York
Dieses Werk ist urheb errecht lich geschützt. Die dadurch b eg rü ndeten Rechte, in sbe-
sondere die der Übersetzung, des Nachdrucks, des Vortrags, der Entnahme von Abbildungen und
Tabellen, der Funksendung, der Mikroverfilmung oder der Ver vielfältigung auf anderen Wegen und der
Speicherung in Datenverarbeitungsanlagen bleiben, auch bei nur auszugsweiser Verwertung,
vorbehalten. Eine Vervielfältigung dieses Werkes oder von Teilen dieses Werkes ist auch im Einzelfall nur
in den Grenzen der gesetzlichen Bestimmungen des Urheberrechtsgesetzes der Bundesrepublik
Deutschland vom 9. September 1965 in der jeweils geltenden Fassung zulässig. Sie ist grundsätzlich
vergütungspflichtig. Zuwiderhandlungen unterliegen den Strafbestimmungen des Urheberrechtsgesetzes.
Springer ist nicht Urheber der Daten und Programme. Weder Springer noch die Autoren übernehmen
die Haftung für die CD-ROMs und das Buch, einschließlich ihrer Qualität, Handels- und
Anwendungseignung. In keinem Fall übernehmen Springer oder die Autoren Haftung für direkte,
indirekte, zufällige oder Folgeschäden, die sich aus der Nutzung der CD-ROMs oder des Buches ergeben.
Das vorliegende Werk steht in einer langen Tradition. Die erste Auflage erschien 1974 mit
den Autoren Gisela Engeln-Müllges und Fritz Reutter unter dem Titel Formelsamm-
”
lung zur Numerischen Mathematik mit FORTRAN-Programmen“. Der Inhalt des Buches
orientierte sich an Vorlesungen über Numerische Mathematik für Studierende der Inge-
nieurwissenschaften.
Bei den nachfolgenden Auflagen, an denen Fritz Reutter bis zu seiner schweren Erkran-
kung 1983 mitarbeitete, nahm der Umfang des Buches stark zu, und die Programm-
anhänge wurden um die Sprachen Turbo Pascal, C, Modula 2 und Quick-Basic erweitert.
Dies hing zusammen mit der raschen Entwicklung der Computertechnologie und ihren
Impulsen auf die Numerische Mathematik, die zu einer enormen qualitativen wie auch
quantitativen Zunahme numerischer Verfahren führte.
Grundlage des vorliegenden Buches ist die achte Auflage des Buches von Gisela Engeln-
Müllges und Fritz Reutter mit dem Titel Numerik-Algorithmen“ (VDI-Verlag, 1996),
”
die erstmals eine CD-ROM mit Fortran 77/90-, Turbo Pascal- und ANSI C-Programmen
enthielt.
Wie bisher wird auch in dieser neunten Auflage von den unten genannten Verfassern
besonderer Wert gelegt auf die Erläuterung der Prinzipien der behandelten Verfahren
der Numerischen Mathematik, auf die exakte Beschreibung leistungsfähiger Algorith-
men und die Entwicklung und Dokumentation zugehöriger Programme. Die ausführli-
che Darstellung der mathematischen Grundlagen, ergänzt durch viele Abbildungen und
durchgerechnete Beispiele, soll den Zugang zur Numerischen Mathematik erleichtern und
das Verständnis des algorithmischen Vorgehens fördern. Zahlreiche Beispiele aus inge-
nieurwissenschaftlichen Anwendungen belegen die Notwendigkeit der behandelten nume-
rischen Verfahren.
Das Buch wendet sich in erster Linie an Ingenieure sowie an Naturwissenschaftler und
Informatiker in Studium und Beruf, ferner an Lehrkräfte für mathematisch-naturwissen-
schaftlichen Unterricht. Deshalb ist besonders auf die Verwendbarkeit der behandelten
Themen in der Praxis geachtet worden.
Beigefügt ist dem Buch eine CD-ROM, die umfassend getestete C-Programme zu nahezu
allen angegebenen Algorithmen enthält, sowie weitere Software, zu der Informationen auf
den folgenden Seiten zu finden sind.
VIII Vorwort
Die inhaltliche Erweiterung des Buches erfordert es, die bisher in ihm enthaltenen nu-
merischen Verfahren für Anfangs- und Randwertprobleme bei gewöhnlichen Differential-
gleichungen in einen weiteren Band zu verlagern. Dieser wird, einem häufig geäußerten
Wunsch folgend, zusätzlich Numerik partieller Differentialgleichungen sowie wichtige sto-
chastische Methoden, insbesondere statistische Schätz- und Prüfverfahren enthalten und
ebenfalls bei Springer erscheinen. Verfasser dieses Buches sind Gisela Engeln-Müllges,
Klaus Niederdrenk und Wieland Richter.
Ganz besonders herzlich danken wir den Autoren der Programme und ebenso Doris und
Uli Eggermann, die mit Hilfe des Satzprogramms LATEX das reproduktionsreife Manu-
skript sowie eine Vielzahl der Abbildungen mit äußerstem Engagement und sehr großem
Geschick fertig gestellt haben. Darüber hinaus hat uns Uli Eggermann tatkräftig beim
Korrekturlesen und beim Rechnen der Beispiele zur Kubatur unterstützt.
Ein herzlicher Dank gilt nicht zuletzt Herrn Hermann Engesser vom Springer-Verlag für
die hervorragende und effektive Zusammenarbeit.
Informationen
zur CD-ROM mit C-Programmen u.a. (CD-1)
Auf dieser CD befinden sich Ansi-C-Quellen von Unterprogrammen zu den meisten der im
Buch angegebenen Algorithmen. Diese Quellen können compilerunabhängig in eigenen
C-Programmen verwendet werden. Mittels beigefügter Makefile-Datei ist man in der
Lage,
Systemvoraussetzungen
Die C-Programme der CD-ROM sind bei der FH-Aachen auch in C++ erhältlich. Infor-
mationen über Lizenzen zu den C++ -Programmen sowie über Campuslizenzen zu den
C- und C++ -Programmen können per e-Mail angefordert werden:
[email protected]
X Informationen zur beigefügten Software
Systemvoraussetzungen
Betriebssystem: Windows (98, 2000 oder XP)
Arbeitsspeicher: mindestens 256 MB RAM
Browser: Netscape ab 6.2, Internet Explorer ab 5.5, Mozilla ab 1.4
Java: aktivieren (ab Version 1.2); wenn von CD gestartet
(aus ReadMe.htm), wird mitgeliefertes Java verwendet.
Informationen
zum Expertensystem CurveTrac“ (CD-1)
”
CurveTrac ist ein Baukastensystem von numerischen Anwendungen, die auf den C-Pro-
grammen aus diesem Buch basieren. Die Benutzer-Oberfläche ist in Java programmiert.
Dadurch ist gewährleistet, dass die Anwendung auf den gängigsten Betriebsystemen ge-
nutzt werden kann. Getestet wurde CurveTrac unter Windows und Linux. In der hier
vorliegenden Ausbaustufe ist das Modul Höhenlinien“ eingebunden. Informationen über
”
weitere Module können Sie bei [email protected] anfordern. Dies sind im
Einzelnen:
− Expertensystem zur Numerischen Mathematik
− Berechnung der Schnittkurve zweier Flächen im Raum
− Splinemodul (Flächensplines und Kurvensplines)
− u. a. m. . . .
− Flächensplines
− Nullstellenverfahren
− Differentialgleichungen
− Matrizen
− Statistik
− u. a. m. . . .
− Quadratur
− Nullstellen
− Splines
− Differentialgleichungen
Informationen
zur Demo-CD-ROM NUMAS“ (CD-2)
”
Die beiliegende Demo-CD-ROM enthält das Lernfeld Kurvensplines“ aus dem multime-
”
dialen Lehr- und Lernsystem NUMAS zur Numerischen Mathematik und Statistik. Der
XII Informationen zur beigefügten Software
Inhalt der Demo-CD-ROM entspricht den Kapiteln über Kurvensplines in diesem Buch
und informiert auch über den Gesamtinhalt des Lernsystems.
NUMAS bietet didaktisch aufbereitetes Wissen zur Numerischen Mathematik und Sta-
tistik. Es werden Inhalte angeboten, die für die Hochschulausbildung vieler Fachrichtun-
gen grundlegend sind und die sich gleichzeitig hervorragend dazu eignen, in Formen des
angeleiteten Selbststudiums aufgenommen zu werden. Lernerinnen und Lerner werden
zeitgemäß in ihren Selbstlern- und Selbstorganisationskompetenzen gefördert.
Für die Vollversion NUMAS wurden in großen Teilen die Inhalte dieses Buches und des
im Vorwort angekündigten weiteren Bandes mit Differentialgleichungen und Statistik
verwendet.
NUMAS ist aus einem vom BMBF und vom Land NRW geförderten Projekt innerhalb
der Ausschreibungen Neue Medien in der Bildung“ und Neue Medien in der Hoch-
” ”
schullehre“ hervorgegangen. An dem Projekt waren die Fachhochschule Aachen (Prof.
Dr. Engeln-Müllges), die Freie Universität Berlin (Prof. Dr. Martus), die Fachhochschule
Münster (Prof. Dr. Niederdrenk) und die Fachhochschule Südwestfalen (Prof. Dr. Rich-
ter) beteiligt.
NUMAS ist auch im Internet unter http://www.numas.de/ zu erreichen.
Der Umgang mit der CD-ROM selbst ist denkbar einfach: Nach dem Einlegen der CD-
ROM in das Laufwerk erscheint selbstständig ein Auswahlbildschirm mit den im Ge-
samtsystem zur Verfügung stehenden Lernfeldern. Von hier aus kann man sich den In-
halt ansehen und bearbeiten. Bei der Demo-CD-ROM ist nur das Lernfeld Kurvensplines
freigeschaltet.
Die CD-ROM stellt nur einen Teil der Funktionalität des Lernsystems zur Verfügung.
Zum Beispiel sind die Kommunikationswerkzeuge, die an das Internet gebunden sind,
nur im Online-System verfügbar.
Systemvoraussetzungen
Ihr Computersystem sollte folgende Bedingungen erfüllen, um einen reibungslosen Be-
trieb von NUMAS erreichen zu können:
Betriebssystem: Windows (98, 2000 oder XP)
Browser: Netscape ab 6.2, Internet Explorer ab 5.5, Mozilla ab 1.4
Javascript: aktivieren (ab Version 1.2)
Cookies: aktivieren
Um auch die interaktiven und multimedialen Elemente von NUMAS nutzen zu können,
müssen zusätzlich folgende Plug-Ins installiert sein:
Macromedia Flash Player und SUN Java (JRE) ab Version 1.4 (mit Ausnahme der Ver-
sion 1.4.1 03)
Weitere Informationen zu NUMAS können Sie unter der Mail-Adresse [email protected]
anfordern.
Bezeichnungen
√
i imaginäre Einheit i := −1
e Eulersche Zahl = 2.718 281 828 459 . . .
a für a ≥ 0
|a| Betrag von a mit |a| :=
−a für a < 0
· Norm von ·
{ak } Folge von ak
lim ak Limes von ak für k → ∞
k→∞
max {f (x)|x ∈ [a, b]} Maximum aller Funktionswerte f (x) für x ∈ [a, b]
min {|Mi | i = 1, 2, . . . , m} Minimum aller |Mi | für i = 1, 2, . . . , m
(x, y) geordnetes Paar
(x1 , x2 , . . . , xn ) geordnetes n-Tupel
f : I → IR Abbildung f von I nach IR
x → f (x), x ∈ D x wird f (x) zugeordnet für x ∈ D
f , f , f , f (4) , . . . , f (n) erste, zweite, dritte, vierte, . . . ,
n-te Ableitung von f
C[a, b] die Menge der auf [a, b] stetigen Funktionen
C n [a, b] die Menge der auf [a, b] n-mal stetig differenzierbaren
Funktionen
IRn n-dimensionaler euklidischer Raum
A = O(hq ) |A/hq | ≤ C für h → 0, C = const.
i = m(1)n m, n ∈ ZZ, m ≤ n, i = m, m + 1, . . . , n
⎧
⎨ 1 für a > 0
sign (a), sgn (a) 0 für a = 0
⎩
−1 für a < 0
x, y, z , . . . Vektoren
A, B, C , . . . Matrizen
0 Nullmatrix bzw. Nullvektor
x ×y Vektorprodukt bzw. Kreuzprodukt
E Einheitsmatrix
AT transponierte Matrix von A ⎛ ⎞
x1
x T = (x1 , x2 , . . . , xn ) transponierter Vektor zu x = ⎝ ... ⎠
xn
A−1 inverse Matrix von A
|A|, det(A) Determinante von A
[Verweis] s. im Literaturverzeichnis unter [Verweis]
o. B. d. A. ohne Beschränkung der Allgemeinheit
Ende eines Beispiels oder eines Beweises
Inhaltsverzeichnis
Sachwortverzeichnis 669
Kapitel 1
∆a = a − A
|∆a | = |a − A|
Sehr oft wird in der mathematischen Literatur ∆a bereits als absoluter Fehler und |∆a |
als Absolutbetrag des Fehlers bezeichnet. In ingenieurwissenschaftlichen Anwendungen
ist allerdings die Schreibweise in Definition 1.1 häufiger anzutreffen.
In den meisten Fällen ist die Zahl a nicht bekannt, so dass weder der wahre noch der
absolute Fehler eines Näherungswertes A angegeben werden können. Daher versucht man,
für den absoluten Fehler |∆a | von A eine möglichst kleine obere Schranke εa > 0 anzu-
geben, so dass |∆a | ≤ εa gilt.
2 1. Darstellung von Zahlen und Fehleranalyse
Streng genommen müsste man δa = ∆a /a als relativen Fehler von A bezeichnen. Das
Vorzeichen von δa gibt dann eine zusätzliche Information über die Richtung des Fehlers,
d. h. für eine positive Zahl a hat δa = (a − A)/a < 0 demnach a < A, also einen zu
großen Näherungswert A für a zur Folge.
relativer Fehler von A genannt. Da dann auch ∆a nicht exakt angebbar ist, wird man sich
wieder mit der Angabe einer möglichst guten oberen Schranke für den relativen Fehler
behelfen müssen.
Beispiel 1.6.
Für die Zahl x = π sind X = 3.14 ein Näherungswert und εx = 0.0016 eine Schranke für
den absoluten Fehler |∆x |. Also gilt nach (1.1)
3.1384 = 3.14 − 0.0016 ≤ π ≤ 3.14 + 0.0016 = 3.1416 ,
und der relative Höchstfehler ergibt sich zu
0.0016
|δx | ≤ x = ≈ 0.00051 .
3.14
Für den prozentualen Fehler folgt nach Definition 1.5
100 · |δx | = 0.051% .
1.2 Zahlensysteme
1.2.1 Darstellung ganzer Zahlen
Für jede ganze Zahl a gibt es genau eine Potenzentwicklung zur Basis 10 (Zehnerpotenzen)
der Gestalt
n
a = v · (an 10n + an−1 10n−1 + . . . + a1 101 + a0 100 ) = v · ak 10k
k=0
Es gibt keinen triftigen Grund, nur das uns gewohnte und vertraute Dezimalsystem zu
benutzen. Wegen der einfachen technischen Realisierbarkeit benutzen digitale Rechenan-
lagen durchweg das Dualsystem, das auf den zwei verschiedenen Ziffern 0 und 1 beruht.
Jede ganze Zahl a kann damit in eindeutiger Weise als Potenzentwicklung zur Basis 2
(Zweierpotenzen)
n
a =v· ak · 2 k mit ak ∈ {0, 1}
k=0
und einem Vorzeichen v ∈ {−1, 1} geschrieben werden. Und wiederum gibt nur die Stel-
lung einer Ziffer ak Auskunft über ihren Stellenwert, so dass die Angabe der dualen
Ziffern in der richtigen Reihenfolge zur Charakterisierung ausreicht:
a = v · (an an−1 . . . a1 a0 )2
Der Index 2 soll hierbei auf die Darstellung als Dualzahl verweisen; weggelassen wird nur
der Index 10 für die üblichen Dezimalzahlen. So hat man beispielsweise
2004 = (+1) · (2004)10
= (+1) · (1 · 210 + 1 · 29 + 1 · 28 + 1 · 27 + 1 · 26 + 1 · 24 + 1 · 22 )
= (+1) · (11111010100)2 .
Die Wahl β = 10 führt auf das geläufige Dezimalsystem und β = 2 auf das Dualsy-
stem. Anwendungen finden immer wieder auch die Basiswahlen β = 8 (Oktalsystem) und
β = 16 (Hexadezimalsystem mit den Ziffern 0, 1, . . . , 9, A, B, C, D, E, F)
a = v · (an an−1 . . . a1 a0 )β ,
Rechnerintern wird eine ganze Zahl a = v · |a| als Dualzahl durch eine Reihe von Bits
(“binary digits”) abgespeichert, gewöhnlich in zwei oder vier Bytes (1 Byte = 8 Bits).
Bei einer 2-Byte-Hinterlegung bedeutet das:
1. Byte 2. Byte
.
µ .. α14 α13 α12 α11 α10 α9 α8 α7 α6 α5 α4 α3 α2 α1 α0
Damit eine Subtraktion auf eine Addition über a − b = a + (−b) zurückgeführt werden
kann, werden negative Zahlen rechnerintern durch ein Komplement dargestellt, d. h. im
Fall µ = 1 bzw. v = −1 stimmen die Hinterlegungen αk nicht mehr mit den Dualzif-
fern ak der Zahl a = (−1) · (a14 a13 . . . a1 a0 )2 überein.
Vorsicht ist geboten, wenn bei einer Addition oder Subtraktion das Ergebnis nicht mehr
im darstellbaren Bereich liegt: Addiert man beispielsweise bei einer 1-Byte-Zahl [größte
darstellbare positive Zahl: 27 −1 = 127] 97 = (+1)·(1100001)2 und 43 = (+1)·(0101011)2 ,
so folgt
.
0 .. 1 1 0 0 0 0 1
.
+ 0 .. 0 1 0 1 0 1 1
.
1 .. 0 0 0 1 1 0 0
d. h. ein Übertrag beeinflusst letztendlich das Vorzeichenbit und führt auf ein völlig
inplausibles negatives Resultat! Solche Effekte treten leider häufiger auf, ohne dass vom
Rechnersystem eine Fehlermeldung oder zumindest eine Warnung ausgesprochen wird.
Deshalb ist bei Rechnerergebnissen stets eine Prüfung auf Plausibilität
erforderlich!
mit v ∈ {+1, −1} und ak , bk ∈ {0, 1, . . . , β−1}; der Index β verweist wieder auf das
zugrunde liegende Stellenwertsystem und entfällt nur im vertrauten Fall β = 10. Der
zwischen den Ziffern a0 und b1 gesetzte Punkt heißt nun β-Punkt , und man nennt
(an an−1 . . . a1 a0 )β – manchmal auch unter Berücksichtigung des Vorzeichens v – den
ganzzahligen Anteil und (. b1 b2 b3 . . .)β den gebrochenen oder fraktionierten Anteil von a.
Der Sonderfall β = 2 führt nun zur Dualbruchdarstellung einer reellen Zahl. Gibt es nur
endliche viele Nachkomma“-Stellen, so spricht man von einem endlichen β-Bruch, an-
”
dernfalls von einem unendlichen β-Bruch.
Es gilt der
Hilfssatz 1.10.
Jede rationale Zahl p/q mit teilerfremden p ∈ ZZ und q ∈ IN wird durch einen endlichen
oder durch einen unendlichen periodischen Dezimal- oder Dualbruch dargestellt, jede
irrationale Zahl durch einen unendlichen nicht periodischen Dezimal- oder Dualbruch.
Zur Umwandlung reeller Zahlen in ein β-System reicht es, den ganzzahligen Anteil – siehe
Abschnitt 1.2.1 – und den gebrochenen Anteil unabhängig voneinander zu konvertieren.
Beispielsweise folgt mit β = 2 aus
0.625 = (+1) · (. b1 b2 b3 . . .)2 mit bk ∈ {0, 1}
−1 −2 −3
= b1 2 + b2 2 + b3 2 + ...
Beispiel 1.12.
Algorithmus 1.11 führt mit a = 0.1 auf die folgenden Größen für das Dualsystem:
b1 = int (0.2) = 0 , c1 = 0.2
b2 = int (0.4) = 0 , c2 = 0.4
b3 = int (0.8) = 0 , c3 = 0.8
b4 = int (1.6) = 1 , c4 = 0.6
b5 = int (1.2) = 1 , c5 = 0.2
b6 = int (0.4) = 0 , c6 = 0.4
... ...
Man erkennt, dass sich die Dualziffernfolge periodisch wiederholt:
Beispiel 1.14.
Für die letzte Zahl ergibt sich im Dualsystem a = (0.00011)2 , eine Zahl mit unendlich
vielen Dualstellen und somit unendlich vielen tragenden Ziffern.
für ein k ∈ ZZ dargestellt werden, wobei v ∈ {+1, −1} das Vorzeichen von a angibt
und d1 = 0 gilt. (1.4) heißt normalisierte β-Gleitpunktdarstellung von a,
m = (. d1 d2 d3 . . .)β ihre β-Mantisse, und k ihr β-Exponent . Besitzt die Mantisse
s tragende Ziffern, s ∈ IN, so heißt sie s-stellig.
Einem Computer stehen für Berechnungen nur endlich viele in ihm darstellbare Zahlen,
die Maschinenzahlen, zur Verfügung. Die Mantissen m dieser Maschinenzahlen
haben gewöhnlich eine feste Anzahl von Ziffern. Ferner ist der Exponent k ∈ ZZ durch
−k1 ≤ k ≤ k2 mit k1 , k2 ∈ IN begrenzt. Eine weltweit gültige Norm nach ANSI (Ame-
rican National Standards Institute) und IEEE (Institute of Electrical and Electronics
Engineers) schreibt für reelle Zahlen, als Dualbruch mit 4 Bytes (= 32 Bits) im Rechner
hinterlegt, das Format
vor, wobei
→µ das Vorzeichen der Zahl wiedergibt,
→ e 1 . . . e8 zur Darstellung des Exponenten k in (1.4) gehört; als nicht negative ganze
Zahl, auch Charakteristik genannt, wird hiervon stets die Konstante 127
abgezogen, d. h. k = (e1 . . . e8 )2 − 127,
→ d2 . . . d24 die Mantisse (bei negativen Zahlen wieder komplementär genommen) der
Zahl darstellt, wobei d1 = 1 zu ergänzen ist (als normalisierte Gleitpunkt-
darstellung bleibt für d1 = 0 im Dualsystem nur d1 = 1 übrig!).
Hinzu kommen zusätzliche Kennungen für Fehlerfälle wie ”Overflow” (etwa bei Division
durch 0) und anderes mehr.
Für die Darstellung von 0.1 = (+1) · (. 00011)2 = (+1) · (. 1100)2 · 2−3 bedeutet das
. .
0 .. 01111100 .. 100 1100 1100 1100 1100 1100 1100 11 . . .
| |
=124 d2 d23
Als 2-Exponent ergibt sich 124−127 = −3, und d1 = 1 ist unterdrückt worden. Mithin ist
klar, dass 0.1 im Rechner nicht exakt wiedergegeben werden kann! Rundet die Rechner-
arithmetik auf, so wird d24 zu 1 gewählt, und die interne Maschinenzahl ist größer als 0.1
(genauer: 0.1000000015 . . .); wird abgerundet, also d24 = 0 gesetzt, dann hat die interne
Maschinenzahl einen Wert unterhalb von 0.1 (genauer: 0.09999999404 . . .). Deshalb wird
ein Rechner – übrigens unabhängig vom Speicherformat – für das Produkt 10 · 0.1 nie
genau den Wert 1 ermitteln.
Dass mit Gleitpunktzahlen auf Rechenanlagen überhaupt gearbeitet wird, liegt daran,
dass man damit bei gleicher Speichergröße einen enorm viel größeren Zahlenbereich über-
streicht.
Beispiel 1.16.
Nimmt man – abgesehen vom Vorzeichen – an, dass 8 dezimale Speichereinheiten zur
Verfügung stehen, so ließen sich damit in den einzelnen Systemen die folgenden positiven
Zahlen darstellen:
1) Ganzzahlsystem
− kleinste darstellbare positive Zahl 0000 0001
− größte darstellbare positive Zahl 9999 9999
Da alle Versionen von demselben Speicheraufwand ausgehen, lassen sich jeweils gleich
viele verschiedene Maschinenzahlen darstellen. Die zugehörigen Zahlenwerte sind aller-
dings auf dem reellen Zahlenstrahl ungleich verteilt:
• Bei dem Ganzzahlsystem lassen sich im positiven Bereich alle ganzen Zahlen zwi-
schen 1 und 99 999 999 (=100 Mio. −1) hinterlegen.
• Bei dem angegebenen Festpunktsystem lassen sich positive Zahlen zwischen 0.0001
(= 1 Zehntausendstel) und 9999.9999 (= 1 Zehntausendstel unter 10 000) wiederum
gleichabständig, diesmal im Abstand von 1 Zehntausendstel, als Maschinenzahlen
abspeichern.
• In dem betrachteten normaliserten Gleitpunktsystem liegen zwischen jeder zulässi-
gen Zehnerpotenz – also zwischen 10−50 und 10−49 genauso wie zwischen 10−49 und
10−48 usw. – gleich viele Maschinenzahlen (genauer: jeweils 899 999). Das heißt das
Intervall zwischen 0.1 und 1 enthält genauso viele darstellbare Zahlen wie das In-
tervall zwischen 1 und 10 oder das Intervall zwischen 10 und 100. Mithin sind diese
Zahlen höchst uneinheitlich dicht und unregelmäßig verteilt: Grob gesprochen lie-
gen zwischen 0 und 1 genauso viele Gleitpunktzahlen wie oberhalb von 1.
Daraus folgt:
Der absolute Fehler |a −A| gibt Auskunft über die Anzahl gültiger Dezimalen.
12 1. Darstellung von Zahlen und Fehleranalyse
Beispiel 1.18.
Sei a = 180.1234567. Gerundet auf 4 Dezimalen erhält man A = 180.1235, so dass gilt
|a − A| = 0.0000433 ≤ 0.5 · 10−4 ;
damit besitzt A also 4 gültige Dezimalen.
Beispiel 1.19.
Jede der folgenden Zahlen ist auf vier tragende Ziffern (d. h. auf eine dezimale Gleit-
punktzahl mit 4-stelliger Mantisse) nach Definition 1.17 zu runden und als normalisierte
dezimale Gleitpunktzahl (Definition 1.15) darzustellen.
Beispiel 1.21.
Jede der folgenden Zahlen ist nach Definition 1.20 auf vier tragende Ziffern (4-stellige
Mantisse) zu runden:
3.2355 ⇒ 3.236 ,
3.2345 ⇒ 3.234 ,
3.234500 ⇒ 3.234 ,
3.234501 ⇒ 3.235 .
Die letzte Ziffer einer Näherungszahl sollte immer eine gültige (sichere) Ziffer sein.
1.3 Rechnung mit endlicher Stellenzahl 13
Beispiel 1.23.
1. Sei X eine Näherungszahl von x. Wird keine Angabe über den absoluten Fehler
|∆x | von X gemacht, so sollte die letzte tragende Ziffer der Zahl eine sichere Ziffer
sein. Die Schreibweise X = 3.14 sollte bedeuten, dass für den absoluten Fehler von
X gilt |∆x | ≤ 0.5 · 10−2 , d. h. dass X zwei sichere Dezimalen bzw. drei sichere
Ziffern besitzt.
2. Die Angabe X = 0.004534 ± 0.000004 bedeutet, dass X wegen |x − X| ≤ 0.5 · 10−5
5 sichere Dezimalen und 3 sichere Ziffern besitzt. Dagegen heißt X = 0.004534 ±
0.000006, dass X wegen |∆x | > 0.5 · 10−5 , aber |∆x | ≤ 0.5 · 10−4 nur vier sichere
Dezimalen und zwei sichere Ziffern hat.
Die Anzahl sicherer Dezimalen liefert somit eine Abschätzung des absoluten Fehlers.
Die Anzahl der sicheren Ziffern hingegen liefert eine grobe Schätzung des relativen
Fehlers.
Es gilt der
Satz 1.24.
a habe die Darstellung (1.4) im Dezimalsystem (β = 10). A sei eine s-stellige dezimale
Gleitpunktzahl, die aus a durch Rundung auf eine s-stellige Mantisse entstanden ist.
Dann gelten folgende Abschätzungen:
Beweis.
Sei a in der normalisierten Darstellung (1.4) gegeben
a = v · (. d1 d2 . . . ds ds+1 . . .) · 10k , v ∈ {+1, −1}
und werde auf s Mantissenstellen gerundet
A = v · (. d1 d2 . . . Ds ) · 10k .
Dann gilt
1 1
|A − a| ≤ · 10−s · 10k = 10k−s
2 2
und weiter
A − a 0.5 · 10k−s 0.5 · 10−s 0.5 · 10−s
= 5 · 10−s .
a ≤ (0.d1 d2 . . .) 10k = .d1 d2 . . . ≤ 0.1
Daraus folgt:
Der relative Fehler gibt Auskunft über die Anzahl gültiger Ziffern einer Zahl.
14 1. Darstellung von Zahlen und Fehleranalyse
Beispiel 1.25.
Gegeben ist die Zahl a = 180.1234567, die in normalisierter dezimaler Gleitpunktdarstel-
lung lautet
a = .1801234567 · 103 (also k = 3) .
Bei Rundung auf 6-stellige Mantisse (s = 6) erhält man
A = .180123 · 103 .
Beispiel 1.26.
Es sollen die Zahlen a = 0.054 320 69 und b = 999.964 88 in einer 5-stelligen dezimalen
Gleitpunktarithmetik addiert werden. Die Zahlen a und b werden zunächst im gültigen
Format dargestellt:
A = gl (a) = 0.54321 · 10−1
B = gl (b) = 0.99996 · 103
gl (·) bezeichne hierbei die zulässige Gleitpunkthinterlegung und -rechnung. Damit ergibt
sich .
A 0.00005 ... 4321 ·103
+ B .
0.99996 . ·103
..
1.00001 .. ·103
.
= 0.10000 . 1 ·104 (normalisiert)
Exakt wären a + b = 1 000.019 200 69 und A + B = 1 000.014 321. Nach Satz 1.24 liegen
die absoluten und relativen Fehler zu Beginn bei
|∆a | = |a − A| ≤ 0.5 · 10−6 , |δa | = ∆aa ≤ 5 · 10−5 ,
∆
|∆b | = |b − B| ≤ 0.5 · 10−2 , |δb | = b ≤ 5 · 10−5 .
b
1.3 Rechnung mit endlicher Stellenzahl 15
Bezeichnet x = a + b die exakte Summe und X deren Näherung in der benutzten Arith-
metik, so erhält man entsprechend:
∆x
|∆x | = |x − X| = 1.92 . . . · 10−2 , |δx | = = 1.92 · 10−5 .
x
Obwohl das Resultat X dem auf die benutzte Gleitpunktarithmetik gerundeten exakten
Ergebnis gl (x) entspricht, hat sich der absolute Fehler in Relation zu den Eingangsfehlern
doch merklich vergrößert.
Beispiel 1.27.
Soll das Produkt der Zahlen a = 0.030 121 48 und b = 109.9761 in einer 5-stelligen dezi-
malen Gleitpunktarithmetik berechnet werden, so ergibt sich mit den gerundeten Werten
A =gl (a) und B =gl (b) über
.
A 0.30121 ... ·10−1
· B .
0.10998 . ·103
.
0.03312 ... 707 . . . ·102
= 0.33127 .. 07 . . . ·101 (normalisiert)
Exakt hätten sich a · b = 3.312 707 58 . . . und A · B = 3.312 707 58 . . . ergeben; wiederum
entspricht das Resultat X = gl (A · B) dem exakten Ergebnis x = a · b, gerundet auf die
verwendete Arithmetik: X = gl (x).
Man muss beachten, dass mathematische Gesetze wie das Assoziativ- und das Distri-
butivgesetz bei einer Rechnung in einer Gleitpunktarithmetik ihre Allgemeingültigkeit
verlieren. Beispielsweise ist mit 5-stelliger Mantissenrechnung
16 1. Darstellung von Zahlen und Fehleranalyse
gl (0.12345 · 10−2 + 0.22232 · 102 ) − 0.22223 · 102 = 0.10000 · 10−1 ,
hingegen
gl 0.12345 · 10−2 + (0.22232 · 102 − 0.22223 · 102 ) = 0.10235 · 10−1 .
Beispiel 1.29.
Es sei a = 1, die Darstellung erfolge mit 6-stelliger dezimaler Mantisse. Dann gilt
gl (1) = .100 000 · 101 ,
= 5 · 10−6 .
Denn: Bei der Addition 1 + in 6-stelliger Gleitpunktarithmetik ergibt sich
1 = .100 000 · 101
−6
+5 · 10 = .000 0005 · 101
1 + 5 · 10−6 = .100 001 · 101 = 1 .
Jede positive Zahl, die kleiner als 5 · 10−6 ist, wäre bei der Addition unberücksichtigt
geblieben.
Für einen Rechner, der reelle Zahlen a auf Gleitpunktzahlen A = gl (a) mit s-stelliger
Mantisse zur Basis β rundet, gilt
A − a β
−s
a ≤ 2 · β = ,
wobei die Maschinengenauigkeit nach Definition 1.28 bedeutet. Im gewohnten Dezi-
malsystem führt das auf = 5 · 10−s . Für Rechner mit Dualzahlarithmetik ergibt sich
demnach = 2−s , also für den zuvor angegebenen ANSI- bzw. IEEE-Standard mit s = 24
dann = 2−24 ≈ 5.96 · 10−8, was nach Satz 1.24 über den relativen Fehler mit einer etwa
7-stelligen dezimalen Genauigkeit gleichzusetzen ist.
Kombinierter Test
AbsErr und RelErr seien Schranken für den absoluten bzw. relativen Fehler. Für den
Einsatz in Programmen ist der folgende kombinierte Test zweckmäßig, der wahlweise
eine Abfrage auf den absoluten oder den relativen Fehler erlaubt:
a) RelErr = 0 und AbsErr > 0
|a − A | ≤ |A| RelErr + AbsErr mit (1.5)
b) RelErr > 0 und AbsErr = 0
1.4 Fehlerquellen 17
(1.5) ist mit a) eine Abfrage auf den absoluten Fehler und mit b) eine Abfrage auf den
relativen Fehler.
Es macht unter Umständen Sinn, zugleich beide Fehlerschranken verschieden von Null
zu wählen, etwa dann, wenn die Größenordnung des Ergebnisses nicht bekannt ist und
in Abhängigkeit davon eine gewisse Genauigkeit erreicht werden soll (bei betragsmäßig
großen Zahlen überwiegt so die relative Fehlerschranke, bei betragsmäßig kleinen die
absolute Fehlerschranke).
1.4 Fehlerquellen
Bei der numerischen Behandlung eines Problems treten verschiedene Fehlerquellen auf.
Der gesamte Fehler einer Berechnung von der Eingabe bis zur Ausgabe setzt sich im
Allgemeinen zusammen aus:
• Eingabefehlern (Eingangsfehlern)
• Verfahrensfehlern (Abbruchfehlern, Diskretisierungsfehlern)
• Fortpflanzungsfehlern
• Rechnungsfehlern
1.4.1 Eingabefehler
Eingabefehler sind die Fehler, mit denen die Eingabedaten behaftet sind, z. B. wegen
fehlerhafter Messungen oder Rundung. Diese unvermeidbaren Fehler wirken sich auf die
Ausgabedaten eines Algorithmus aus. Daher müssen numerische Algorithmen so konzi-
piert werden, dass der Einfluss von Eingabefehlern möglichst begrenzt wird (siehe auch
Abschnitt 1.4.4).
Beispiel 1.30.
Für das Integral
1
In = (1 − x)n · sin x dx
0
In keinem Rechner wird der Wert cos(1) = 0.540 302 3 . . . exakt darstellbar sein. Ver-
schiedene Genauigkeiten bei diesem Eingabedatum haben schon nach wenigen Schritten
erstaunliche Konsequenzen:
1.4.2 Verfahrensfehler
Verfahrensfehler sind solche, die durch die verwendete numerische Methode verursacht
werden. Das Ersatzproblem für eine zu lösende Aufgabe muss so formuliert werden, dass
• es numerisch gelöst werden kann und
• seine Lösung nicht wesentlich von derjenigen des gegebenen Problems abweicht.
Bei einem geeignet formulierten Ersatzproblem wird dessen Lösung eine Näherungslösung
für das gegebene Problem sein. Die Differenz zwischen diesen beiden Lösungen heißt der
Verfahrensfehler. Dieser Verfahrensfehler hängt in hohem Maße vom gegebenen Problem
und von dem ausgewählten Ersatzproblem ab. Der Verfahrensfehler berücksichtigt weder
Eingabe- noch Rechnungsfehler im Verfahren selbst.
Beispiele:
1. Die Berechnung eines bestimmten Integrals wird ersetzt durch die Berechnung einer
endlichen Summe. Dann ist der Verfahrensfehler die Differenz
b
n
f (x) dx − A f (x
k .
)
k
k=0
a
2. Die Berechnung der 1. Ableitung einer Funktion f wird ersetzt durch die Berech-
nung des vorderen Differenzenquotienten. Dann ist der Verfahrensfehler
f (x) − f (x + h) − f (x) .
h
1.4 Fehlerquellen 19
Fehler der Ausgabedaten eines Problems, die durch Fehler der Eingabedaten erzeugt wer-
den, heißen Fortpflanzungsfehler. So lassen sich die Ergebnisse in Beispiel 1.30 dadurch
erklären, dass sich ein Eingabefehler ε0 , wenn statt I0 ein Wert I0 mit | I0 − I0 | ≤ ε0
genommen wird, fortlaufend vervielfacht: Statt I2 = 1 − 2 · 1 · I0 berechnet man dann
I2 = I2 −2 ε0 , anstelle von I4 = 1−4·3·I2 dann I4 = 1−4·3· I2 = I4 +24 ε0 , und entspre-
chend folgt I6 = I6 −720 ε0 sowie schließlich I8 = I8 +40 320 ε0 für die berechneten Werte.
Satz 1.31.
Es sei f : G ⊂ IRn → IR eine auf einem Gebiet G stetig differenzierbare Funktion,
ferner seien x = (x1 , . . . , xn )T ∈ G und X = (X1 , . . . , Xn )T ∈ G.
Dann gibt es ein x = (x1 , . . . , xn )T ∈ G mit xi zwischen xi und Xi für i = 1(1)n, so
dass für den Fortpflanzungsfehler gilt:
n
∂f (x )
∆y = y − Y = f (x ) − f (X ) = ∆xi ,
i=1
∂xi
wobei ∆xi = xi − Xi der Eingabefehler von xi ist.
Praktisch werden die Eingabefehler nicht exakt bekannt sein, sondern nur obere Schran-
ken εxi mit | xi − Xi | ≤ εxi dafür. Ebenso lässt sich, abgesehen von Ausnahmefällen, die
Zwischenstelle“ x nicht genau angeben. Deshalb wendet man die Formel aus Satz 1.31
”
gewöhnlich wie folgt an: n
∂f (X )
| ∆y | ≤
∂xi · εxi (1.6)
i=1
Für kleine Eingabefehler εxi liegt auch x nahe bei X , und wegen der Stetigkeit der
partiellen Ableitungen von f bleibt dann die in (1.6) auf der rechten Seite erfolgte In-
korrektheit in überschaubaren Grenzen.
Das Ergebnis aus Satz 1.31 lässt sich leicht auf den relativen Fehler übertragen:
∆y n
xi ∂f (x )
δy = = · · δx i , (1.7)
y i=1
f (x ) ∂xi
20 1. Darstellung von Zahlen und Fehleranalyse
∆xi xi − Xi
wobei δxi = = der relative Eingabefehler von xi ist.
xi xi
Für die praktische Anwendung heißt das wiederum
n
∆y Xi ∂f (X )
≤ ·
y f (X ) ∂xi · xi
i=1
xi − Xi
mit den oberen Schranken ≤ xi für die relativen Eingabefehler.
xi
Beispiel 1.32.
Man gebe eine Abschätzung des Eingangsfehlers an, der bei der Bestimmung des spezifi-
schen Gewichtes γ eines Zylinders vom Gewicht G, dem Radius r und der Höhe h entsteht,
wenn ∆G, ∆r die wahren Messfehler von G, r sind und ∆π der wahre Rundungsfehler
für π ist; h sei genau gemessen.
Lösung: Es gilt
G G
γ= = 2 = γ(G, r, π, h) .
V πr h
γ(G, r, π, h) = γ ∗ (G, r, π) .
Mit
f (x1 , x2 , x3 ) = γ ∗ ,
x1 = G, x2 = r, x3 = π,
Bei kleinen Radien wirken sich also alle Messfehler stärker aus als bei großen Radien.
als Verstärkungsfaktoren der relativen Eingabefehler δxi anzusehen. Man erkennt, dass
sich Funktionen f mit betragsmäßig kleinen partiellen Ableitungen fxi bezüglich der
1.4 Fehlerquellen 21
Für große Probleme ist eine Abschätzung von Fortpflanzungsfehlern und Konditionszah-
len sehr kompliziert und selten möglich. In solchen Fällen ist der Einsatz statistischer
Fehlerabschätzungen sinnvoll, siehe [HENR1972], Bd.2, S.381.
Es sei
y = f (x1 , x2 ) = x1 + x2 .
Wenn die Eingabedaten x1 und x2 dasselbe Vorzeichen haben, ist die Addition eine
gut konditionierte Operation (wegen |Ki | = |xi /(x1 + x2 )| < 1). Wenn jedoch x1
und x2 verschiedene Vorzeichen haben und dem Betrage nach nahezu gleich sind,
ist |x1 + x2 | sehr klein und somit |Ki | = |xi |/|x1 + x2 | sehr groß; die Addition ist
dann schlecht konditioniert. In diesem Fall gehen bei der Addition tragende Ziffern
verloren, ein Effekt, den man Auslöschung nennt.
Beispiel 1.33.
Gegeben seien die drei Werte
x1 = 123.454 Mrd. C,
x2 = 123.446 Mrd. C und
x3 = 123.435 Mrd. C.
Die jeweils in der fünften tragenden Stelle liegenden Rundefehler nehmen schon
nach einer Rechenoperation Einfluss auf die Größenordnung des Ergebnisses!
Der Ausdruck x − (x − y) entspricht korrekt dem Wert y. Mit x = 1030 und y = 106
würde allerdings ein Rechner stattdessen das Ergebnis 0 produzieren, da in der
vorgegebenen Auswertungsreihenfolge die Gleitpunktarithmetik den Klammeraus-
druck zu x werden lässt. In bestimmten Fällen lassen sich durch geschickte Umfor-
mungen derartige Gegebenheiten vermeiden.
Beispiel 1.34.
x
Es soll der Ausdruck sin (1+x)−sin (1) für x nahe 0 auswertet werden. Setzt man in
(Taschen-)Rechnern für x Werte wie 10−5 , 10−10 , 10−15 , 10−20 usw. ein, so wird
man schon bald eine Fehlermeldung erhalten ( Division durch Null“). Dies liegt
”
natürlich daran, dass die beiden Werte im Nenner rechnerintern nicht mehr unter-
schieden werden können und deshalb die Differenz auf eine 0 führt.
Die Auswertung des umgeformten Ausdrucks erweist sich als numerisch stabil: Für
auf 0 zugehende Werte von x liefert auch ein (Taschen-)Rechner Werte in der kor-
rekten Größenordnung von
x 1
lim = = 1.850 815 7 . . .
x→0 sin (1 + x) − sin (1) cos (1)
1.4 Fehlerquellen 23
Beispiel 1.35.
Die Funktion
f (x, y) = 9 x4 − y 4 + 2 y 2
Die letzten Beispiele lassen evident erscheinen, dass man sich keinesfalls auf je-
des Rechnerergebnis verlassen darf; man sollte stets eine Plausibilitätsprüfung
vornehmen! Selbst wenn das Resultat vertrauenswürdig ist, kann man nie allen an-
gegebenen Stellen Glauben schenken.
Es sei
y = f (x1 , x2 ) = x1 x2 .
Es sei x1
y = f (x1 , x2 ) = .
x2
Hier sind fx1 = 1/x2 und fx2 = −x1 /x22 , und somit gilt (wegen K1 ≈ 1 und
K2 ≈ −1)
x1 1 x2 −x1
δy = · · δx 1 + · · δx 2 ≈ δx 1 − δx 2 .
x1 /x2 x2 x1 /x2 x22
Also ist auch die Division gut konditioniert.
24 1. Darstellung von Zahlen und Fehleranalyse
Es sei
y = f (x1 ) = xp1 für p > 0 , x1 > 0 .
Hier ist fx1 = px1p−1 und (wegen K1 ≈ p) δy ≈ pδx1 ; damit sind Wurzeln im All-
gemeinen gut konditioniert und Potenzen mäßig schlecht konditioniert für große p.
Dies ist einer der Gründe dafür, dass man keine vernünftigen Resultate erwarten
kann, wenn versucht wird, ein Polynom
Pn (x) = a0 + a1 x + a2 x2 + . . . + an xn
auszuwerten, indem Terme der Form ak xk summiert werden; man sollte stattdes-
sen immer das besser konditionierte Horner-Schema benutzen, in dem abwechselnd
multipliziert und addiert, aber nicht potenziert wird (vgl. Abschnitt 3.2).
Zur Durchführung eines numerischen Verfahrens muss ein Algorithmus formuliert werden.
In diesem Sinne wurden bereits zuvor die Algorithmen“ 1.8, 1.9 und 1.11 formuliert.
”
Während der Ausführung der Rechenoperationen ergibt sich durch Anhäufung lokaler
Rechnungsfehler ein akkumulierter Rechnungsfehler. Die√ lokalen Rechnungsfehler entste-
hen z. B. dadurch, dass irrationale Zahlen wie π, e, 2 durch endliche Dezimalbrüche
(Maschinenzahlen) ersetzt werden; dadurch werden Abbruch- oder Rundungsfehler er-
zeugt. Hinreichend kleine Größen, die durch Unterlauf entstehen, werden vernachlässigt.
Führende genaue Stellen können bei der Subtraktion fast gleich großer Zahlen ausgelöscht
werden. Mit der Anzahl der Operationen in einem Algorithmus wächst somit die Gefahr,
dass völlig falsche Ergebnisse entstehen.
Algorithmen, die eine Verstärkung und Anhäufung von Rundungsfehlern vermeiden, wer-
den numerisch stabil genannt. Es gibt unterschiedliche Definitionen für den Begriff der
numerischen Stabilität:
3. Eine sehr häufig benutzte Definition für numerische Stabilität fordert die Existenz
eines x∗ nahe bei x, so dass f ∗ (x∗ ) = f (x) ist, vgl. [WILK1969].
Beispiel 1.37.
Hat man eine Gleichung der Form
yn+1 = a yn + b yn−1 (n ∈ IN) ,
auch Differenzengleichung genannt, mit reellen Konstanten a, b und n ∈ IN, so lässt sich
jeder Wert yn hierüber nach Angabe von zwei konkreten Startwerten“ y0 und y1 be-
”
rechnen. Eine allgemeine Lösung der Differenzengleichung kann mit Hilfe des Ansatzes
yn = ξ n (1.8)
ξ n+1 = a ξ n + b ξ n−1
erfüllt ist. Dann ist ein über obige Differenzengleichung gegebener Algorithmus numerisch
stabil.
Die Verwendung instabiler Algorithmen ist praktisch sinnlos. Aber selbst dann, wenn
ein stabiler Algorithmus verwendet wird, hat es natürlich keinen Sinn, zwar mit einer
exakten Prozedur zu arbeiten, aber die Berechnung mit nur geringer Genauigkeit aus-
zuführen; es werden dann natürlich große Rechnungsfehler auftreten.
26 1. Darstellung von Zahlen und Fehleranalyse
Gewöhnlich sind numerisch stabile Algorithmen und gut konditionierte Probleme notwen-
dig, um überhaupt zufriedenstellende Resultate erreichen zu können; diese Bedingung ist
jedoch keineswegs hinreichend, da das Instrument Computer“, mit dem das Ergebnis
”
erzeugt wird, mit begrenztem Speicherplatz und begrenzter Zeit arbeitet. Computer-
Operationen setzen sich zusammen aus dem Schreiben in den Speicher, dem Lesen aus
dem Speicher, Overhead (z. B. zusätzlicher Verwaltungsaufwand) und den arithmetischen
Operationen.
f (x) = 0 , (2.1)
n
f (x) ≡ Pn (x) = aj xj , aj ∈ IR , an = 0 , n ∈ IN
j=0
ist, heißt die Gleichung (2.1) algebraisch, und die natürliche Zahl n heißt der Grad des
Polynoms bzw. der algebraischen Gleichung. Jede Gleichung (2.1), die nicht algebraisch
ist, heißt transzendent (z. B. ln x − 1/x = 0; x − sin x = 0).
In diesem Kapitel werden Verfahren zur Bestimmung einfacher und mehrfacher Nullstel-
len ξ ∈ I von f vorgestellt. Dabei wird zwischen den klassischen Iterationsverfahren (all-
gemeines Iterationsverfahren, Newton-Verfahren, Sekantenverfahren) und den sogenann-
ten Einschlussverfahren (Bisektion, Pegasus-Verfahren, Verfahren von Anderson-Björck,
Verfahren von King) unterschieden. Die angegebenen Einschlussverfahren benötigen zwei
Startwerte, in denen die Funktion f unterschiedliche Vorzeichen hat. Die Startwerte
schließen dann (mindestens) eine Nullstelle ungerader Ordnung von f ein. Dieser Ein-
schluss bleibt im Laufe der Rechnung erhalten. Diese Verfahren sind jedoch grundsätzlich
unter Anwendung von Satz 2.3 auch im Falle von Nullstellen gerader Ordnung anwend-
bar. Einschlussverfahren höherer Konvergenzordnung sind im Allgemeinen den klassi-
schen Iterationsverfahren vorzuziehen (vgl. Abschnitt 2.9). Verfahren zur Bestimmung
sämtlicher Nullstellen algebraischer Polynome ohne Kenntnis von Startwerten werden im
Kapitel 3 behandelt.
28 2. Lösung nichtlinearer Gleichungen
auf. Bei der Berechnung des Druckverlustes in einer Rohrströmung muss man den Rohr-
reibungsbeiwert λ für hydraulisch glatte Rohre bei turbulenter Strömung nach dem
universellen Prandtlschen Widerstandsgesetz bei vorgegebener Reynoldszahl Re aus der
transzendenten Gleichung
1 √
f (λ) = √ − 2 lg Re λ + 0.8 = 0
λ
ermitteln. Ebenfalls auf eine transzendente Gleichung stößt man bei der Betrachtung
der Ausstrahlung eines vollkommenen schwarzen Körpers“. Wird ein schwarzer Körper
”
erhitzt, so sendet er elektromagnetische Wellen verschiedener Wellenlänge λ (Wärme-
strahlen, Licht usw.) aus.
.... ...
.. .. ...
... ....... .
.... ..
................
.. λ 2 .... λ2 ......
.
λ4 ..
... ..............
...
...
..........................
....
...
.
........
....... ....
. λ3 ......
.... ................
..
λ 2
... ... ... ..
... .............. ... ....... ... .......... ..
.
......... ... .... .. ... .............
... .. .. ..
...... ..... ..... .....
... ......
.......... ..........
........ .............. ............... .............. ....... ........................................................................................
....... ......
......
λ ..... λ1
3 ..
..
.....................
........... .
........................
.....
...
.... ...... ......... .
...
...
..... .. ... ... ....... ...
..
...
...... ... ....
.
.
...... ...
... ............ ....
............ .....
...... .. .. .....
...
..
.. .... ... .
... .... ...
....
..............
... λ 2 .... .... .....
.... λ 2
... ..................... .......... ........... ......
. .
.................................................. .....
......... ..
.
.
.
...
....
....
....
Bei einer bestimmten Temperatur T liegt das Maximum der emittierten Strahlung E(λ, T )
bei λmax . Mit steigender Temperatur verschiebt sich die Stelle maximaler Emission nach
kürzeren Wellenlängen; es gilt (Wiensches Verschiebungsgesetz)
α
λmax (T ) =
z·T
mit α = 14.3881 · 10−3 m K. Die Konstante z ist die Lösung der transzendenten Gleichung
1 1
e−z = 1 − z bzw. f (z) = e−z − 1 + z = 0 , z ∈ IR , z = 0.
5 5
Mit dem so ermittelten z ≈ 4.9651 lässt sich die Formel für λmax (T ) angeben. Aus dieser
Formel kann dann z. B. ein Näherungswert für die absolute Temperatur der Sonnenober-
2.2 Definitionen und Sätze über Nullstellen 29
fläche bestimmt werden, wobei für das Sonnenspektrum nach Langley das Maximum der
Strahlung bei λmax = 5 · 10−7 m liegt.
Definition 2.1.
Eine Nullstelle ξ einer Funktion f ∈ C[a, b] heißt j-fache Nullstelle oder Nullstelle der
Ordnung j (j ∈ IN), falls f sich auf [a, b] in der Form
f (x) = (x − ξ)j h(x)
darstellen lässt mit einer stetigen Funktion h, für die h(ξ) = 0 ist.
Im Fall j = 1 heißt ξ einfache Nullstelle, für j ≥ 2 mehrfache Nullstelle.
Satz 2.2.
Die Funktion f sei im Intervall I j-mal stetig differenzierbar (j ∈ IN). Dann ist ξ ∈ I
genau dann eine j-fache Nullstelle von f , wenn gilt
f (ξ) = f (ξ) = . . . = f (j−1) (ξ) = 0 , f (j) (ξ) = 0 .
Mit der Aussage des folgenden Satzes 2.3 lassen sich Einschlussverfahren auch zur Be-
rechnung von mehrfachen Nullstellen einsetzen, insbesondere auch von Nullstellen gerader
Ordnung (also ohne Vorzeichenwechsel).
Satz 2.3.
Ist ξ ∈ I eine j-fache Nullstelle von f , j ≥ 2, und ist f (j+1)-mal stetig differenzierbar,
so ist ξ eine einfache Nullstelle von g mit
f (x)
g(x) = .
f (x)
30 2. Lösung nichtlinearer Gleichungen
Beweis. Nach Satz 2.2 gelten für j ≥ 2 f (ξ) = f (ξ) = . . . = f (j−1) (ξ) = 0 und
f (j) (ξ) = 0. Zu zeigen ist, dass für
f (x) f (x) f (x)
g(x) = und g (x) = 1 −
f (x) f 2 (x)
nach Satz 2.2 gelten
g(ξ) = 0 und g (ξ) = 0 .
Für den Nachweis werden die Taylorentwicklungen von f , f und f an der Stelle ξ
benötigt (sie finden auch in Abschnitt 2.5.3 Verwendung).
(x − ξ)j (j)
f (x) = f (ξ) + O (x − ξ)j+1 = (x − ξ)j h0 (x) (2.2)
j!
1
mit h0 (x) = f (j) (ξ) + O(x − ξ) und
j!
1 (j)
h0 (ξ) = f (ξ) = 0 ; (2.3)
j!
(x − ξ)j−1 (j)
f (x) = f (ξ) + O (x − ξ)j = (x − ξ)j−1 h1 (x) (2.4)
(j − 1)!
1
mit h1 (x) = f (j) (ξ) + O(x − ξ) und
(j − 1)!
1
h1 (ξ) = f (j) (ξ) = 0 ; (2.5)
(j − 1)!
(x − ξ)j−2 (j)
f (x) = f (ξ) + O (x − ξ)j−1 = (x − ξ)j−2 h2 (x) (2.6)
(j − 2)!
1
mit h2 (x) = f (j) (ξ) + O(x − ξ) und
(j − 2)!
1
h2 (ξ) = f (j) (ξ) = 0 . (2.7)
(j − 2)!
x = ϕ(x) (2.9)
eine in einem abgeschlossenen Intervall I stetige Funktion, und ξ ∈ I heißt eine Lösung
von (2.9) bzw. ein Fixpunkt der Abbildung ϕ, wenn ϕ(ξ) = ξ ist; darum heißt (2.9) auch
Fixpunktgleichung.
Die Untersuchung von Gleichungen der Form x = ϕ(x) bedeutet keine Beschränkung der
Allgemeinheit, denn es gilt der
Hilfssatz 2.5.
Sind f und g stetige Funktionen in einem abgeschlossenen Intervall I und ist g(x) = 0
für alle x ∈ I, dann besitzen die Gleichungen f (x) = 0 und x = ϕ(x) mit
Beweis. Ist ξ ∈ I Lösung von f (x) = 0, so folgt wegen f (ξ) = 0 aus (2.10) ϕ(ξ) = ξ. Ist
umgekehrt ξ ∈ I Lösung von x = ϕ(x), so folgt wegen ξ = ϕ(ξ) aus (2.10) f (ξ)g(ξ) = 0;
wegen g(ξ) = 0 ist also f (ξ) = 0 .
Jede geeignete Wahl von g liefert eine zu f (x) = 0 äquivalente Gleichung x = ϕ(x).
Häufig kann eine Gleichung f (x) = 0 auf die Form x = ϕ(x) gebracht werden, indem
irgendeine Auflösung nach x vorgenommen wird.
32 2. Lösung nichtlinearer Gleichungen
Beispiel 2.6.
f (x) = x2 + x − 2 = 0
(I) x = 2 − x2 = ϕ(x) ,
√
(II) x = 2 − x = ϕ(x) , x ≤ 2 ,
(III) x= 2
x − 1 = ϕ(x) , x = 0.
Bei der Angabe eines Intervalls I, in dem die Gleichungen äquivalent zur gegebenen sind,
müssen die Einschränkungen für x berücksichtigt werden.
Nun sei eine Gleichung der Form x = ϕ(x) mit dem zugehörigen Intervall I gegeben.
Dann konstruiert man mit Hilfe eines Startwertes x(0) ∈ I eine Zahlenfolge {x(ν) } nach
der Vorschrift
x(ν+1) := ϕ(x(ν) ) , ν = 0, 1, 2, . . . (2.11)
Diese Folge lässt sich nur dann sinnvoll konstruieren, wenn für ν = 0, 1, 2, . . .
x(ν+1) = ϕ(x(ν) ) ∈ I
ist, da ϕ nur für x ∈ I erklärt ist. Durch ϕ muss also eine Abbildung des Intervalls I in
sich gegeben sein, d. h. der Graph von y = ϕ(x) muss im Quadrat
Q = {(x, y)|x ∈ I, y ∈ I}
......................................................
....................................
....................................
..................
a
-
a b x
Abb. 2.2.
2.3 Allgemeines Iterationsverfahren 33
Wenn die Folge {x(ν) } konvergiert, d. h. wenn die Zahlen x(1) , x(2) , x(3) , . . . gegen ξ stre-
ben, somit
lim x(ν) = ξ
ν→∞
ist, dann ist ξ eine Lösung der Gleichung (2.9). Es gilt wegen der Stetigkeit von ϕ
Die Iterationsschritte (2.11) für ν = 0(1)N bilden zusammen mit dem Startwert x(0) das
algorithmische Schema des Iterationsverfahrens:
⎧ (0)
⎪
⎪ x = Startwert ,
⎪ (1)
⎪
⎪
⎪ x = ϕ(x(0) ) ,
⎪
⎪
⎨ x (2)
= ϕ(x(1) ) ,
· (2.12)
⎪
⎪
⎪
⎪ ·
⎪
⎪
⎪
⎪ ·
⎩ (N +1)
x = ϕ(x(N ) ) .
und es sei I = [−50, 0]. Das algorithmische Schema (2.12) des Iterationsverfahrens lautet
mit dem Startwert x(0) = −3:
x(0) = −3 ∈ I ,
x(1) = ϕ(x(0) ) = 2 − 32 = −7 ∈ I ,
x(2) = ϕ(x(1) ) = 2 − 72 = −47 ∈ I ,
x(3) = ϕ(x(2) ) = 2 − 472 = −2207 ∈ / I.
Der Verlauf der Rechnung zeigt, dass die so konstruierte Folge {x(ν) } nicht gegen die
Lösung ξ2 = −2 ∈ I konvergiert.
x(0) = −3 ∈ I ,
−3 − 1 = − 3 = −1.6666667 ∈ I ,
2 5
x(1) =
x(2) = −2.2000000 ∈ I ,
x(3) = −1.9090909 ∈ I ,
x(4) = −2.0476190 ∈ I .
Diese vier Iterationsschritte zeigen bereits, dass sich die so konstruierte Folge der Lösung
ξ2 = −2 immer mehr nähert, d. h. gegen die gesuchte Lösung konvergiert.
Es stellt sich also die Frage nach Bedingungen für die Konvergenz einer Iterationsfolge.
Die folgenden Aussagen über die Existenz einer Lösung der Gleichung (2.9) und deren
Eindeutigkeit dienen der Beantwortung dieser Frage. Es wird sich zeigen, dass die Be-
dingungen für die Existenz und Eindeutigkeit auch hinreichend für die Konvergenz der
mit (2.11) konstruierten Iterationsfolge sind.
Zum Nachweis der Existenz einer Lösung ξ ∈ I = [a, b] der Gleichung x = ϕ(x) folgen
nun lediglich geometrische Überlegungen; zum analytischen Nachweis s. [HENR1972]
Bd.1, S. 85/87. Man erhält ξ als Abszisse des Punktes, in dem die Graphen der
Gerade y = x und der Funktion y = ϕ(x) sich in dem Quadrat Q schneiden. Es wird
also davon ausgegangen, dass für x ∈ I auch ϕ(x) ∈ I ist. Wenn
a = ϕ(a) und
b = ϕ(b)
sind, müssen ϕ(a) > a und ϕ(b) < b sein. Die Punkte a, ϕ(a) und b, ϕ(b) liegen also
auf verschiedenen Seiten der Diagonale y = x im Quadrat Q. Ist ϕ stetig, so garantiert
der stetige Verlauf des Graphen von y = ϕ(x) in Q die Existenz von mindestens einem
Schnittpunkt mit der Gerade y = x Abb. 2.3(a) . Ist ϕ dagegen nicht stetig, so existiert
nicht notwendig ein solcher Schnittpunkt Abb. 2.3(b) .
a) y6 b) y6
y = ϕ(x)
b r
...................................................................
..............................................
b @ ..
........................................................................
@
.................................................
............................................... ..................................................
................................................................ ................................................
. . . . . . . . . . . . . . .. . . . . . . . . . . . . . . .... . . . . . . . . . . . . .
.............................................................. .............................................................
............................... ...............................
................................ y =x
................................. . . . . . . .
...............................
................................. . . . . . . .
y =x
............................................. ..............................................
...............................................
. ... ... ... ... .. ... ... .
.........................................
..............................................
. . . . .. . . . . . . . . . . . . . . . . . r
..............................................................................
y = ϕ(x)
..................................... r .. .
..................................................................................
PP y = ϕ(x)
.............................
.................................................. ........................................................................
................................................ ...............................................
.
r
................................................ . . . . . . . . . . . . . . . ... . . . . .
..................................................................... ......................................................................
..............................................
.................................................................. .............................................
. . . . . . . . . . . . . ... . . . . . . . .
a ......................
a
- -
a b a x b x
Abb. 2.3. (a) ϕ stetig; (b) ϕ nicht stetig
2.3 Allgemeines Iterationsverfahren 35
Anhand eines Beispiels wird gezeigt, dass die Stetigkeit (ii) für die Existenz einer Lösung
ξ ∈ I allein nicht hinreicht.
Beispiel 2.9.
Gegeben: Die Gleichung ex = 0, I = [−a, 0], a > 0, von der bekannt ist, dass sie keine
endliche Lösung besitzt.
x = ϕ(x) = x − ex .
Die Funktion ϕ ist stetig in [−a, 0]. Wegen ϕ(−a) = −a − 1/ea < −a ist ϕ(−a) ∈
/ I und
somit die Bedingung (i) des Satzes 2.8 nicht erfüllt.
Die Frage nach der Eindeutigkeit einer Lösung der Gleichung x = ϕ(x) kann man be-
antworten, wenn ϕ in I einer sogenannten Lipschitzbedingung genügt. Eine Funktion ϕ
heißt lipschitzbeschränkt, wenn es eine Konstante L mit 0 ≤ L < 1 gibt, so dass
gilt. L wird Lipschitzkonstante genannt, und (2.13) heißt eine Lipschitzbedingung für die
Funktion ϕ. Eine differenzierbare Funktion ϕ ist sicher lipschitzbeschränkt, wenn
Abbildungen ϕ, für die eine Lipschitzbedingung (2.13) bzw. (2.14) gilt, werden auch als
kontrahierende Abbildungen bezeichnet, weil der Abstand |ϕ(x)−ϕ(x )| der Bilder kleiner
ist als der Abstand |x − x | der Urbilder.
36 2. Lösung nichtlinearer Gleichungen
und daraus
(1 − L) |ξ1 − ξ2 | ≤ 0 .
Wegen 1 − L > 0 ist |ξ1 − ξ2 | ≤ 0, also kann nur |ξ1 − ξ2 | = 0 sein und damit ξ1 = ξ2 .
Die Gleichung x = ϕ(x) besitzt also höchstens eine Lösung ξ ∈ I.
Da eine Funktion ϕ, die in I einer Lipschitzbedingung genügt, überall in I stetig ist (die
Umkehrung gilt nicht, d. h. nicht jede stetige Funktion genügt einer Lipschitzbedingung)
und da die Stetigkeit von ϕ zusammen mit der Bedingung (i) des Satzes 2.8 hinreichend
ist für die Existenz einer Lösung ξ in I, gilt mit Satz 2.10 weiter der
Falls eine Lösung existieren würde, wäre sie nach Satz 2.10 eindeutig bestimmt. Im vor-
liegenden Falle existiert aber keine Lösung, da die Bedingung (i) in Satz 2.8 bzw. 2.11
nicht erfüllt ist.
2.3 Allgemeines Iterationsverfahren 37
Die Funktion ϕ : I → IR sei in I differenzierbar. Dann ergibt sich anschaulich, dass die
durch die Iterationsvorschrift x(ν+1) = ϕ(x(ν) ) definierte Folge {x(ν) } konvergiert, wenn
ϕ die Bedingungen des Satzes 2.11 für die Existenz und Eindeutigkeit einer Lösung der
Gleichung x = ϕ(x) erfüllt:
(i) ϕ(x) ∈ I für alle x ∈ I,
(ii) |ϕ (x)| ≤ L < 1 für alle x ∈ I.
Die Konvergenz ist für 0 ≤ ϕ (x) < 1 monoton (Abb. 2.4) und für −1 < ϕ (x) ≤ 0
alternierend (Abb. 2.5).
y 6 y 6
. ...
...
y =x ...
.................................................................................................................
.... - .. y =x
...........
..
........... 6
.. ....
.. ....
.. ....
..
..
..
.......... .. ... ..
..
...
............
. ..
...
.... ..
.
......... .. .... ..
........ ϕ(x) ....
........ .. ..
r
....
.......
- ...
.. ....
.
...
...
. ....... .. ... ..
.
-
...........................................
....... ... .. ..........
....... .. ..
-6
............. ...
...... ..
..
.. 6 - ... ..........................
......
.. ..
r?
.
..................... ...
..
.. ................. ..
.
...
..
66
. .
. ..
.
.
... . .. ..
...... .. .... .. ............................ ..
-6
.. .. ..
..
.
..... ..
......
.
........................
.
.
.... ..
..
..
.. ? ..
..
.. .. ..............
......................................
.
..
..
..
.. .
. ... .. ... ... .. ... ................... ..
. .
.
.
.... ..
... .. 6 .
.
. ..
.
.
...
...
..
.. .. . . . . ............
?
.. .
...............................................................................................................................
... ... .... .... .. .. .. .. .. . .. ..
...
. . .. .. . ..
.. ... ... ... ... .. ... ..... ... . ..... ϕ(x) ...
.. . . .. . .. .. .
.. . . ... .. .. .. .. . .. ..
.. ... ... .. .. .. .. .. .... .. ..
.. ... .. ..
.
.. .... .... .. .. .... .... .. .. .... ....
.. .. .
... ..
..
.
...
.
....
. ... - ..
.
....
.
.... ... ...
. . .
....
. - ....
.
ξ x ξ x
x(0) x(1) x(2) x(0) x(2) x(4) x(5) x(3) x(1)
Je kleiner der Betrag der 1. Ableitung der Schrittfunktion ϕ ist, desto schneller konver-
giert die Iterationsfolge gegen ξ.
38 2. Lösung nichtlinearer Gleichungen
Für |ϕ (x)| > 1 divergiert das Verfahren, wie aus den Abbildungen 2.6 und 2.7 zu erken-
nen ist.
y 6 ..
...
...
y 6 ...
...
.... ...
ϕ(x) ..
... y =x ...
...
ϕ(x) y =x
... ...
.... ...
...
.. ...
...
-
...
... ...................................................
.... ..... .
..
...
... r
......
.. ...6
.. ....
..
..
..
....
. .. .... ..
... .. .. ... ..
.............. ...
... ....... ..
... .. ..
..
...
r
...
....
..
..
..
6
.. .....
.. ... .. .. .. ..
.. .... .. .. ... ..... ..
.. .
.. ... .. .. ................................. ..
..... .. .. ... ..
?
. .
.... ... ......
...................................
.
...
..
.. ...
.. ..
.. .. ....
..
...
6 .
. .....
.
..
..
.. .
... .
. .... ...
.. .
. ... .. .. ... ... .... .
... ...
?
.. .. ... .. .. .... ..
.. .... . .. .. . . . ......
.. ..
. .
. .. .. .....................................................................................
.. .
... .
. .. .. .... ... . ... .
.... ........
.. ..... . .. .. ..
. . .....
... ..... ... .. .. .... .... .. ... .... .....
.....
.. . .. .. .. .....
? - -
. . ..
.. .... .. .. ... .... .. ... ...
ξ x ξ x
x(2) x(1)x(0) x(3) x(1) x(0) x(2)
Abb. 2.6. Divergenz, 1 < ϕ (x) Abb. 2.7. Divergenz, ϕ (x) < −1
y 6
...........................
......... ......
.....
..
....
...
.....
....
....
y =x
.
... ...
y =2 − x 2 ..... ...
...
.. ...
.
. ...
r
.
... ..
...
.
.
..
.. 1 ...
...
... ...
... ...
...
..
. ..
.
... ..
..
... ..
.
...
...
.
..
..
..
..
-
−2 .
...
.
.
−1 1
..
..
.. 2 x
.. ..
..
... ...
.
.
. ...
.... ...
...
.
.
...
... −1 ...
...
...
.... ...
... ..
...
... ..
... ..
... ..
..
r ...
...
... −2
..
..
..
..
... ..
..
... ..
x(ν+1) = ϕ(x(ν) ), ν = 0, 1, 2, . . .
erzeugte Folge {x(ν) } mit einem beliebigen Startwert x(0) ∈ I gegen den Fixpunkt ξ
der Abbildung ϕ, d. h. es gilt
lim x(ν+1) = ξ .
ν→∞
Beweis. Nach Satz 2.11 hat die Gleichung x = ϕ(x) in I genau eine Lösung ξ, also gilt
ξ = ϕ(ξ). Mit x(ν+1) = ϕ(x(ν) ) und der Lipschitzbedingung (2.13) erhält man
(LBed)
|x(ν+1) − ξ| = |ϕ(x(ν) ) − ϕ(ξ)| ≤ L|x(ν) − ξ|
2
ϕ(x) = −1 für x = 0
x
2
x(ν+1) = ϕ(x(ν) ) = −1, ν = 0, 1, 2, . . .
x(ν)
a) Wegen ϕ (x) = − x22 < 0, x = 0 , ist ϕ streng monoton fallend (Abb. 2.9).
Als Intervall wird I = [−3, −1] gewählt. Wegen ϕ(−3) = − 53 ∈ I und
ϕ(−1) = −3 ∈ I gilt: ϕ erfüllt in I die Bedingung (i) des Satzes 2.11.
40 2. Lösung nichtlinearer Gleichungen
2 2 2
|ϕ (x)| = ≤ = =2 > 1.
x2 min(x2 ) 1
x∈I
in I ist die Voraussetzung (ii) des Satzes 2.11 erfüllt. Damit kann die Lösung
ξ2 = −2 iterativ mit der Schrittfunktion aus der Umformung (III) berechnet wer-
den.
y6
-
x
−4 −3 −2 −1
.... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ..
. .. −1
.. ..
.. ..
...........................
.................... .......................
.
. .. ..
............... .. ..
.. .................
−1 r
2 ........ .. ..
y= x
..
.
........
....... .
...... .
.
.
..
.. −2
.. ...... ..
.. ...... .
.. .. .... ...
.. .. ..... ..
.. .. .... ..
.. ..
.. ... ... ... ... ... ... ... ... ... ... ........ ... ...........
...
...
−3
...
...
...
...
...
y =x −4
...
...
...
..
Abb. 2.9. Intervall [−3, −1] wird verkleinert auf [−3, −1.5]
Zur Bestimmung der Lösung ξ1 = 1 ist die Umformung (III) ungeeignet. Wählt
man nämlich z. B. I = [0.5, 1.2], so gilt dort |ϕ (x)| > 1, so dass die Bedingung (ii)
des Satzes 2.11 verletzt ist.
1−L |x
= L 1−L − x(ν−1) | .
(ν)
L
|x(ν+m) − x(ν) | ≤ |x(ν) − x(ν−1) |. (2.18)
1−L
Setzt man in der letzten Zeile von (2.17) ν = 1, m = ν − 1 und schreibt anschließend
wieder ν statt ν, so erhält man |x(ν) − x(ν−1) | ≤ Lν−1 |x(1) − x(0) |. Geht man damit in
(2.18) ein, so ergibt sich
L Lν
|x(ν+m) − x(ν) | ≤ |x(ν) − x(ν−1) | ≤ |x(1) − x(0) | . (2.19)
1−L 1−L
1. die a posteriori-Fehlerabschätzung
L
|∆(ν) | = |x(ν) − ξ| ≤ |x(ν) − x(ν−1) | = ε1 (2.20)
1−L
und
2. die a priori-Fehlerabschätzung
Lν
|∆(ν) | = |x(ν) − ξ| ≤ |x(1) − x(0) | = ε2 (2.21)
1−L
mit ε1 ≤ ε2 .
42 2. Lösung nichtlinearer Gleichungen
Die a priori-Fehlerabschätzung (2.21) kann bereits nach dem ersten Iterationsschritt vor-
genommen werden. Sie dient vor allem dazu, bei vorgegebener Fehlerschranke ε die An-
zahl ν der höchstens erforderlichen Iterationsschritte abzuschätzen, denn aus der Forde-
rung
Lν !
|x(ν) − ξ| ≤ |x(1) − x(0) | ≤ ε
1−L
ε(1−L)
ergibt sich mit Lν ≤ |x(1) −x(0) |
=: K und mit log L < 0 für 0 < L < 1 die Ungleichung
log K
ν≥ . (2.22)
log L
Die a posteriori-Fehlerabschätzung (2.20) kann erst im Verlauf oder nach Abschluss
der Rechnung genutzt werden, da sie x(ν) als bekannt voraussetzt; sie liefert eine bes-
sere Schranke als die a priori-Fehlerabschätzung und wird deshalb vorzugsweise zur
Abschätzung des Fehlers verwendet. Um rasche Konvergenz zu erreichen, sollten die
Schrittfunktion ϕ und das zugehörige Intervall I so gewählt werden, dass L < 0.2 gilt.
Wenn 1−LL
≤ 1 ist, also L ≤ 1/2, folgt aus (2.20) |x(ν) − ξ| ≤ |x(ν) − x(ν−1) |, d. h. der
absolute Fehler von x(ν) ist kleiner ( für L < 1/2 ) oder höchstens gleich (für L = 1/2)
der absoluten Differenz der letzten beiden Näherungen. Jeder Iterationsschritt bringt
demnach eine bessere Annäherung der gesuchten Lösung. Für 1/2 < L < 1 kann jedoch
der absolute Fehler von x(ν) größer sein als |x(ν) − x(ν−1) |, so dass hier die Iterationsfolge
noch nichts über den Fehler aussagen kann.
Im Falle monotoner Konvergenz (0 ≤ ϕ (x) < 1) folgt aus der Abschätzung |x(ν) − ξ| ≤ ε
(bzw. x(ν) − ε ≤ ξ ≤ x(ν) + ε) eine schärfere Eingrenzung der Lösung: x(ν) ≤ ξ ≤ x(ν) + ε,
falls die x(ν) von links gegen ξ streben, x(ν) − ε ≤ ξ ≤ x(ν) , falls die x(ν) von rechts gegen
ξ streben.
x(ν) ξ x(ν+1)
(ν+1) 1 (ν+1)
x − ξ ≤ x − x(ν) = α; (2.23)
2
damit hat man eine Möglichkeit, den Fehler ohne Kenntnis der Lipschitzkonstante
abzuschätzen.
2.3 Allgemeines Iterationsverfahren 43
Sei x(ν) eine iterativ bestimmte Näherung für die Nullstelle ξ ungerader Ordnung
von f (Abb. 2.10) und gelte für ein vorgegebenes ε > 0
so folgt daraus nach dem Zwischenwertsatz von Bolzano (Satz 2.4), dass im Intervall
(x(ν) − ε, x(ν) + ε) eine Nullstelle ξ liegen muss. Damit gilt die Fehlerabschätzung
|x(ν) − ξ| < ε.
y
6 y = f (x)
.....................
............
........
.......
......
......
r
.....
......
.
....
.
.... ..
..... ..
x (ν)
−ε x (ν) ξ ....
.
.... .
... x (ν)
+ε
..
....
.... .. -
.
..
.. ..
.. ....
....
x
.. ....
.. ε ...
...... ε
.. ...
.....
.. .....
.. .
...
......
..
.. ......
r . ........
.....................................
Abb. 2.10. f (x(ν) −ε) · f (x(ν) +ε) < 0 ⇒ |x(ν) −ξ| < ε
Praktisch geht man bei der Fehlerabschätzung mit Bolzano wie folgt vor: Man setzt
zunächst ein ε fest, welches sich über das Abbruchkriterium für die Iteration sinn-
voll festlegen lässt (z. B. ε = 10−k , k ∈ IN). Für dieses ε prüft man die Bedingung
(2.24), wobei eine Rechnung mit doppelter Genauigkeit zu empfehlen ist. Ist (2.24)
erfüllt, so ist ε eine obere Schranke für den absoluten Fehler. Um eine möglichst
kleine obere Schranke zu erhalten, führt man die Rechnung noch einmal mit einem
kleineren ε durch (z. B. mit ε1 = 10−k−1 ). Ist (2.24) auch für ε1 erfüllt, so ist ε1 für
|x(ν) − ξ| eine kleinere obere Schranke als ε. Analog fährt man so lange fort, bis sich
(2.24) für ein εj nicht mehr erfüllen lässt (z. B. εj = 10−k−j ). Dann ist εj−1 (z. B.
εj−1 = 10−k−j+1 ) die genaueste Fehlerschranke, die man auf diese Weise erhalten
hat.
Es wurde bisher vorausgesetzt, dass für diese Art der Fehlerabschätzung eine Null-
stelle ξ ungerader Ordnung vorliegen muss. Unter Verwendung des Satzes 2.3 kann
man sie aber auch für Nullstellen ξ gerader Ordnung einsetzen, indem man anstelle
der Funktion f die Funktion g = f /f für die Fehlerabschätzung verwendet, da
ξ dann eine einfache Nullstelle von g ist. Statt (2.24) ergibt sich hier die analoge
Bedingung
g(x(ν) − ε) · g(x(ν) + ε) < 0.
44 2. Lösung nichtlinearer Gleichungen
Rechnungsfehler
Es sei ε(ν) der lokale Rechnungsfehler des ν-ten Iterationsschrittes, der bei der Berech-
nung von x(ν) = ϕ(x(ν−1) ) entsteht. Gilt |ε(ν) | ≤ ε für ν = 0, 1, 2, . . ., so ergibt sich für
den akkumulierten Rechnungsfehler des ν-ten Iterationsschrittes
ε
|r(ν) | ≤ , 0 ≤ L < 1.
1−L
Die Fehlerschranke ε/(1 − L) ist also unabhängig von der Anzahl ν der Iterationsschritte;
der Algorithmus (2.11) ist somit stabil (vgl. Anmerkung nach Definition 36).
Da sich der Gesamtfehler aus dem Verfahrensfehler und dem Rechnungsfehler zusam-
mensetzt, sollten Rechnungsfehler und Verfahrensfehler von etwa gleicher Größenordnung
sein. Dann ergibt sich aus (2.21) bei bekanntem L die Beziehung
Lν ε
|x(1) − x(0) | ≈ ;
1−L 1−L
mit Lν ≈ |x(1) −x
ε
(0) | und 0 < L < 1 folgt für die Anzahl der höchstens auszuführenden
Iterationsschritte
ε
ν ≥ log (1) / log L .
|x − x(0) |
4. Fehlerabschätzung unter Verwendung des Satzes von Bolzano, also ohne Verwen-
dung einer Lipschitzkonstante. Gegeben ist die Gleichung f (x) = x2 + x − 2 = 0.
Untersucht wird der absolute Fehler von x(16) = −2.0000114.
1. Wahl: ε1 = 10−3 . Dann folgt aus f (x(16) − ε1 ) > 0 und f (x(16) + ε1 ) < 0, dass
für den absoluten Fehler gilt: |x(16) − ξ| < 10−3 .
2. Wahl: ε2 = 10−4 . Hier ergibt sich ebenfalls f (x(16) − ε2 ) > 0, f (x(16) + ε2 ) < 0,
also |x(16) − ξ| < 10−4 .
3. Wahl: ε3 = 10−5 . Man erhält f (x(16) − ε3 ) > 0 und f (x(16) + ε3 ) > 0, also ist
ε3 zu klein gewählt. Ein weiterer Versuch wird unternommen mit der
4. Wahl: ε4 = 0.5·10−4 . Man erhält f (x(16) − ε4 ) > 0 und f (x(16) + ε4 ) < 0; also
gilt für den absoluten Fehler |x(16) − ξ| < 0.5·10−4.
46 2. Lösung nichtlinearer Gleichungen
Tatsächlich ist |x(16) − ξ| = 0.0000114. Man sieht, dass diese Fehlerabschätzung nach
Bolzano am besten geeignet ist, um den absoluten Fehler möglichst genau anzugeben.
Zudem ist nur sie praktikabel, denn meist ist es bedeutend schwieriger als beim vorlie-
genden Beispiel, eine möglichst kleine Lipschitzkonstante anzugeben.
Algorithmus 2.17.
Gesucht ist eine Lösung ξ der Gleichung f (x) = 0.
1. Schritt: Festlegung eines Intervalls I, in welchem mindestens eine Nullstelle von
f liegt.
2. Schritt: Äquivalente Umformung von f (x) = 0 in eine Gleichung der Gestalt
x = ϕ(x).
3. Schritt: Prüfung, ob die Funktion ϕ für alle x ∈ I die Voraussetzungen des Satzes
2.11 erfüllt.
4. Schritt: Aufstellung der Iterationsvorschrift x(ν+1) := ϕ(x(ν) ), ν = 0, 1, 2, . . . und
Wahl eines beliebigen Startwertes x(0) ∈ I.
5. Schritt: Berechnung der Iterationsfolge {x(ν) }, ν = 1, 2, . . . Die Iteration ist so
lange fortzusetzen, bis mit einer Schranke δ1 > 0 für den relativen Fehler
Beispiel 2.18.
√
Gegeben: f (x) = cos x + 1 − x, x ≥ 0.
Gesucht: Eine Lösung ξ der Gleichung f (x) = 0 mit dem allgemeinen Iterations-
verfahren.
2.3 Allgemeines Iterationsverfahren 47
1. Schritt: Überblick über Lage und Anzahl der Lösungen von f (x) = 0. Mit
√
f1 (x) := cos x + 1 , f2 (x) := x
gilt
f (x) = f1 (x) − f2 (x) = 0 ⇐⇒ f1 (x) = f2 (x) .
Daher liefert die Abszisse eines Schnittpunktes der Graphen von f1 und f2
eine Näherung für eine Lösung ξ der Gleichung f (x) = 0 (Abbildung 2.11).
y
6 y = f (x) 2 ..........
...............
y = f (x) 1 ..............
..............
...............
.
2 ...............
........
...... ............
........ .........
. .....
...........
.......
...... ............ ......
.....
....................
. ..
..
......
..... .
..... .......... .....
..... .................. ...
.
r ..........
..
........ ..........
. ..
....
....
....
1 ......
..
.......
.....
. ..
.. ....... ....
....
.
..
ξπ 3π
x
1 2 2 π 2 2π
Abb. 2.11.
Aus Abbildung 2.11 liest man ab ξ ∈ [1, π2 ]. Mit f (1.3) = 0.127 und f (1.5) =
−0.154 ergibt sich I = [1.3, 1.5] als ein engeres Einschlussintervall für die
Lösung ξ.
erweist sich wegen |ϕ (x)| > 2 für x ∈ I als unbrauchbar. Daher wird die
Auflösung
√
x = arccos( x − 1) =: ϕ(x)
gewählt mit
−1
ϕ (x) = √ " √ , x > 0.
2 x 1 − ( x − 1)2
48 2. Lösung nichtlinearer Gleichungen
3. Schritt: Genügt ϕ für x ∈ I den Voraussetzungen (i) und (ii) des Satzes 2.11?
zu (i):
Wegen ϕ (x) < 0 für x > 0 ist ϕ streng monoton fallend. Mit ϕ(1.3) =
1.43 ∈ I und ϕ(1.5) = 1.34 ∈ I gilt wegen der Monotonie auch ϕ(x) ∈ I für
alle x ∈ I = [1.3, 1.5]. Also ist die Bedingung (i) erfüllt.
zu (ii):
Um eine Lipschitzkonstante L zu ermitteln, wird ϕ (x) in I grob abgeschätzt:
|ϕ (x)| √ √
1
= √
2 x 1−( x−1)2
1√
≤ √ √
2 min x·min 1−( x−1)2
x∈I x∈I
= √ √ 1√ = 0.450 = L < 1
2 < 1.
2 1.3 1−( 1.5−1)2
5. Schritt: Wegen ϕ (x) < 0 in I liegt alternierende Konvergenz vor, so dass mit der
zugehörigen Fehlerabschätzung eine Abbruchbedingung konstruiert werden
kann, die eine bestimmte Anzahl genauer Dezimalen garantiert. Für drei ge-
naue Dezimalen erhält man mit (2.23)
1 (ν) !
|x(ν) − ξ| ≤ |x − x(ν−1) | ≤ 0.5·10−3
2
und damit die Abbruchbedingung
|x(ν) − x(ν−1) | ≤ 1·10−3 .
Iteration:
ν x(ν) |x(ν) − x(ν−1) |
0 1.3 .
1 1.430 157 740 .
2 1.373 629 308 .
3 1.397 917 137 .
4 1.387 435 119
5 1.391 950 063 0.004 514 944
6 1.390 003 705 0.001 946 358
7 1.390 842 462 0.000 838 757
6. Schritt: Fehlerabschätzungen
Annahme: ε1 = 0.5·10−3
f (x(7) − ε1 ) = f (1.390 342 462) = 0.000 348 . . . > 0
f (x(7)
+ ε1 ) = f (1.391 342 462) = −0.001 059 . . . < 0
=⇒ |x(7) − ξ| < 0.5·10−3 d. h. 3 genaue Dezimalen
Annahme: ε2 = 0.5·10−4
f (x(7) − ε2 ) = f (1.390 792 462) = −0.000 285 . . . < 0
f (x(7) + ε2 ) = f (1.390 892 462) = −0.000 426 . . . < 0
=⇒ ε2 zu klein!
|x(ν+1) − ξ|
lim = M. (2.27)
ν→∞ |x(ν) − ξ|p
Das Iterationsverfahren x(ν+1) = ϕ(x(ν) ) heißt dann ein Verfahren von mindestens
p-ter Ordnung; es besitzt genau die Ordnung p, wenn M = 0 ist.
Durch (2.27) wird also ausgedrückt, dass der Fehler der (ν+1)-ten Näherung ungefähr
gleich M -mal der p-ten Potenz des Fehlers der ν-ten Näherung ist. Die Konvergenzge-
schwindigkeit wächst mit der Konvergenzordnung. Bei p = 1 spricht man von linearer
Konvergenz, bei p = 2 von quadratischer Konvergenz und allgemein bei p > 1 von super-
linearer Konvergenz. Es gilt der
50 2. Lösung nichtlinearer Gleichungen
Satz 2.20.
Die Schrittfunktion ϕ sei für x ∈ I p-mal stetig differenzierbar. Gilt dann mit
lim x(ν) = ξ
ν→∞
Wegen ϕ(ξ) = ξ, ϕ (ξ) = . . . = ϕ(p−1) (ξ) = 0, ϕ(p) (ξ) = 0 erhält (2.28) die Form
(x(ν) − ξ)p (p)
x(ν+1) − ξ = ϕ (ξ) + O (x(ν) − ξ)p+1
p!
bzw.
x(ν+1) − ξ 1
= ϕ(p) (ξ) + O (x(ν) − ξ) .
(x − ξ)
(ν) p p!
Im Fall p = 1 ist das Iterationsverfahren x(ν+1) = ϕ(x(ν) ) mit ϕ (ξ) = 0 ein Verfahren
erster Ordnung, und wegen der Lipschitzbedingung (2.14) ist M = |ϕ (ξ)| ≤ L < 1.
Es gilt außerdem der folgende in [COLL1968], S.231 bewiesene Satz, der zur Konstruk-
tion von Iterationsverfahren beliebig hoher Konvergenzordnung eingesetzt werden kann.
2.5 Newtonsche Verfahren 51
Satz 2.21.
Sind x(ν+1) = ϕ1 (x(ν) ) und x(ν+1) = ϕ2 (x(ν) ) zwei Iterationsverfahren der Konver-
genzordnungen p1 bzw. p2 , so ist
x(ν+1) = ϕ1 ϕ2 (x(ν) )
ein Iterationsverfahren, das mindestens die Konvergenzordnung p1 · p2 besitzt.
ϕ1 (x) = ϕ(x),
ϕs (x) = ϕ ϕs−1 (x) für s = 2, 3, . . .
Die Funktion f sei im Intervall [a, b] stetig differenzierbar und besitze in (a, b) eine ein-
fache Nullstelle ξ; also sind f (ξ) = 0 und f (ξ) = 0. Ferner sei f (x) = 0 für alle x ∈ [a, b].
x(1) wird als eine gegenüber x(0) verbesserte Näherung für die Nullstelle ξ angesehen.
Mit x(1) wird das Verfahren in derselben Weise fortgesetzt. Damit ergibt sich die nächste
Näherung
f (x(1) )
x(2) = x(1) − .
f (x(1) )
Für dieses geometrisch plausible Verfahren werden im folgenden Satz hinreichende Bedin-
gungen für die Konvergenz der Folge x(0) , x(1) , x(2) , . . . gegen die Nullstelle ξ angegeben.
Satz 2.22.
Die Funktion f sei im Intervall [a, b] zweimal stetig differenzierbar und besitze in (a, b)
eine einfache Nullstelle ξ; also sind f (ξ) = 0 und f (ξ) = 0.
Dann gibt es ein Intervall I = (ξ − r, ξ + r) ⊂ [a, b], r > 0, so dass die Iterationsfolge
f (x(ν) )
x(ν+1) = x(ν) − , ν = 0, 1, 2, . . . , (2.29)
f (x(ν) )
für das Verfahren von Newton mit jedem Startwert x(0) ∈ I gegen die Nullstelle ξ
konvergiert, und zwar von mindestens zweiter Ordnung.
Beweis. Vergleicht man (2.29) mit (2.4) und (2.2), so ergibt sich für die Schrittfunktion
ϕ des Newton-Verfahrens
f (x)
ϕ(x) = x − . (2.30)
f (x)
(2.3) zeigt, dass hier g(x) = 1/f (x) gesetzt wurde.
Um den Fixpunktsatz 2.14 auf die Schrittfunktion (2.30) anwenden zu können, muss ein
Intervall I angegeben werden, so dass ϕ für alle x ∈ I den Voraussetzungen (i) und (ii)
des Satzes 2.11 genügt.
Da f stetig ist, gibt es wegen f (ξ) = 0 eine Umgebung der Nullstelle ξ, in der f (x) = 0
ist. Diese Umgebung
sei Is = (ξ − s, ξ + s) ⊂ [a, b], s > 0. Somit ist die Schrittfunktion
ϕ vgl. (2.30) stetig im Intervall Is .
f (x) f (x)
ϕ (x) = . (2.31)
f 2 (x)
Auch ϕ ist stetig in Is .
Aus
| ϕ (x) | ≤ L für alle x ∈ I (2.34)
mit 0 < L < 1 folgt nach (2.14), dass die Schrittfunktion ϕ lipschitzbeschränkt ist, und
damit gilt (ii).
Um (i) nachzuweisen, muss gezeigt werden, dass mit x ∈ I, also | x−ξ | < r, auch ϕ(x) ∈ I
ist, d. h. | ϕ(x) −ξ | < r.
Wegen ϕ(ξ) = ξ vgl. (2.30) ist mit (ii)
Nach dem Fixpunktsatz 2.14 konvergiert also die Iterationsfolge (2.29) für jeden Start-
wert x(0) ∈ I gegen die Nullstelle ξ.
Wegen (2.32) konvergiert das Newtonsche Verfahren (2.29) nach Satz 2.20 von minde-
stens zweiter Ordnung, d. h. mindestens quadratisch.
Aus dem Beweis des Satzes 2.22 geht hervor, dass das Newtonsche Verfahren nur für
solche Startwerte x(0) konvergiert, die genügend nahe bei der Nullstelle ξ liegen.
Wegen der speziellen Schrittfunktion (2.30) kann eine gegenüber der a posteriori-Fehler-
abschätzung (2.20) in Abschnitt 2.3.4 verbesserte Fehlerabschätzung für das Newtonsche
Verfahren angegeben werden.
Satz 2.23.
Die Newtonsche Iterationsfolge (2.29) besitze gemäß Satz 2.22 das Konvergenzinter-
vall I. Dann gilt unter Verwendung der a posteriori-Fehlerabschätzung (2.20) mit
max |f (x)|
1 x∈I
≤ M1
2 min |f (x)|
x∈I
f (x(ν) ) 1 f (
x(ν) )
− − ξ + x(ν)
= (ξ − x(ν) )2
f (x(ν) ) 2 f (x(ν) )
Wird (2.20)
(ν) L (ν)
x −ξ ≤ x − x(ν−1) (2.40)
1−L
auf der rechten Seite der Ungleichung (2.39) eingesetzt, so ergibt sich die Fehlerab-
schätzung (2.35).
(ν+m) 1 2m
x −ξ ≤ M1 | x(ν) − ξ | (2.41)
M1
gilt nach (2.39) für m = 1. Mit der Annahme, (2.41) gelte für m, wird (2.41) für m + 1
bewiesen (vollständige Induktion). Aus (2.41) für m + 1 folgt mit (2.39)
(ν+1+m) 1 2m
x −ξ ≤ M1 | x(ν+1) − ξ |
M1
1 1 2 2m
≤ M1 M1 | x(ν) − ξ |
M1 M1
1 2m+1
= M1 | x(ν) − ξ | .
M1
Wenn auf der rechten Seite der Ungleichung (2.41) | x(ν) − ξ | mittels (2.40) ersetzt wird,
entsteht die Fehlerabschätzung (2.36).
2.5 Newtonsche Verfahren 55
Mit (2.35) kann mit zwei Näherungen x(ν−1) und x(ν) der absolute Fehler | x(ν+1) − ξ |
der nächsten (noch nicht berechneten) Näherung x(ν+1) im Voraus abgeschätzt werden.
Eine Lipschitzkonstante L erhält man durch Abschätzung von (2.31) auf dem Inter-
vall I:
2
|ϕ (x)| = | f (x) f (x)/f (x) | ≤ L < 1 .
Beispiel 2.24.
1 (ν) a
x(ν+1) = x + (ν) , ν = 0, 1, 2, . . . (2.42)
2 x
Ferner ist
1 a
ϕ (x) = 1− 2 .
2 x
√
Als numerisches Beispiel wird a = 5 gewählt. Wegen 2 < 5 < 3 sei I = [2, 3]. Dann ist
für alle x ∈ I
1
ϕ (x) = 1 − 5 ≤ 1 1 − 5 ≤ 0.23 = L < 1 .
2 x2 2 9 2
Die Folge (2.42) mit a = 5 konvergiert also mit jedem Startwert x(0) ∈ I = [2, 3].
Mit der Lipschitzkonstante L = 0.23 ist L/(1 − L) = 0.2987, und die a posteriori-
Fehlerabschätzung (2.40) lautet
2
benötigt. Damit ist M1 = 2·4 = 0.25. Die Fehlerabschätzung (2.35) lautet
√
ν x(ν) | x(ν) − x(ν−1) | | x(ν) − 5 | (2.43) (2.44)
0 3.00000000 0.76393202
1 2.33333333 0.66666667 0.09726536 0.19913333
2 2.23809524 0.09523810 0.00202726 0.02844762 0.00991352
3 2.23606890 0.00202634 0.00000092 0.00060527 0.00020232
4 2.23606798 0.00000092 0.00000000 0.00000027 0.00000009
5 2.23606798 0.00000000 0.00000000 0.00000000 1.88 · 10−14
Die mit den Fehlerabschätzungen (2.43) und (2.44) ermittelten√Schranken für den absolu-
ten Fehler zeigen im Vergleich mit dem exakten Fehler | x(ν) − 5 |, dass die Abschätzung
(2.44) kleinere Schranken liefert.
Zur Berechnung von Nullstellen beliebiger transzendenter Funktionen eignet sich besser
ein Einschlussverfahren. Siehe dazu Abschnitt 2.9 Entscheidungshilfen.
2.5 Newtonsche Verfahren 57
Analog zum gedämpften Newton-Verfahren für nichtlineare Systeme (Abschnitt 6.3.1, Al-
gorithmus 9) lässt sich das gedämpfte Newton-Verfahren für Einzelgleichungen angeben.
Man führt für ν = 0, 1, 2, . . . folgende Schritte durch:
(i) Berechne ∆x(ν+1) := x(ν+1) − x(ν) = −f (x(ν) )/f (x(ν) ) Newtonschritt (2.29) .
(ii) Berechne für i = 0, 1, . . .
(ν+1) 1
xi := x(ν) + ∆x(ν+1) .
2i
(ν+1) (ν+1)
Wenn |f (xi )| < |f (x(ν) )| gilt, wird x(ν+1) := xi gesetzt.
Andernfalls wird der Schritt mit dem nächsten i wiederholt.
Wenn mit einem vorgegebenen imax ∈ IN für alle
i ≤ imax die obige Bedingung nicht erfüllt ist, wird (mit i = 0)
Die Funktion f sei im Intervall [a, b] genügend oft stetig differenzierbar und besitze in
(a, b) eine Nullstelle ξ der Vielfachheit j, j ≥ 2. Nach Satz 2.2 sind also
Für die Ableitung (2.31) der Schrittfunktion (2.30) des Newtonschen Verfahrens für ein-
fache Nullstellen ergibt sich mit (2.2), (2.4) und (2.6)
Das Newtonsche Verfahren, das für einfache Nullstellen mindestens quadratisch konver-
giert (Satz 2.22), konvergiert nach Satz 2.20 in der Umgebung einer mehrfachen Nullstelle
nur linear.
kann die quadratische Konvergenz auch in der Umgebung einer mehrfachen Nullstelle ξ
(j ≥ 2) beibehalten werden. Mit (2.2) und (2.4) lautet (2.47)
h0 (x)
ψ(x) = x − j (x − ξ) ,
h1 (x)
und es ist
ψ(ξ) = ξ .
Für die Ableitung der Schrittfunktion (2.47) ergibt sich mit (2.45)
f (x) f (x) h0 (x) h2 (x)
ψ (x) = 1 − j 1 − = 1 − j 1 − .
f 2 (x) h21 (x)
Mit (2.46) ist
1 1
ψ (ξ) = 1 − j 1− 1− = 1−j = 0,
j j
und daraus folgt nach Satz 2.20 die mindestens quadratische Konvergenz mit der Schritt-
funktion (2.47).
Wegen ψ (ξ) = 0 und wegen der Stetigkeit von ψ in einer Umgebung von ξ gilt: Es gibt
ein Intervall I = (ξ − r, ξ + r) ⊂ [a, b], r > 0, so dass für alle x ∈ I |ψ (x)| ≤ L < 1 ist.
Somit ist ψ in I lipschitzbeschränkt. Ferner ist wegen ψ(ξ) = ξ für alle x ∈ I
| ψ(x) − ξ | = | ψ(x) − ψ(ξ) | ≤ L | x − ξ | ≤ | x − ξ | < r ,
also ψ(x) ∈ I. Somit genügt die Schrittfunktion (2.47) den Voraussetzungen (i) und (ii)
des Fixpunktsatzes 2.14. Damit folgt der
f (x(ν) )
x(ν+1) = x(ν) − j , ν = 0, 1, 2, . . . ,
f (x(ν) )
mit jedem Startwert x(0) ∈ I von mindestens zweiter Ordnung gegen die j-fache
Nullstelle ξ konvergiert.
Die Anwendung des Newtonschen Verfahrens für mehrfache Nullstellen setzt die Kennt-
nis der Vielfachheit j der gesuchten Nullstelle ξ voraus; j wird allerdings nur in speziellen
Fällen bekannt sein.
Nach Satz 2.3 ist eine j-fache Nullstelle ξ (j ≥ 2) einer Funktion f eine einfache Nullstelle
der Funktion g mit g(x) = f (x)/f (x). Insofern kann man die Nullstelle von g mit dem
Newtonschen Verfahren für einfache Nullstellen ermitteln. Dieses Verfahren beschreibt
der Satz
2.5 Newtonsche Verfahren 59
lautet
g(x) f (x)
ϕ(x) = x − = x − J(x) .
g (x) f (x)
Mit den Taylorentwicklungen (2.2), (2.4) und (2.6) an der Stelle ξ sind
h0 (x)
ϕ(x) = x − J(x) (x − ξ) und (2.48)
h1 (x)
1
J(x) = h0 (x) h2 (x)
.
1− h21 (x)
(2) Zur Bestimmung der Konvergenzordnung wird ϕ (ξ) berechnet. Dabei werden (2.49)
und (2.8) verwendet. Mit (2.48) und ϕ(ξ) = ξ ist
ϕ(x) − ϕ(ξ)
ϕ (ξ) = lim
x→ξ x−ξ
1 h0 (x)
= lim x − J(x) (x − ξ) −ξ
x→ξ x − ξ h1 (x)
h0 (x) 1
= lim 1 − J(x) = 1−j = 0.
x→ξ h1 (x) j
(3) Ähnlich wie beim Satz 2.25 für (2.47) kann gezeigt werden, dass die Schrittfunktion
ϕ in einem Intervall I = (ξ−r, ξ+r), r > 0, den Voraussetzungen des Fixpunktsatzes
2.14 genügt.
Beispiel 2.27.
Gegeben: Die Funktion f : f (x) = 1 − sin x. f besitzt bei ξ = π
2 eine doppelte Nullstelle.
Gesucht: Die Nullstelle ξ mit dem Newtonschen Verfahren für einfache Nullstellen, mit
dem Newtonschen Verfahren für mehrfache Nullstellen und mit dem modi-
fizierten Newtonschen Verfahren. Es wird mit 15-stelliger Mantisse gerech-
net und die Rechnung abgebrochen, wenn die Bedingung |x(ν) − x(ν−1) | ≤
0.5·10−14 erfüllt ist. Als Startwert wird x(0) = 2 gewählt.
2. Newtonsches Verfahren für doppelte Nullstellen (Satz 2.25 mit Vielfachheit j = 2):
f (x(ν) ) 1 − sin(x(ν) )
x(ν+1) = x(ν) − 2 = x(ν)
− 2
f (x(ν) ) − cos(x(ν) )
Hier ist die Abfrage |x(ν) − x(ν−1) | ≤ 0.5·10−14 bereits nach vier Iterationsschritten
erfüllt, da das Verfahren quadratisch konvergiert. Aber auch hier gilt für den abso-
luten Fehler des Näherungswertes x(4) |x(4) − π2 | ≤ 0.13·10−11; die Ursache dafür
ist dieselbe wie im ersten Falle.
62 2. Lösung nichtlinearer Gleichungen
Für den absoluten Fehler des Näherungswertes x(4) gilt |x(4) − π2 | ≤ 0.25·10−12,
obwohl für die iterierten Werte die Abfrage |x(4) − x(3) | ≤ 0.5·10−14 erfüllt ist. Die
Ursache liegt wie oben darin, dass sin(x(3) ) in den mitgeführten Stellen exakt 1 ist,
d. h. 1 − sin(x(3) ) verschwindet, aber cos(x(3) ) = 0 ist, so dass x(4) = x(3) gilt.
Man sieht hier außerdem, dass J(x(2) ) der tatsächlichen Vielfachheit j = 2 schon
sehr nahe ist, J(x(3) ) jedoch gleich 1 wird. Dies liegt daran, dass
f (x(3) ) = 1 − sin(x(3) ) in den mitgeführten Stellen verschwindet, f (x(3) ) = 0
ist und damit f f /f 2 verschwindet. Würde hier f 2 vor f Null, so würde J(x(ν) )
von j in beliebiger Weise abwandern.
Wegen f (x) = 1 − sin x kann in diesem Fall J(x(ν) ) vereinfacht dargestellt werden
(Satz 2.26):
1
J(x(ν) ) = = 1 + sin(x(ν) ).
1−sin(x(ν) ) sin(x(ν) )
1− cos2 (x(ν) )
Wegen seiner geringen Effizienz (vgl. Abschnitt 2.9) sollte man mit dem modifizierten
Newton-Verfahren nur so lange iterieren, bis die Vielfachheit der Nullstelle geklärt ist,
dann aber mit dem Verfahren von Newton für mehrfache Nullstellen weiterrechnen. Die
Vielfachheit ist dann geklärt, wenn sich von einem gewissen ν an entweder J(x(ν) ) = 1
ergibt oder
|J(x(ν) ) − J(x(ν−1) )| > |J(x(ν−1) ) − J(x(ν−2) )|
gilt, d. h. J(x(ν) ) sich wieder von der Vielfachheit j entfernt. In beiden Fällen ist die zu
J(x(ν−1) ) nächste ganze Zahl die gesuchte Vielfachheit j. Dieses Verhalten der J(x(ν) ) ist
bedingt durch die beschränkte Stellenzahl der Maschinenzahlen und den für x(ν) → ξ un-
2
bestimmten Ausdruck f f /f im Nenner von J(x(ν) ). Wird nämlich für ein x(ν) wegen
der beschränkten Stellenzahl f (x(ν) ) identisch Null, während f (x(ν) ) noch verschieden
von Null ist, so erhält man J(x(ν) ) = 1, obwohl eine mehrfache Nullstelle vorliegt.
2.6 Das Sekantenverfahren 63
Generelle Empfehlung. Da eine doppelte Nullstelle ξ von f nach Satz 2.3 eine einfache
Nullstelle ξ von g(x) = f (x)/f (x) ist, kann auch zur Berechnung einer doppelten Null-
stelle ein Einschlussverfahren verwendet werden. Für Beispiel 2.27 würde das bedeuten,
dass man anstelle der doppelten Nullstelle von f (x) = 1 − sin x die einfache Nullstelle
der Funktion
f (x) 1 − sin x
g(x) = =
f (x) − cos x
mit einem Einschlussverfahren ermittelt, das mit größerer Effizienz arbeitet als das
Newton-Verfahren (siehe Abschnitt 2.9). Die oben erwähnten numerischen Probleme we-
gen g(ξ) = 00 treten natürlich auch hier auf.
Beispielsweise hat die quadratische Funktion f (x) = x2 +px+q genau dann eine doppelte
Nullstelle, wenn exakt p2 − 4q = 0 ist. Wenn p und q mit Messfehlern behaftet sind, wird
| p2 − 4q | = ε > 0 sein. Für p2 − 4q > 0 gibt es dann zwei benachbarte einfache Nullstellen
und für p2 − 4q < 0 keine.
Eine doppelte Nullstelle ξ von f mit f (ξ) = f (ξ) = 0, f (ξ) = 0 ist eine Stelle, an der die
Funktion f ein lokales Minimum oder Maximum besitzt. Wegen f (ξ) = 0 ist ξ einfache
Nullstelle von f und kann mit einem Einschlussverfahren berechnet werden. Wenn dann
|f (ξ)| ≤ ε ist, kann ξ als doppelte Nullstelle von f akzeptiert werden.
Das Sekantenverfahren wird hier behandelt, weil es in fast jeder Numerik-Vorlesung als
Standardverfahren vorkommt. Für die praktische Anwendung ist es allerdings nicht zu
empfehlen, man sollte dem Sekantenverfahren die Einschlussverfahren mit höherer Effi-
zienz unbedingt vorziehen.
Die Funktion f sei in I = [a, b] stetig und besitze in (a, b) eine einfache Nullstelle ξ.
Zur näherungsweisen Bestimmung von ξ mit Hilfe des Verfahrens von Newton ist die
Berechnung der Ableitung f von f erforderlich, so dass die Differenzierbarkeit von f
vorausgesetzt werden muss. Das Sekantenverfahren ist ein Iterationsverfahren, das ohne
Ableitungen arbeitet und zwei Startwerte x(0) , x(1) erfordert (siehe auch Abschnitt 2.7
Einschlussverfahren).
64 2. Lösung nichtlinearer Gleichungen
x(1) − x(0)
x(2) = x(1) − f (x(1) ).
f (x(1) ) − f (x(0) )
y6 y = f (x)............................
.......... (1)
.
... ......
...r
.......
x , f (x(1) )
..
............
.
... . . .... .. ..
..... ...
..... .... ...
...... ...
.
........ ...... ...
.
...... .... ..
..... ......
..... ....
..
...... ....
........
. ..... ...
..... ..... ..
(2).............. ..............
x (0) x ..
... r r .
......
..
.
....
.....
. ..
......
.
....
...
...... .. -
... ..... .............
a ......
... .......... ................
ξ
.
x (1)
b x
r
.
.... ........
. .
...................
.................. (0) (0)
.. x , f (x )
Abb. 2.13. Sekantenverfahren
Falls f (x(ν+1) ) = 0 ist, wird das Verfahren mit ξ = x(ν+1) abgebrochen. Wesentlich für
die Konvergenz des Verfahrens ist, dass die Startwerte x(0) , x(1) hinreichend nahe an der
Nullstelle ξ liegen. Es gilt der folgende Konvergenzsatz :
Satz 2.28.
Falls die Funktion f für alle x ∈ (a, b) zweimal stetig differenzierbar ist und mit zwei
positiven Zahlen m, M den Bedingungen
|f (x)| ≥ m, |f (x)| ≤ M, x ∈ (a, b),
genügt, gibt es immer eine Umgebung Ir = [ξ − r, ξ + r] ⊂ (a, b) , r > 0 , so dass ξ in
Ir die einzige Nullstelle von f ist und das Verfahren für jedes Paar von Startwerten
x(0) , x(1) ∈ Ir , x(0) = x(1) , gegen die gesuchte Nullstelle ξ konvergiert.
Da die Überprüfung der Voraussetzungen dieses Satzes meist nicht praktikabel ist, bleibt
nur die Empfehlung, die Startwerte x(0) und x(1) möglichst nahe bei der Nullstelle ξ zu
wählen.
Wenn das Startintervall den Bedingungen des Satzes 2.28 nicht genügt, kann beim Sekan-
tenverfahren die Konvergenz von der Bezeichnung der Startwerte x(0) und x(1) abhängen.
2.6 Das Sekantenverfahren 65
...........................
...............
r
.......
.........
.......
.............................................................. .......
......
.
....... .
. ......
.r
................................
....
-
....... ...... ....................................
..
.....
.
.. .
......... ...
..... ....... ...
f (x) ..
...... .............. ...
. x (3)
.... .......
..
..... ..
..........
.... .
......
r r r -
.. ...
..
.... ......
.... .......
.... .......
.... ....... x
(0) ... ....... (2) (1)
x ..
....
.. ...
.......
..
...
. x x
... .......
....
. ...
..........
.
..
... .......
... ......
.......
... .......
.
....
. ...
.........
. .
... .......
... .......
... ......
.. r..........
......
...
..
...
...
...
Abb. 2.14. Startwerte, für die das Sekantenverfahren divergiert
.............................................
r
...........
........
......... ......
.....
.
............
....
.......... r ....
.
.......
....... ......
.
.. ......
...... .....
. ...... ...
... ..... .... ...
f (x) ...
..... .... ......
.......
..
..... ...
...... ......
. ......
..... ....
. ..........
.
.... ... .......
r .... ...
.... ....
.... ...
r .......
r ......
....... r -
(1) ..
..... ....... (3) ..
.......... (2) (0) x
x .. ..
.... ....
.... ...
x ....
.......
.......
x x
.... ......
. .
...
........
.. .. ......
... .... .......
... .... .......
... ... .......
.
.......... ............
.
... ... ...... .
...... .......
........ ........
...............
..
...
.
..r
...............
...
..
...
.
...
Abb. 2.15. Startwerte, für die das Sekantenverfahren konvergiert
Beispiel 2.29.
Gegeben: Die Funktion f (x) = 18 x2 − x + 32 und die Startwerte x(0) = 1, x(1) = 5 , die
wegen f (1) = 0.625, f (5) = −0.375 eine Nullstelle von f einschließen.
Gesucht: Die zwischen x(0) und x(1) liegende Nullstelle mit dem Sekantenverfahren.
Lösung: Mit den gegebenen Startwerten konvergiert das Sekantenverfahren gegen die
Nullstelle 6, die nicht im Einschlussintervall [1, 5] liegt. Mit den Startwerten
x(0) = 5, x(1) = 1 liefert das Sekantenverfahren die Nullstelle 2.
Ein Einschlussverfahren (siehe Abschnitt 2.7) ermittelt in beiden Fällen die Null-
stelle 2.
Beispiel 2.30.
√
Gegeben: Die Funktion f (x) = ln x − x + 1.5, die für x > 0 reell und stetig ist, und
die Startwerte x = 0.2, x = 2, die wegen f (0.2) ≈ −0.557, f (2) ≈ 0.779
(0) (1)
Die Konvergenzordnung
√ des Sekantenverfahrens mit der Iterationsvorschrift (2.50) ist
p = (1 + 5)/2 ≈ 1.62 . Prinzipiell kann die Vorschrift (2.50) auch zur näherungsweisen
Berechnung mehrfacher Nullstellen verwendet werden, dann geht jedoch die hohe Kon-
vergenzordnung verloren. Das modifizierte Sekantenverfahren (Abschnitt 2.6.2) besitzt
auch bei mehrfachen Nullstellen die Konvergenzordnung p ≈ 1.62. Zur Effizienz der Ver-
fahren siehe Abschnitt 2.9.
2.7 Einschlussverfahren
Eine stetige Funktion f : [a, b] → IR mit f (a) · f (b) < 0 besitzt nach dem Zwischenwert-
satz im Intervall (a, b) mindestens eine Nullstelle; ein solches Intervall [a, b] nennt man
ein Einschlussintervall. Im Folgenden sei das Intervall [a, b] so gewählt, dass es genau eine
Nullstelle ξ der Funktion f einschließt.
Bei den in diesem Abschnitt behandelten Einschlussverfahren wird ein gegebenes Ein-
schlussintervall in zwei Teilintervalle zerlegt; eines von diesen ist wieder ein Einschluss-
intervall. So können fortlaufend kleinere Einschlussintervalle erzeugt werden (Abschnitt
2.7.1).
2.7 Einschlussverfahren 67
Das einfachste dieser Verfahren ist das Bisektionsverfahren (Abschnitt 2.7.2), das stets,
aber nur linear, gegen die Nullstelle konvergiert.
Bei der Regula falsi (Abschnitt 2.7.3), die ebenfalls nur linear konvergiert, kann es vor-
kommen, dass die erzeugten Grenzen der Einschlussintervalle sich der Nullstelle nur von
einer Seite her nähern (Abbildung 2.17).
Deshalb wurden einige Einschlussverfahren (Abschnitte 2.7.4 bis 2.7.6) entwickelt, die
dieses Verhalten vermeiden und zudem eine höhere Konvergenzordnung besitzen. Aller-
dings wird diese höhere Konvergenzordnung erst in einer hinreichend kleinen Umgebung
der Nullstelle wirksam.
Daher ist es zweckmäßig, ein gegebenes Einschlussintervall zunächst mit dem Bisektions-
verfahren zu verkleinern und dann ein Verfahren höherer Ordnung einzusetzen (siehe
dazu das Beispiel 2.41).
Falls eine genügend oft stetig differenzierbare Funktion f im Intervall [a, b] eine Nullstelle
ξ gerader Ordnung besitzt, ist diese wegen Satz 2.3 eine einfache Nullstelle der Funktion
g(x) = f (x)/f (x) und kann daher mit einem Einschlussverfahren, angewandt auf g,
ermittelt werden. Siehe dazu auch die Anmerkungen zu mehrfachen Nullstellen am Ende
des Abschnitts 2.5.
erzeugt, die zwischen x(1) und x(2) liegt. Die Verfahren unterscheiden sich in der Wahl
von q. Mit x(3) wird der Funktionswert f3 = f (x(3) ) berechnet. Wenn f3 = 0 ist, ist
ξ = x(3) Nullstelle von f .
Wenn f3 = 0 ist, hat f entweder zwischen x(2) und x(3) oder zwischen x(1) und x(3) einen
Vorzeichenwechsel. Eines der beiden durch x(3) erzeugten Teilintervalle ist also wieder
ein Einschlussintervall, mit dem das Verfahren fortgesetzt werden kann. Das neue Ein-
schlussintervall wird wie folgt ermittelt.
Wenn f2 · f3 < 0 ist, liegt ξ zwischen x(2) und x(3) , und die Intervallgrenzen und Funk-
tionswerte werden umbenannt:
x(1) := x(2) , f1 := f2 ,
x(2) := x(3) , f2 := f3 .
Wenn f2 · f3 > 0 ist und somit f1 · f3 < 0 ist, liegt ξ zwischen x(1) und x(3) . Dann werden
gesetzt:
x(2) := x(3) , f2 := f3 .
68 2. Lösung nichtlinearer Gleichungen
Nun gilt wieder f1 · f2 < 0, und mit x(1) und x(2) kann das Verfahren fortgesetzt werden.
Dabei ist x(2) die zuletzt berechnete Intervallgrenze.
Da ξ zwischen x(1) und x(2) liegt, gelten für die absoluten Fehler dieser Intervallgrenzen
Wenn mit einer positiven Schranke RelErr für den relativen Fehler mit der zuletzt be-
rechneten Intervallgrenze x(2) (= 0)
Von den beiden Intervallgrenzen x(1) , x(2) kann man diejenige als die beste Näherung für
die Nullstelle ξ wählen, für die der Funktionswert dem Betrage nach kleiner ist.
Die Fehlerschranken AbsErr und RelErr müssen größer als die Maschinengenauigkeit
sein (etwa 2 bis 3 ).
Bemerkung. Für die praktische Anwendung ist es sinnvoll, die Abbruchbedingung (2.54)
für den relativen Fehler zu verwenden, weil so eine Aussage über die Anzahl der gültigen
Ziffern einer Näherungszahl gemacht werden kann (vgl. Satz 1.24). Der absolute Fehler
macht eine Aussage über die Anzahl der gültigen Dezimalen (vgl. Definition 17).
| x(2) − ξ |
≤ RelErr = 5 · 10−m
| x(2) |
liefert x(2) mit m genauen Ziffern, beginnend mit der ersten von 0 verschiedenen Ziffer
von x(2) .
2.7 Einschlussverfahren 69
Bei diesem einfachsten Einschlussverfahren wird in (2.52) q = 0.5 gesetzt. Dann wird
wegen
| x(3) − x(2) | = 0.5 | x(1) − x(2) |
die Länge des Einschlussintervalls halbiert. Mit fortgesetzter Intervallhalbierung konver-
giert das Bisektionsverfahren linear gegen die Nullstelle.
y6
.........
.....................
...............
.............
f (x) ......................................... ..
.
.
......... ...
........
........
..
...
........ ..
..
.
.......
....... ...
......
........
..
..
...... . ....
.. .
..... ...
..... ...
..... .
......
.... ...
.
.. ... ....
.... . .
..... .. .. ....
.... .. .
r
... ....... ..
... ....
...... .
..
..
.
..
. -
x (1)
=a x (4) .
.....
.
ξx x
...
(6) (5)
x (3)
x (2)
=b x
... ...
....
. .....
... ..
.......
.
. ......
.....
... ......
......
.
...
........
.................
.
.........
............
...........................
......................................................................................................................................................................................................................................................................................
...........................................................................................................................................
. . .
...........................................................................
......................................
. .
|b − a| 1 −k
≤ ·10 .
2n 2
Dabei sei 0.5 · 10−k < |b − a|. Dann folgen
2n ≥ 2 |b − a| 10k ,
n lg 2 ≥ lg 2 + lg |b − a| + k ,
k + lg |b − a|
n ≥ 1+ . (2.55)
lg 2
Beispiel 2.33.
Gegeben: Die Funktion f (x) = sin x + 1 − 1/x für x ∈ [a, b] mit a = 0.6, b = 0.7 und die
Schranke AbsErr = 0.5·10−6 für den absoluten Fehler.
Gesucht: Eine Nullstelle ξ der Funktion f im Intervall (0.6, 0.7) mit Hilfe des Bisektions-
verfahrens auf 6 Dezimalen genau.
Lösung: Wegen f (0.6) ≈ −0.102 und f (0.7) ≈ 0.216 liegt mindestens eine Nullstelle in
(0.6, 0.7). Mit k = 6, a = 0.6 und b = 0.7 erhält man nach der Formel (2.55)
2.7 Einschlussverfahren 71
In dieser Tabelle steht jeder Schritt in einer Zeile. Die neue Näherung ist x(3) . Nach Um-
speichern ergeben sich rechts die Grenzen x(1) und x(2) des neuen Einschlussintervalls.
In der letzten Spalte kann man die Abbruchbedingung prüfen.
Dieses Verfahren verwendet die Sekante, die die Punkte (x(1) , f1 ) und (x(2) , f2 ) verbin-
det, und bestimmt deren Schnittpunkt (x(3) , 0) mit der x-Achse (Sekantenschritt). Mit
der Gleichung
f1 − f2
y = f2 + (1) (x − x(2) )
x − x(2)
dieser Sekante ergibt sich für y = 0
f2
x(3) = x(2) + (x(1) − x(2) ) .
f2 − f1
Hier ist also mit (2.52)
f2
q= ,
f2 − f1
und wegen f1 · f2 < 0 ist 0 < q < 1.
72 2. Lösung nichtlinearer Gleichungen
Wie die Abbildung 2.17 zeigt, kann der Fall eintreten, dass die Grenzen der Einschluss-
intervalle sich der Nullstelle nur von einer Seite her nähern, während die Grenze x(1) auf
der anderen Seite der Nullstelle unverändert bleibt.
dem Betrage nach immer kleiner. Um von einer Intervallgrenze x(2) nahe der Nullstelle
auf die andere Seite der Nullstelle zu gelangen und damit die Grenze x(1) loszuwerden,
kann man wie folgt vorgehen:
Es sei
gesetzt. Wenn dann x(2) und x(3) = x(2) + ∆x die Nullstelle einschließen, ist wegen
| ∆x | < tol auch die Abbruchbedingung
erfüllt. Der Faktor 0.9 verhindert, dass die Abbruchbedingung infolge von Rundungsfeh-
lern evtl. nicht erfüllt ist.
Bemerkung. Dieser zusätzliche Schritt wird auch in die folgenden Algorithmen zum
Pegasus-Verfahren und zum Verfahren von Anderson-Björck aufgenommen.
.......... .....
r
.............................................
......
........ ......... .....
f (x)
...
.. ...... .
.... .....
...... . ....
.....
...
...... ......... ...
...
... ... ...
r
..
. .
..
r .... . ..... ...
.
.. ...
. ...
...... .
.....
. ......... .....
.
...
.
..... ...
..... ...
......... ..
..........
. .
.....
. ... .
. .......
(2).............. . .
..... ........
r x r .
......
.......
r .
.....
........
. r
........
................ r -
(1) ........ ..... ........ (2) x
x ....... ......... (2)
x ........
........
x x (2)
....... .. ........
.. .......... ......... ....
...
.
......... ....... .
...
........
.
..... ....
..... ..... ........
...... ..... ........
........
...... ...... ........
........
.................... .
. .....
...... .... ........
........... ........
............. ..........
......................
.
...
. r
...................
.
...
...
....
..
Berechne
tol := | x(2) | RelErr + AbsErr
und
f2
∆x := (x(1) − x(2) ) .
f2 − f1
Wenn | ∆x | ≤ tol ist, setze
∆x := 0.9 · tol · sgn (x(1) − x(2) ) .
Berechne x(3) := x(2) + ∆x.
2. Berechnung des neuen Funktionswertes f3 := f (x(3) ).
Falls f3 = 0 ist, wird die Iteration mit ξ := x(3) abgebrochen, andernfalls geht es
mit 3. weiter.
3. Festlegung des neuen Einschlussintervalls:
Falls f2 · f3 < 0 ist, liegt ξ zwischen x(2) und x(3) , und es wird gesetzt
x(1) := x(2) , x(2) := x(3) , f1 := f2 , f2 := f3 ;
falls f2 · f3 > 0 ist, liegt ξ zwischen x(1) und x(3) , und es wird gesetzt
x(2) := x(3) , f2 := f3 .
(1)
In beiden Fällen liegt jetzt ξ zwischen x und x(2) , und x(2) ist der zuletzt
berechnete Näherungswert.
4. Prüfung der Abbruchbedingung:
Falls
| x(1) − x(2) | ≤ tol
ist, erfolgt Abbruch. Dann wird gesetzt
ξ := x(2) , falls |f2 | ≤ |f1 | ist , und sonst ξ := x(1) .
Andernfalls, also mit | x(2) − x(1) | > tol, wird die Iteration mit 1. fortgesetzt.
74 2. Lösung nichtlinearer Gleichungen
Gesucht: Die Nullstelle ξ der Funktion f im Intervall (0.6, 0.7) mit der Regula falsi mit
7 gültigen Ziffern.
Lösung:
⇒ ξ ≈ x(2) = 0.6294464
7 genaue Ziffern
Das Pegasus-Verfahren ist ein gegenüber der Regula falsi verbessertes Verfahren, das in
einer hinreichend kleinen Umgebung einer einfachen Nullstelle die Konvergenzordnung
p = 1.642 besitzt. Falls nach einem Sekantenschritt die älteste Intervallgrenze x(1) nicht
ersetzt werden kann, wird der Funktionswert f1 mittels g · f1 , 0 < g < 1, modifiziert.
Damit verbessert sich die Chance, die Grenze x(1) nach dem nächsten Schritt ersetzen
zu können. Siehe dazu die geometrische Interpretation und Abbildung 2.18.
Basierend auf den Originalarbeiten [DOWE1971] und [DOWE1972] wurde hier der fol-
gende Algorithmus entwickelt.
2.7 Einschlussverfahren 75
Für den modifizierten Schritt wird nachfolgend eine geometrische Interpretation angege-
ben.
geometrisch wie folgt (Abbildung 2.18): Die Verbindungsgeraden der Punkte (x(2) , f2 +f3 )
und (x(1) , f1 ) sowie (x(2) , f2 ) und (x(1) , f1∗ ) schneiden sich in einem Punkt S der x-Achse.
Mit dem Punkt (x(1) , f1∗ ) verbessert sich die Chance, das Einschlussintervall auf der Seite
von x(1) zu verkürzen. Das zeigt die nächste Näherung x(4) in Abbildung 2.18.
...
..............
q (x(2) , f2 +f3 )
........... ....
...............
..
.. ... . . ....
..... ....... ...
..... ......... ...
..... .
....
.. ...
..... .
... ...
....
. . .
..... ...
..... ...
..... ...
..... ...
.... ...
..... .
...... ....
..... .......
. ..
. ..... ..
...... ...
.... ...
.. .....
... .................... ...... ...... ..
. ..
. ...
.
.
. .
.
...
.
.
.... ....
..........................
................
f (x) ...
....
...
............. .. ... ....
... .......... ..
....
. .
...
.......... . .
....
. . . . . .......... ... (2)
q (x , f2 )
...... ... ... ........ ...
(3) ...... .... .
(x ,f ) 3 ............. .
..
.. .
........ . . .............
q ... .. . ..
... ... . ..... .. .......
. . .
...... ..... .............. .... ...
...... .. ..... ...............
. .
.. ............. .... .
.......
. ........................... ..
...
. .
........... .... .....
. ........ ....... ...
.................
. .
. ... ...... ..........................
. ...
...
.......... .
.
. ...
.... ..
.......... .......... ...
. ..
.......... .
.
. ...
..... . .....
. ..... .......... ...
. . ........... .
.
. ..
...... ....
. ...... ..........
. . .
. . .......... .
.
. .
..
.... ............ .......... .
.
.
............ .. .........
. ..... ....... .. .
....... .
.
.
. . ......... ... .......
. ..... ....... .....
..... .
.
.
... .......
.
.. ..
........ ...... ........ ......
...
......... ...... ........ ....... ....
..... .. .......
......
.......... ..... .... ........ ....... ...
. ..
.. ..
.. ..... .
. .
...
..........
. .
...
........ ...
.. ..... .... ... ............
.
.. ... ..
..... ...
...... ..
..... . . . .. ......... .......... ...
.
....... . .
... ............. ... .....
.
(4) ............ ... . ..
... ............ . ....... . ...
x ....
.................
.
.
. ......... .
...
.. ...
...
...
q .
..
.
.qr
..
.. ...
.........
.. ....
...... .
............
.... q ...
........
.
. .. .
. .... ... q .
. .
........ ...
. q -
... .
... .... . .
.............. . .... .
.. . .. ...
x(1) ...
... .
.. .
.
.
.
... ....
. ξ ... . . . .
... ...
..... .......
.....
. S .......
........... ..
x (3) (2)
x x
..... ..... ................ ..... .......
... ..... . .... ..... .......
... ..... .................. ..... .......
... .... . ..... ......
... ............................ .... .....
..... ......
(x(1) , f ) q
. .. ... .
..... .......
∗ ..................... ....... ..... .......
1 ..... . ..... .......
..... ..
......
. .....
........
... .
... ..... ......
. ..... .......
..... ..... .......... .............
... .... ...... .........
... ... ....................
.......................
(x(1) , f ) q
....
.......
1 .....
..
.
...
..
...
Gegeben: Die Funktion f (x) = sin x + 1 − 1/x und die Schranke RelErr = 5 · 10−7 für
den relativen Fehler (AbsErr = 0).
Gesucht: Die Nullstelle ξ von f im Intervall (0.6, 0.7) mit Hilfe des Pegasus-Verfahrens.
2.7 Einschlussverfahren 77
Lösung:
Zur Darstellung dieser Tabelle vgl. Beispiel 2.33. Wegen |x(2) − x(1) |/|x(2) | < 5 · 10−7 für
ν = 4 ist die Nullstelle auf 7 Ziffern genau, also ξ ≈ x(1) = 0.629 446 5 .
Das Verfahren von Anderson-Björck arbeitet ähnlich wie das Pegasus-Verfahren; lediglich
im 3. Schritt des Algorithmus (bei der Festlegung des neuen Einschlussintervalls) wird die
Modifikation des Funktionswertes f1 an der Stelle x(1) auf andere Weise vorgenommen.
Die Konvergenzordnung des Verfahrens in der Umgebung einer einfachen Nullstelle liegt
zwischen 1.682 und 1.710.
oder, falls g ≤ 0 ist, mit g = 0.5 zugeordnet (modifizierter Schritt). Dann wird
gesetzt
x(2) := x(3) , f2 := f3 .
Jetzt liegt ξ zwischen x(1) und x(2) , und x(2) ist die zuletzt berechnete Näherung.
4. Prüfung der Abbruchbedingung:
Falls
| x(1) − x(2) | ≤ tol
ist, erfolgt Abbruch. Dann wird gesetzt
Die Originalarbeit zu dem Verfahren von Anderson-Björck ist [ANDE1973]. Auf dieser
Basis wurden hier zusätzlich der Algorithmus formuliert und die geometrische Interpre-
tation entwickelt.
Die Ersetzung von f1 durch gf1 im 3. Schritt des Algorithmus lässt sich wie folgt geo-
metrisch interpretieren (Abbildung 2.19): Durch die drei Punkte
wird die interpolierende quadratische Parabel gelegt, und im (mittleren) Punkt P3 wird
die Parabeltangente konstruiert. Falls diese Tangente die x-Achse zwischen x(3) und x(1)
schneidet, wird dieser Schnittpunkt x(4) als die nächste Näherung für die gesuchte Null-
stelle genommen.
Die Parabeltangente kann wie folgt konstruiert werden. H1 sei der Schnittpunkt der Ge-
rade P2 P3 mit der Gerade x = x(1) und H2 der Schnittpunkt der Gerade P1 P3 mit der
Gerade x = x(2) . Die Parabeltangente durch P3 ist parallel zur Verbindungsgerade H1 H2 .
Im Algorithmus wird die Tangente durch den Punkt P3 und ihren Schnittpunkt (x(1) , f1∗ )
mit der Gerade x = x(1) festgelegt.
2.7 Einschlussverfahren 79
Wenn g > 0 ist, liegen die Punkte P1 = (x(1) , f1 ) und (x(1) , f1∗ ) auf derselben Seite der
x-Achse, und der Schnittpunkt x(4) der Parabeltangente mit der x-Achse liegt zwischen
x(3) und x(1) .
Wegen des vorhergehenden Sekantenschritts sind die Punkte (x(1) , f1 ), (x(2) , f2 ) und
(x(3) , 0) kollinear, und daher gilt
f1 − f2 f2
s12 = = (2) .
x(1) − x(2) x − x(3)
Damit ergibt sich für g wesentlich einfacher
Abb. 2.19. Zum Verfahren von Anderson-Björck. Mit den Geraden P2 P3 und P1 P3
ergeben sich die Punkte H1 und H2 . Die zu H1 H2 parallele Gerade durch P3 ist die
Parabeltangente in diesem Punkt
80 2. Lösung nichtlinearer Gleichungen
Beispiel 2.39.
Gegeben: Die Funktion f (x) = sin x + 1 − 1/x und die Schranke RelErr = 5 · 10−7 für
den relativen Fehler (AbsErr = 0).
Gesucht: Die Nullstelle ξ von f im Intervall (0.6, 0.7) mit Hilfe des Verfahrens von
Anderson–Björck.
Lösung:
Zur Darstellung dieser Tabelle vgl. Beispiel 2.33. Wegen |x(2) − x(1) |/|x(2) | < 5 · 10−7 für
ν = 4 ist die Nullstelle auf 7 Ziffern genau, also ξ ≈ x(1) = 0.629 446 5 .
Das Verfahren von King (s. Originalarbeit [KING1973]) unterscheidet sich vom Pegasus-
Verfahren nur dadurch, dass nie zwei Sekantenschritte nacheinander ausgeführt werden,
sondern auf jeden Sekantenschritt ein modifizierter Schritt folgt. Es besitzt eine etwas
höhere Konvergenzordnung (siehe Abschnitt 2.9).
Das Verfahren von Anderson-Björck-King verläuft ganz analog; es arbeitet nach dem
Verfahren von Anderson-Björck mit der Zusatzbedingung von King, dass nie zwei Sekan-
tenschritte nacheinander erfolgen dürfen. Hier geschieht der modifizierte Schritt nach der
Anderson-Björck-Methode, beim Verfahren von King nach der Pegasus-Methode.
Gemeinsam ist den Einschlussverfahren Pegasus, Anderson-Björck und damit auch King
der Sekantenschritt. Lediglich der modifizierte Schritt wird unterschiedlich realisiert; bei
Pegasus wird g = f2 /(f2 + f3 ) gesetzt, bei Anderson-Björck g = 1 − f3 /f2 und, falls g ≤ 0
ist, wird g = 0.5 gesetzt. Das Illinois-Verfahren verwendet stets g = 0.5.
2.7 Einschlussverfahren 81
In dem folgenden Algorithmus 2.40 wird das gegebene Einschlussintervall [a, b] zunächst
mit dem Bisektionsverfahren so lange verkleinert, bis seine Länge nicht größer ist als
eine vorgegebene Länge LB. Erst dann kommt ein Einschlussverfahren höherer Konver-
genzordnung, das Pegasus-Verfahren, das Verfahren von Anderson-Björck oder auch das
Illinois-Verfahren, zum Einsatz.
Die Formulierung der Algorithmen 2.36 und 2.38 ist günstig für die Herleitung und die
geometrische Interpretation der Verfahren. Der folgende Algorithmus führt dagegen einen
eventuell erforderlichen modifizierten Schritt erst nach einer nicht erfüllten Abbruchbe-
dingung aus; daher weichen die Formeln für g von denen in Abschnitt 2.7.4 und 2.7.5 ab.
Der Algorithmus 2.40 kann in etwas erweiterter Form auch mit den Verfahren von King
und Anderson-Björck-King arbeiten. Zugunsten einer einfachen und übersichtlichen Dar-
stellung wird auf diese Verfahren verzichtet (siehe jedoch Beispiel 2.41).
Algorithmus 2.40.
Gegeben: f ∈ C[a, b] mit f (a) · f (b) < 0, Schranken AbsErr und RelErr für den
absoluten bzw. den relativen Fehler mit AbsErr > 0 und RelErr = 0
oder mit AbsErr = 0 und RelErr > 0, die Länge LB > 0 für die
Entscheidung zwischen Bisektions- und Sekantenschritt, die maximal
zulässige Anzahl nf max der Funktionsauswertungen.
Gesucht: Eine Zahl ξ ∈ (a, b), für die f (ξ) = 0 ist, oder ein Einschlussintervall
[x(1) , x(2) ] bzw. [x(2) , x(1) ] für ξ mit
|x(1) − x(2) | ≤ |x(2) | RelErr + AbsErr .
Vorbereitung: x(1) := a, x(2) := b, f1 := f (x(1) ), f2 := f (x(2) ),
nf := 2, v := x(1) − x(2) .
1. Wenn nf ≥ nf max ist, konnte mit nf max Funktionsauswertungen eine Nullstelle
von f nicht ermittelt werden. Die Iteration wird dann abgebrochen.
6. Prüfung der Abbruchbedingung: Wenn |v| ≤ tol ist, wird die Iteration abgebro-
chen mit
ξ := x(2) , falls |f2 | ≤ |f1 | ist , andernfalls mit ξ := x(1) .
Wenn |v| > tol ist, weiter mit 7.
Bemerkungen.
• Für das kombinierte Verfahren eignet sich eine Länge LB mit 0 < LB ≤ 0.2,
beispielsweise LB = 0.15. Mit LB > |b − a| werden nur Sekantenschritte, mit
LB = 0 nur Bisektionsschritte ausgeführt. Siehe dazu auch Beispiel 2.41.
• Zu 1. Die Vorgabe von nf max (etwa nf max = 100) schützt vor einer Endlos-
schleife, wenn AbsErr oder RelErr zu klein oder wenn LB zu groß gewählt worden
ist.
• Zu 5. Wenn e < 0 ist, liegt ξ zwischen x(2) und x(3) , andernfalls zwischen x(1) und
x(3) .
• Zu 7. Wegen der Vertauschung von f2 mit f3 in 5. sehen die Formeln für g
anders aus als in den Abschnitten 2.7.4 und 2.7.5. Ein Anderson-Björck-Schritt
mit g = 1 − f2 /f3 setzt voraus, dass vorher ein Sekantenschritt (Bis = 0) aus-
geführt wurde (siehe Abschnitt 2.7.5). Deshalb wird er nach einem Bisektionsschritt
(Bis = 1) durch einen Pegasusschritt ersetzt.
2.7 Einschlussverfahren 83
Beispiel 2.41.
Für einen Vergleich verschiedener Einschlussverfahren werden üblicherweise die Nullstel-
len von gewissen Testfunktionen berechnet und die dafür benötigten Funktionsauswer-
tungen gezählt.
Für 12 solcher Funktionen (als eine Auswahl aus einem umfangreicheren Test) sind die
Ergebnisse in der nachfolgenden Tabelle zusammengestellt. Verglichen werden das Bisek-
tionsverfahren, das Zeroin-Verfahren, das Illinois- und das Pegasus-Verfahren sowie die
Verfahren von Anderson-Björck, King und Anderson-Björck-King.
Für die letzten fünf Verfahren wird der Algorithmus 2.40 (mit einer Erweiterung für die
letzten beiden Verfahren) angewendet. Sie werden jeweils einmal ohne eine vorangehende
Verkleinerung des Startintervalls [a, b] (LB > |b − a|) und einmal mit einer solchen Ver-
kleinerung mittels Bisektion (LB < |b−a|) eingesetzt. Mit LB = 0 führt der Algorithmus
das Bisektionsverfahren aus.
Damit beim Test alle Verfahren mit einem Intervall der Länge LB beginnen, haben die
Startintervalle die Längen 4 · LB, 8 · LB oder 16 · LB, so dass nach 2, 3 bzw. 4 Bisekti-
onsschritten die Länge LB erreicht ist. Für die Tabelle ist LB = 0.15 gewählt.
Die Tabelle zeigt, dass bei der Mehrzahl der Funktionen eine vorbereitende Intervallver-
kleinerung zweckmäßig oder sogar unerlässlich ist. — bedeutet, dass die Nullstelle mit
100 Funktionsauswertungen nicht gefunden wurde.
Verfahren/Funktion 1 2 3 4 5 6 7 8 9 10 11 12
Bisektion 40 39 40 41 41 40 41 42 40 38 40 12
Illinois ohne 14 21 63 18 22 22 14 18 19 13 28 25
Bisektion
Illinois mit 12 13 19 14 14 14 13 15 13 12 15 21
Bisektion
Pegasus ohne 12 19 63 17 20 28 11 18 18 12 26 36
Bisektion
Pegasus mit 11 12 16 12 12 11 11 16 11 10 14 29
Bisektion
King ohne 11 17 60 16 19 28 12 18 17 10 25 38
Bisektion
King mit 11 12 15 12 12 11 11 15 11 10 13 29
Bisektion
Anderson-Björck 14 74 — 29 11 22 14 7 19 12 — 30
ohne Bisektion
Anderson-Björck 10 11 16 11 11 11 12 11 12 10 14 24
mit Bisektion
Anderson-Björck- 16 67 — 26 11 22 13 7 18 11 — 29
King ohne
Bisektion
Anderson-Björck- 10 11 15 11 11 11 11 10 11 10 13 24
King mit
Bisektion
Zeroin 12 14 17 10 11 11 13 13 15 12 14 28
2.8 Anwendungsbeispiele 85
2.8 Anwendungsbeispiele
Eingabeparameter:
λ : Rohrreibungsbeiwert
Zu den Eingabeparametern gehören noch Größen aus den Abbruchkriterien für die Ite-
ration, z. B.
Ausgabeparameter:
d : Rohrdurchmesser in m.
............. ...
......................
...................................................................................................
..............................................................................
6
.....................................................
.......
Düse am Rohrleitungsende:
. . . . . . . ... . . ................
....................................................................................
...........................................................
.............................................................. ..
.......................................................................... ..................
..................... ... .................................. ........................................
................................................................................ .... ........
.... ....
6
............
................................................................................................................................... ... ...
... ... .. ...
.......
............
............................ . ...............................
..................................................................................................................................................................... ... ... ... ...
..........
......................................................................................
... ...
... ... d ... ...
. .
....
. 6
.. ... ...
... ... ..... ..... ..... ....... ..... ..... ..... ..... ..... ..... ..... .......
. .
dst
... ... .. .....
...
...
... ...
... ...
... ...
... ...
... ...
. .. .......
. .. ...
.............
......
...
?
....................
.... ... ... .
?
. ....... ...
.. .. .. ...... ... ......
... ...
... .... .... .................................. ...................................
. .
... .... .... ... .................. A
-d
. st
.... . . ... ..
H .. .... ..... AR
... ... ...
... ...
...
dst 2
.... ... ...
.. ..... ..... Turbinenlaufrad
... ... ... Ast
...
...
... ...
.... .... ..
..... .... ......
α= AR = d
...... . ..
... ..
...
a
... ...
... ... ... ....... ..... .......
... ... ...
............ ......
? ... ... ...
. .
...... ........ ...... ...... ...... ...... ...... ......... .......... ...... ...... .................................................................... ..... ..... ...
... ...
. ....
................... .... .
.
.
.
... .
.
... ... .... .... ...
... ... ... ... ..
.... Düse
... ......................................... ...
..............................................
.
Abb. 2.20.
86 2. Lösung nichtlinearer Gleichungen
Vorgehensweise:
(a) Zunächst berechnet man das optimale Flächenverhältnis der Düse αopt aus der
Forderung, dass die Strahlleistung Pst maximal wird.
(b) Für dieses optimierte Flächenverhältnis αopt ist dann derjenige Rohrdurchmesser
d zu bestimmen, der für den vorgegebenen Volumenstrom erforderlich ist.
vR = α vst .
$ %3/2
2gH
= 12 AR α 1+ζD +α2 (λ L
.
d +ζR )
∂Pst
Zu (a): Mit der notwendigen Bedingung ∂α = 0 für eine maximale Strahlleistung ergibt
sich das optimale Flächenverhältnis
#
1 + ζD
αopt = .
2(λ Ld + ζR )
∂ 2 Pst
(Wegen ∂α2 (αopt ) < 0 liegt tatsächlich ein Maximum vor.)
Setzt man αopt für α in die Formel für die Strahlgeschwindigkeit vst ein, so ergibt sich
die optimale Strahlgeschwindigkeit
#
4gH
vst (αopt ) = .
3(1 + ζD )
2.8 Anwendungsbeispiele 87
d2
AR = π
4
ist nach (2.57) der optimale Volumenstrom
#
d2 d2 2gH
V̇ = π vR (αopt ) = π ,
4 4 3(λ Ld + ζR )
und daraus folgt für den optimalen Rohrdurchmesser d die nichtlineare Gleichung
#
4V̇ 3 L
d= λ + ζR = ϕ(d) . (2.58)
π 2gH d
24V̇ 2 L
f (d) = d 4 − (λ + ζR ) = 0 . (2.59)
π 2 gH d
Zahlenbeispiel.
V̇ = 6.2 m3 /s, H = 1130 m, L = 1300 m, ζD = 0.04, ζR = 2.5, g = 9.81 m/s2 ,
= 103 kg/m3 , λ = 0.02.
Schranke für den absoluten Fehler ε = 0.5 · 10−6 .
Mit den Startwerten d(1) = 0.7 und d(2) = 0.8 sind f (0.7) < 0 und f (0.8) > 0; also ist
0.7 < d < 0.8. Für die Lösung der Gleichung (2.59)
26
d − 0.008 432 327
4
+ 2.5 = 0
d
liefert das Pegasus-Verfahren mit 5 Iterationsschritten
d = 0.748 551.
Somit ergibt sich für den gesuchten optimalen Rohrdurchmesser bei Rundung auf
3-stellige Mantisse d ≈ 0.749 m.
Mit dem Quadrat der optimalen Strahlgeschwindigkeit
2 4gH 4 · 9.81 · 1130 m2
vst (αopt ) = = = 14211.9 m2/s2
3(1 + ζD ) 3 · 1.04 s2
erhält man die maximale Strahlleistung
1 2
Pst = V̇ vst (αopt ) = 44 057 kW.
2
88 2. Lösung nichtlinearer Gleichungen
Abb. 2.21.
2P A2
f (V̇ ) := V̇ 3 − 2ghA2 V̇ + = 0.
2.9 Effizienz der Verfahren und Entscheidungshilfen 89
Die praktisch interessierende Lösung dieser Gleichung f (V̇ ) = 0 kann z. B. mit dem
Newton-Verfahren ermittelt werden. Ein Startwert V̇ (0) für die Iteration ergibt sich aus
(2.60) bei Vernachlässigung des Austrittsverlustes und anschließender Auflösung nach V̇
zu
P
V̇ (0) =
·h·g
mit = 103 kg/m3 , g = 9.81 m/s2 . Man kann aber auch mit jedem Einschlussverfahren
arbeiten.
Zur Lösung der kubischen Gleichung f (V̇ ) = 0 könnte aber ebenso ein Verfahren zur Be-
stimmung sämtlicher Lösungen einer algebraischen Gleichung ohne Kenntnis von Start-
werten (z. B. das Verfahren von Muller) verwendet werden (s. Abschnitt 3.3); hier ist
dann aus den drei Lösungen die für den Anwendungsfall sinnvolle Lösung auszuwählen
(siehe dazu das Beispiel 11).
Die folgende Tabelle gibt eine Übersicht über die Konvergenzordnung und den Effizienz-
index bei Verfahren zur Berechnung einfacher und mehrfacher Nullstellen; je größer E,
desto effizienter ist das Verfahren in der Umgebung der Nullstelle.
Es zeigt sich, dass man (nicht zuletzt wegen der sicheren Konvergenz) am effektivsten
mit dem Pegasus-Verfahren, dem Verfahren von Anderson-Björck oder den Verfahren von
King und Anderson-Björck-King arbeitet. Das Bisektionsverfahren wird benutzt, um ein
Startintervall mit Einschluss zu verkleinern, bevor eines der genannten Verfahren einge-
setzt wird.
Will man einzelne Polynom-Nullstellen berechnen, so lässt sich auch das Newton-Ver-
fahren zusammen mit dem Hornerschema zur Berechnung der Funktions- und Ablei-
tungswerte effektiv einsetzen. Für transzendente Gleichungen sind weder das Sekanten-
verfahren noch das Newton-Verfahren zu empfehlen.
Die Ermittlung der Lösung einer Gleichung f (x) = 0 mit dem Algorithmus 2.17 für das
allgemeine Iterationsverfahren erfordert umfangreiche Vorbereitungen. Der größte Auf-
wand steckt im 3. Schritt, in dem zu prüfen ist, ob bei der zu f (x) = 0 äquivalenten
Gleichung ϕ(x) = x die Funktion ϕ den Voraussetzungen des Existenz- und Eindeu-
tigkeitssatzes 2.11 genügt. Dagegen kann schon nach dem 1. Schritt – Festlegung eines
Intervalls I, in dem mindestens eine Nullstelle von f liegt – ein Einschlussverfahren zur
Bestimmung dieser Nullstelle eingesetzt werden. Darum hat das allgemeine Iterations-
verfahren für die praktische Ermittlung von Nullstellen keine Bedeutung.
Für den Fall mehrfacher Nullstellen sind die Anmerkungen zu mehrfachen Nullstellen“
”
in Abschnitt 2.5.3 zu beachten.
3.1 Vorbemerkungen
n
Pn (x) = aj xj , aj ∈ C, n ∈ IN, an = 0 (3.1)
j=0
vom Grad n betrachtet. Gesucht sind Nullstellen eines Polynoms (3.1) und somit Lösun-
gen einer algebraischen Gleichung n-ten Grades
n
Pn (x) = aj xj = 0. (3.2)
j=0
Der Fundamentalsatz der Algebra besagt, dass eine algebraische Gleichung (3.2) genau
n komplexe Lösungen xk besitzt, die entsprechend ihrer Vielfachheit αk gezählt werden.
Jedes algebraische Polynom (3.1) lässt sich daher in n Linearfaktoren zerlegen:
Pn (x) = an (x − x1 )(x − x2 ) . . . (x − xn ).
Kommt der Linearfaktor (x − xk ) genau αk -fach vor, so heißt xk eine αk -fache Lösung
von (3.2), und es ergibt sich die Zerlegung
Wenn alle Koeffizienten des Polynoms (3.1) reell sind, aj ∈ IR, können komplexe Lösun-
gen von (3.2) nur als Paare konjugiert komplexer Lösungen auftreten, d. h. mit x = α+iβ
ist auch x = α − iβ eine Lösung von (3.2), und zwar mit derselben Vielfachheit. Der Grad
n einer Gleichung (3.2) mit reellen Koeffizienten, die keine reellen Wurzeln besitzt, muss
somit gerade sein, und jede derartige Gleichung ungeraden Grades besitzt mindestens
eine reelle Lösung.
92 3. Verfahren zur Lösung algebraischer Gleichungen
und das Polynom Pn−1 vom Grad n−1 besitzt die n−1 restlichen Nullstellen von Pn .
Pn−1 ergibt sich mittels Division von Pn durch x − x0
Pn (x)
Pn−1 (x) =
x − x0
und heißt daher abdividiertes Polynom oder Deflationspolynom. Wenn x0 keine Nullstelle
von Pn ist, ist Pn nicht ohne Rest durch x − x0 teilbar, und aus
Außerdem gelten
Pn (x) − Pn (x0 )
Pn−1 (x) = , x = x0 ,
x − x0
und
Pn (x) − Pn (x0 )
Pn−1 (x0 ) = lim = Pn (x0 ).
x→x0 x − x0
Zur Berechnung von Funktions- und Ableitungswerten eines Polynoms und zur Berech-
nung von abdividierten Polynomen bzw. Deflationspolynomen dient das im Folgenden
erläuterte Horner-Schema. Es wird bei allen Verfahren zur Berechnung von Polynomnull-
stellen als Hilfsmittel eingesetzt, beispielsweise wenn Polynomwerte mit den in Kapitel 2
behandelten Newton-Verfahren ermittelt werden.
Will man sämtliche Lösungen einer algebraischen Gleichung (3.2) berechnen, so ver-
wendet man z. B. das Muller-Verfahren, das Verfahren von Bauhuber oder das Verfahren
von Jenkins und Traub (siehe Abschnitt 3.3).
In dieser Darstellung sind außer Additionen nur Multiplikationen mit dem festen, reellen
Faktor x0 auszuführen. Im Einzelnen sind zu berechnen:
⎧ (1)
⎪
⎪ an := an
⎪
⎪
⎪
⎪
(1)
an−1
(1)
:= an x0 + an−1
⎪
⎪
⎨ a(1) (1)
:= an−1 x0 + an−2
n−2
.. (3.3)
⎪
⎪
⎪
⎪
.
⎪
⎪ (1) (1)
⎪
⎪ a := a2 x0 + a1
⎩ 1(1) (1)
a0 := a1 x0 + a0 = Pn (x0 ).
In der ersten Zeile stehen, geordnet nach absteigenden Potenzen von x, alle Koeffizienten
(auch die mit dem Wert 0) des Polynoms Pn . Die Summe der untereinander stehenden
Elemente der ersten und der zweiten Zeile wird in der dritten Zeile notiert. Dann wird
diese Summe mit x0 multipliziert und das Produkt in der zweiten Zeile rechts daneben
aufgeschrieben usw.
Bei komplexen Koeffizienten sind sowohl der Realteil als auch der Imaginärteil mit x0 zu
multiplizieren:
(1) (1) (1) (1) (1)
aj x0 = (αj + i βj ) x0 = αj x0 + i βj x0 .
94 3. Verfahren zur Lösung algebraischer Gleichungen
Beispiel 3.1.
Gegeben: Das Polynom P3 : P3 (x) = 2x3 − 4x2 + 3x + 15 mit reellen Koeffizienten.
P3 2 −4 3 15
x0 = 2 0 4 0 6
2 0 3 21 = P3 (2)
Bei der Division von Pn (x) durch (x − x0 ) mit x0 ∈ IR entsteht als Quotient ein Polynom
Pn−1 (x) vom Grad n−1 und als Rest eine Zahl Rn , dividiert durch (x − x0 ):
Pn (x) Rn
= Pn−1 (x) + , x = x0 , bzw.
x − x0 x − x0
Pn (x) = (x − x0 )Pn−1 (x) + Rn . (3.4)
Setzt man dies in (3.5) ein, so ergibt sich unter Verwendung von (3.3):
(1) (1) (1) (1) (1)
(x−x0 )Pn−1 (x)+Pn (x0 ) = an xn−1 + an−1 xn−2 +. . . + a2 x + a1 (x − x0 ) + a0
(1) (1) (1) (1) (1)
= an xn + an−1 − an x0 xn−1 +. . .+ ak − ak+1 x0 xk
(1) (1) (1) (1)
+. . .+ a1 − a2 x0 x + a0 − a1 x0 = Pn (x).
3.2 Das Horner-Schema 95
Besitzt das Polynom Pn komplexe Koeffizienten und ist x0 ein komplexer Argumentwert,
so kann man zur Berechnung des Funktionswertes Pn (x0 ) das einfache Horner-Schema
verwenden. Man hat dann lediglich für jeden Koeffizienten eine reelle und eine imaginäre
Spalte zu berechnen; siehe Beispiel 3.1.
Besitzt das Polynom Pn jedoch nur reelle Koeffizienten, so kann man zur Berechnung
des Funktionswertes Pn (x0 ) zu einem komplexen Argumentwert x0 mit dem Horner-
Schema ganz im Reellen bleiben, wenn man das sogenannte doppelreihige Horner-Schema
verwendet. Zunächst nimmt man den zu x0 konjugiert komplexen Argumentwert x0 hinzu
und bildet
(x − x0 )(x − x0 ) = x2 − px − q
mit reellen Zahlen p und q; es gilt p = x0 + x0 , q = −x0 · x0 .
Dividiert man jetzt Pn durch (x2 − px − q), so erhält man die Beziehung
⎧
⎨ Pn (x) (1)
= (x2 − px − q)Pn−2 (x) + b1 x0 + b0
(1)
mit
(3.6)
⎩ P (1) (1) (1)
bn xn−2 + bn−1 xn−3 + . . . + b3 x + b2 .
(1)
n−2 (x) =
(1)
Für die Koeffizienten bk von Pn−2 gelten
⎧ (1)
⎪
⎪ bn = an ,
(0)
⎪
⎪
⎪
⎪
⎪
⎪
(1)
bn−1
(0)
= an−1 + pbn ,
(1)
⎪
⎪
⎪
⎨ ..
. (3.7)
⎪
⎪ bk
(1) (0) (1) (1)
⎪
⎪ = ak + pbk+1 + qbk+2 , k = (n−2)(−1)1,
⎪
⎪
⎪
⎪ ..
⎪
⎪ .
⎪
⎩ (1) (0) (1)
b0 = a0 + qb2 .
Beispiel 3.2.
Gegeben: P5 (x) = x5 − 3x4 − x + 3.
Lösung: Zunächst berechnet man die Werte p und q aus der Beziehung
p = x0 + x0
q = −x0 x0 .
Daraus ergeben sich: p = i − i = 0 und q = (−i)(−i) = −1.
P5 1 −3 0 0 −1 3
q = −1 0 0 (−1)1 (−1)(−3) (−1)(−1) (−1)3
p=0 0 0·1 0·(−3) 0·(−1) 0·3 0·0
1 −3 −1 3 0 0
(1)
Da das Horner-Schema neben dem Funktionswert Pn (x0 ) auch die Koeffizienten aj des
abdividierten Polynoms Pn−1 liefert, ergibt sich die Möglichkeit, die k-ten Ableitungen
(k)
Pn des Polynoms Pn für k = 1(1)n an der Stelle x0 ∈ IR zu berechnen. Aus (3.5) ergibt
sich für die 1. Ableitung
Pn (x0 ) erhält man also, indem man an die 3. Zeile des Horner-Schemas ein weiteres ein-
faches Horner-Schema anschließt.
P3 ...
x0 = 2 ...
P2 2 0 3 21 = P3 (2)
x0 = 2 0 4 8
Verfährt man nun mit Pn−1 (x) analog wie mit Pn (x), so erhält man
1
Pn−1 (x0 ) = Pn−2 (x0 ) = P (x0 ).
2 n
So fortfahrend erhält man schließlich für die abdividierten Polynome
1 (k)
Pn−k (x0 ) = P (x0 ) für k = 1(1)n.
k! n
Mit
ck := Pn−k (x0 ) (3.10)
ergibt sich durch sukzessives Einsetzen der Gleichungen in (3.9) für Pn die Darstellung
98 3. Verfahren zur Lösung algebraischer Gleichungen
(n+1) 1 (n)
an = n! Pn (x0 ) = P0 (x0 )
Die Aufstellung der Taylorentwicklung von Pn an einer Stelle x0 mit Hilfe des vollständi-
gen Horner-Schemas erfordert (n2 + n)/2 Punktoperationen, während der übliche Weg
(Differenzieren, Berechnen der Werte der Ableitungen, Dividieren durch k!, wobei k! als
bekannter Wert vorausgesetzt wird) n2 + 2n − 2, also für n ≥ 3 mehr als doppelt so
viele Punktoperationen erfordert. Durch das Einsparen von Punktoperationen wird das
rundungsfehlergünstige Arbeiten ermöglicht, denn durch hohe Potenzen häufen sich sy-
stematische Rundungsfehler an.
Lösung: Die Vorgehensweise erfolgt laut Rechenschema 3.4. Da nur die Koeffizienten
der Taylorentwicklung interessieren, werden die Bezeichnungen Pn−k der Po-
lynome niedrigeren Grades weggelassen.
P3 2 −4 3 15
x0 = 2 0 4 0 6
P2 2 0 3 21 = P3 (2)
x0 = 2 0 4 8
P1 2 4 11 = P3 (2)
x0 = 2 0 4
1
P0 2 8= 2! P3 (2)
x0 = 2 0
1
2= 3! P3 (2)
100 3. Verfahren zur Lösung algebraischer Gleichungen
3.2.4 Anwendungen
Das Horner-Schema wird verwendet
(1) zur bequemen, schnellen und rundungsfehlergünstigen Berechnung der Funktions-
werte und Ableitungswerte eines Polynoms Pn
(2) zur Aufstellung der Taylorentwicklung eines Polynoms Pn an einer Stelle x0
(3) zum Abdividieren von Nullstellen (Deflation von Polynomen)
Man wird z. B. bei der iterativen Bestimmung einer Nullstelle eines Polynoms nach einem
der Newton-Verfahren (siehe Abschnitt 2.5) Pn , Pn bzw. Pn , Pn , Pn mit dem Horner-
Schema berechnen.
Hat man für eine Nullstelle x1 von Pn iterativ eine hinreichend gute Näherung erhalten,
so dividiert man Pn durch (x − x1 ) und wendet das Iterationsverfahren auf das Deflati-
onspolynom Pn−1 an. So erhält man nacheinander alle Nullstellen von Pn und schließt
aus, eine Nullstelle zweimal zu berechnen. Dabei könnten sich aber die Nullstellen der
abdividierten Polynome immer weiter von den Nullstellen des Ausgangspolynoms Pn ent-
fernen, so dass die Genauigkeit mehr und mehr abnimmt. Wilkinson empfiehlt deshalb
in [WILK1969], S.70-83, das Abdividieren von Nullstellen grundsätzlich mit der betrags-
kleinsten Nullstelle zu beginnen, d. h. mit einer Methode zu arbeiten, die für das jeweilige
Polynom eine Anfangsnäherung so auswählt, dass die Iteration gegen die betragskleinste
Nullstelle konvergiert (s. Verfahren von Muller, Abschnitt 3.3.2). Wird diese Forderung
erfüllt, so ergeben sich alle Nullstellen mit einer Genauigkeit, die im Wesentlichen von
ihrer Kondition bestimmt ist und nicht von der Genauigkeit der vorher bestimmten
Nullstelle. Wilkinson empfiehlt außerdem, nachdem man alle Nullstellen mittels Iterati-
on und Abdividieren gefunden hat, die berechneten Näherungswerte als Startwerte für
eine Nachiteration mit dem ursprünglichen Polynom zu verwenden. Man erreicht damit
eine Erhöhung der Genauigkeit, besonders in den Fällen, in denen das Abdividieren die
Kondition verschlechtert hat.
Beispiel 3.6.
Gegeben: Das Polynom P3 : P3 (x) = x3 − 1 und die Nullstelle x1 = 1 von P3 .
Lösung: P3 1 0 0 −1
x1 = 1 0 1 1 1
P2 1 1 1 0 = P3 (1)
Wenn hinreichend genaue Anfangsnäherungen für die Nullstellen eines Polynoms vorlie-
gen, kann man mit Iterationsverfahren Folgen von Näherungswerten konstruieren, die
gegen die Nullstellen konvergieren. Das Problem liegt in der Beschaffung der Startwerte.
Will man z. B. sämtliche reellen Nullstellen eines Polynoms Pn mit reellen Koeffizienten
mit Hilfe eines der bisher angegebenen Iterationsverfahren berechnen, so muss man:
1. ein Intervall I ermitteln, in dem alle Nullstellen liegen. Das kann z. B. nach dem
folgenden Satz geschehen:
Ist Pn (x) = xn + an−1 xn−1 + . . . + a1 x + a0 das gegebene Polynom und
A = max |ak |, so liegen alle Nullstellen von Pn in einem Kreis um den Null-
k=0(1)n−1
punkt der komplexen Zahlenebene mit dem Radius r = A+1. Also ist I = [−r, +r].
Ist Pn ein Polynom mit lauter reellen Nullstellen (z. B. ein Orthogonalpolynom, s.
Abschnitt 8.2.2, Sonderfälle 2.), so kann der Satz von Laguerre angewandt werden:
Die Nullstellen liegen alle in einem Intervall, dessen Endpunkte durch die beiden
Lösungen der quadratischen Gleichung
gegeben sind.
2. die Anzahl reeller Nullstellen nach den Vorzeichenregeln von Sturm und Descartes
berechnen,
3. die Lage der Nullstellen durch Intervallteilung, Berechnung der Funktionswerte und
Abzählung der Anzahl der Vorzeichenwechsel ermitteln.
Mit 3. ist es möglich, Intervalle Ik ⊂ I anzugeben, in denen jeweils nur eine Nullstelle
xk ungerader Ordnung liegt. Dann lässt sich z. B. das Newtonsche Verfahren zur nähe-
rungsweisen Berechnung der xk anwenden. Dabei sind Pn und Pn (bzw. Pn , Pn , Pn ) mit
Hilfe des Horner-Schemas zu berechnen.
Der soeben beschriebene Weg ist mühsam und für die Praxis uninteressant. Hier braucht
man Verfahren, die in kürzester Zeit und ohne Kenntnis von Startwerten sämtliche reel-
len und komplexen Nullstellen eines Polynoms mit reellen bzw. komplexen Koeffizienten
liefern.
Für Polynome mit reellen Koeffizienten werden diese Anforderungen mühelos vom Ver-
fahren von Muller erfüllt, s. Abschnitt 3.3.2. Das Muller-Verfahren lässt sich auch auf
Polynome mit komplexen Koeffizienten erweitern.
102 3. Verfahren zur Lösung algebraischer Gleichungen
Für Polynome mit komplexen Koeffizienten werden hier zwei Verfahren genannt, das
Verfahren von Jenkins und Traub und das Verfahren von Bauhuber. Das Verfahren von
Jenkins-Traub wird hier nur kurz ohne Formulierung eines Algorithmus beschrieben.
Das Verfahren von Muller [MULL1956] liefert ohne vorherige Kenntnis von Startwerten
sämtliche reellen und konjugiert komplexen Nullstellen eines Polynoms
n
Pn : Pn (x) = a0 + a1 x + a2 x2 + . . . + an xn = aj xj für aj ∈ IR, an = 0.
j=0
Hat man durch Abdividieren und Anwendung des Muller-Verfahrens auf das jeweilige
Deflationspolynom alle Nullstellen von Pn näherungsweise gefunden, so empfiehlt es sich
grundsätzlich, die gewonnenen Näherungswerte als Startwerte für eine Nachbesserung
mit dem Newton-Verfahren, angewandt auf das ursprüngliche Polynom Pn , zu verwen-
den. Man sollte aber auf jeden Fall erst sämtliche Nullstellen von Pn auf dem beschrie-
benen Weg näherungsweise berechnen, bevor man sie verbessert. Nach Untersuchungen
von Wilkinson ist es nicht notwendig, direkt jede Nullstelle noch vor dem Abdividieren
mit dem ursprünglichen Polynom zu verbessern.
verwendet. In einem
nächsten Schritt wird eine neue Parabel y = Φ(x) durch die Punkte
x(k) , Pn (x(k) ) für k = 1, 2, 3 gelegt und die zu x(3) nächstliegende Nullstelle von Φ als
neue Näherung x(4) für x1 bestimmt. Analog fortfahrend bricht man das Verfahren ab,
wenn sich für einen festen Index N die Näherungen x(N ) und x(N −1) hinreichend wenig
voneinander unterscheiden, x(N ) ist dann der gesuchte Näherungswert für die Nullstelle
x1 von Pn .
y6
...........
.... ........................
... .....
.. ...... b
....
b
....
........... .....
..... ...
.
..
........... ....
.....
... . .
. ...
..
.
... .. ..
.. .
.... ... ...
... ..
... .
... ..
...
..
y = Pn (x)
..
.. .
.... ... . ... ... .
...
.. ... .. .. ... ... ...
.. ... .. ... ... ..
...
.. ... .. ... .. ..
..
.. ... .. .. .. .. ..
.. .... . .
. .. ..
... .. .. (0)... (2) ... .. (1) ...
...
...
..
.. ..
..
x r x
.. r .. ..
... . . rx ... -
.. ... .. .. ..
... ... .. ... .... .. ...
−2 −1
.
...
...
....
.
...
... .
.. ... . .. .
1
... .. . ... 2 .
...
. x
..... .... .. ...... .. .. ...
........ ........... .. .. .. . ..
.......
.. .... ..... .....
....
..
..
b .. ........
. ....
(3) ... ....................................
..
.. x ... .
..
.. ..
.. ..
.. ..
.. ..
.. y = Φ(x)
..
Um die Nullstelle x2 von Pn zu erhalten, wird der Linearfaktor (x − x(N ) ) von Pn mit
Hilfe des Horner-Schemas abdividiert, der Rest vernachlässigt und das Muller-Verfahren
nunmehr auf Pn−1 angewandt. Man erhält so im Allgemeinen die betragskleinste Null-
stelle von Pn−1 , d. h. die Nullstelle x2 von Pn mit |x1 | ≤ |x2 |.
Mit hν := x(ν) − x(ν−1) , h := x − x(ν) erhält man für (3.13) die Darstellung
(h + hν )(h + hν + hν−1 )
Φ(x) = Φ(x(ν) + h) = f (x(ν) )
hν (hν + hν−1 )
h(h + hν + hν−1 )
− f (x(ν−1) )
hν hν−1
h(h + hν )
+ f (x(ν−2) ).
(hν + hν−1 )hν−1
104 3. Verfahren zur Lösung algebraischer Gleichungen
Zur Bestimmung der Schnittpunkte der Parabel mit der x-Achse setzt man
Aν q 2 + Bν q + Cν = 0.
und damit auf die Lösungen x1/2 der Gleichung Φ(x(ν) + qhν ) = 0 :
Eine der beiden Lösungen wird als neue Näherung x(ν+1) verwendet, es wird diejenige
ausgewählt, die näher an x(ν) liegt.
" Dies wird dadurch realisiert, dass im Nenner von
(3.15) das Vorzeichen der Wurzel Bν2 − 4Aν Cν so gewählt wird, dass der Nenner den
größeren Betrag erhält.
Falls der Nenner von (3.15) verschwindet (dies ist dann der Fall, wenn
f (x(ν) ) = f (x(ν−1) ) = f (x(ν−2) ) gilt), schlägt Muller vor, statt (3.15) für q1/2 = 1
zu setzen und damit weiterzurechnen. In diesem Fall entartet die Parabel in eine zur
x-Achse parallele Gerade, die keinen Schnittpunkt mit der x-Achse hat.
Der ausgewählte Wert von q1/2 wird dann qν+1 genannt, und es wird gesetzt
Automatischer Startprozess
Als Startwerte für die Iteration werden fest vorgegeben
Als Funktionswerte an den Stellen x(0) , x(1) , x(2) (und nur an diesen!) werden im All-
gemeinen nicht die Funktionswerte des jeweiligen Polynoms Pn genommen, sondern die
Werte
a0 − a1 + a2 statt f0 := Pn (x(0) ),
a0 + a1 + a2 statt f1 := Pn (x(1) ), (3.19)
a0 für f2 := Pn (x(2) ).
Die Verwendung dieser künstlichen Funktionswerte wurde von Muller selbst empfohlen
(vgl. Algorithmus 3.7). Man kann aber ebenso an den Startstellen x(0) , x(1) , x(2) die wirk-
lichen Polynomwerte f0 := Pn (x(0) ), f1 := Pn (x(1) ), f2 := Pn (x(2) ) benutzen, wie es in
Abbildung 3.1 gemacht wurde.
Abbruchbedingung
Die Iteration (3.17) wird abgebrochen, falls zu vorgegebenem ε > 0 die Bedingung
(ν+1)
x − x(ν)
x(ν+1) <ε (3.20)
(N )
erfüllt ist. Ist dies für ein ν = N−1 der Fall, so ist x(N ) = x1 der gesuchte Näherungs-
wert für x1 .
Falls der Radikand der Wurzel in (3.15) negativ ausfällt, so kann dies zwei Ursachen
haben:
(1) Eine reelle Lösung der Gleichung f (x) ≡ Pn (x) = 0 wird durch eine Folge konju-
giert komplexer Zahlen approximiert. Die Imaginärteile der Folge {x(ν) } sowie die
Imaginärteile der zugehörigen Polynomwerte streben dann gegen Null.
(2) x1 ist eine komplexe Nullstelle. Mit x1 ist dann auch x1 Nullstelle von Pn . In diesem
(N ) (N )
Fall liefert die Division Pn /[(x− x1 )(x− x1 )] unter Vernachlässigung des Restes
ein Polynom Pn−2 vom Grad n−2 (s. Abschnitt 3.2.2), für das wiederum die Muller-
Iteration eine Näherung für die im Allgemeinen betragskleinste Nullstelle liefert.
1. Schritt: Benutzung der Startwerte (3.18) x(0) = −1, x(1) = 1, x(2) = 0 und Be-
rechnung von fν := Pn (x(ν) ) für ν = 0, 1, 2 bzw. der künstlichen Funk-
tionswerte (3.19).
2. Schritt: Berechnung der Näherungen x(ν+1) , ν = 2, 3, 4 ... für x1 gemäß Itera-
tionsvorschrift (3.17) mit (3.16),(3.15),(3.14).
Nach jedem Iterationsschritt ist die Bedingung (3.21) zu prüfen. Ist sie
nicht erfüllt, so ist qν+1 zu halbieren, und es sind hν+1 , x(ν+1) und fν+1
neu zu berechnen und damit die Iteration fortzusetzen. Die Iteration ist
(N )
abzubrechen, falls (3.20) für ein ν = N−1 erfüllt ist. Dann ist x(N ) = x1
der gesuchte Näherungswert für x1 .
3. Schritt: (a) Berechnung des abdividierten Polynoms Pn−1 im Falle einer reel-
(N )
len Näherung x1 mit dem Horner-Schema (siehe Abschnitt 3.2.1)
(N )
unter Vernachlässigung des Restes, der entstehen kann, weil x1 ein
Näherungswert für x1 ist.
(b) Berechnung des abdividierten Polynoms Pn−2 im Falle einer kom-
(N )
plexen Lösung x1 mit dem doppelreihigen Horner-Schema (siehe
(N ) (N )
Abschnitt 3.2.2); es gilt Pn−2 (x) = Pn (x)/[(x − x1 )(x − x1 )],
(N ) (N )
wobei x1 die zu x1 konjugiert komplexe Lösung von Pn ist. Auch
hier ist der Rest zu vernachlässigen.
Alle weiteren Nullstellen werden nach dem gleichen Algorithmus berechnet, jeweils
angewandt auf die abdividierten Polynome f := Pn−1 , f := Pn−2 ,. . . . Man erhält so
sämtliche Lösungen in etwa dem Betrage nach geordnet.
Konvergenzordnung
Sei ξ eine einfache Nullstelle des Polynoms f (x) := Pn (x), so lässt sich zeigen, dass zwi-
schen den Fehlern
∆x(ν) := x(ν) − ξ
aufeinander folgender Muller-Iterationen die Beziehung
⎧ (3)
⎪
⎨ ∆x(ν+1) = ∆x(ν) · ∆x(ν−1) · ∆x(ν−2) · f (ξ) + O(ε)
6f (ξ) (3.22)
⎪
⎩ mit ξ = max |∆x |
(i)
i=ν−2,ν−1,ν
besteht. Heuristisch betrachtet ergibt sich aus (3.22) eine Beziehung der Form
Durch Vergleich zugehöriger Exponenten in (3.23) kann man wiederum folgern, dass die
Gleichungen
−2− 1 1+ 2 + 1
M = KM p p2 ⇒ K = M p p2
und
1 1
p=1+ + 2 ⇒ p3 − p2 − p − 1 = 0 (3.24)
p p
gelten müssen. Als einzige reelle Nullstelle der kubischen Gleichung (3.24) erhält man
nun für die Konvergenzordnung des Muller-Verfahrens im Falle einfacher Nullstellen den
Wert
p ≈ 1.84 .
Im Falle doppelter Nullstellen ermittelt man analog p ≈ 1.23.
Beispiel 3.8.
Gegeben: Das Polynom P3 mit
Gesucht: Eine (im Allgemeinen die betragskleinste) Nullstelle mit dem Muller-Verfah-
ren. Für die Abbruchbedingung wähle man die sehr grobe Schranke ε = 0.1,
damit nur wenige, mit dem Taschenrechner nachvollziehbare Schritte durch-
geführt werden müssen.
f0 := a0 − a1 + a2 = 281 ,
f1 := a0 + a1 + a2 = 209 ,
f2 := a0 = 252 .
|x(ν+1) − x(ν) |
= K (ν+1) < ε = 0.1
|x(ν+1) |
erhält man folgende Werte in der Tabelle bei 6-stelliger Mantisse:
108 3. Verfahren zur Lösung algebraischer Gleichungen
P2
(6)
1.00000 −1.04518 −42.2224 0.564652 = P3 (x1 )
Der Rest wird vernachlässigt, im Allgemeinen ist er vernachlässigbar klein, wenn die
Abbruchschranke ε entsprechend klein gesetzt wird, also mehr Iterationsschritte durch-
geführt werden. Das Deflationspolynom P2 lautet hier:
P2 (x) = x2 − 1.04518x − 42.2224.
Die Nullstellen von P2 würde man hier nicht mehr mit der Muller-Iteration, sondern mit
der Lösungsformel für die quadratische Gleichung bestimmen.
Da jedoch das Deflationspolynom P2 hier wegen der geringen Anzahl der Iterationen in
der Muller-Iteration zur Bestimmung von x1 so ungenau ist, wird hier auf die Lösung
der quadratischen Gleichung P2 (x) = 0 verzichtet. Stattdessen wird das (exakte) Defla-
tionspolynom P2 zu x1 = 6 bestimmt, es ergibt sich aus
P3 1 −7 −36 252
x1 = 6 0 6 −6 −252
P2 1 −1 −42 0
zu P2 (x) = x2 − x − 42.
&
7
Aus P2 (x) = 0 ⇒ x2,3 = 1
± 1
+ 42 = 1
± 13
= ⇒ x2 = −6, x3 = 7
2 4 2 2 −6
(exakte Lösungen).
3.3 Bestimmung von Lösungen algebraischer Gleichungen 109
Beispiel 3.9.
Gegeben: Tschebyscheff-Polynom 20. Grades. Alle ungeraden Koeffizienten (ai mit
ungeradem i) sind Null, die geraden haben folgende Werte:
i ai i ai i ai
0 1 8 549120 16 5570560
2 −200 10 −2050048 18 −2621440
4 6600 12 4659200 20 524288
6 −84480 14 −6553600
Lösung: Alle Nullstellen sind reell (d. h. die Imaginäranteile sind alle Null):
Als Iterationsverfahren wird das Verfahren von Newton verwendet. Die Iteration wird
abgebrochen, wenn z. B. die Abfrage |Pn (x(ν+1) )| < ε zu vorgegebenem ε > 0 erfüllt ist.
Gilt für ein festes ν
|Pn (x(ν) )| ≤ |Pn (x(ν+1) )|, (3.25)
so muss x(ν+1) aus der Folge der {x(ν) } ausgeschlossen werden. Mit einem zweidimen-
sionalen Suchprozess, der als Spiralisierung bezeichnet und komplex durchgeführt wird,
wird dann ein neues x(ν+1) ermittelt, für das
gilt; damit wird die Iteration fortgesetzt. Die Folgen der Näherungswerte werden durch
Extrapolation verbessert. Ist x(N ) der ermittelte Näherungswert, so wird er als Nullstelle
von Pn bezeichnet; man berechnet das Deflationspolynom Pn (x)/(x − x(N ) ) mit dem
Horner-Schema, vernachlässigt den Rest und wendet das eben beschriebene Verfahren
auf das Restpolynom Pn−1 vom Grad n−1 an. Analog fortfahrend erhält man alle Null-
stellen des Polynoms Pn .
Es sei x(ν+1) derjenige Wert der Folge {x(ν) }, für den erstmals (3.25) gilt. Dann muss
innerhalb eines Kreises um x(ν) mit dem Radius r = |x(ν+1) − x(ν) | ein xs+1 existieren
mit
|Pn (xs+1 )| < |Pn (x(ν) )|,
welches durch Absuchen des Kreisgebietes von außen nach innen mit einer Polygonspirale
ermittelt wird. Dazu wird mit einem komplexen Faktor q = q1 + iq2 , q1 , q2 reell, |q| < 1
gearbeitet, den Bauhuber q = 0.1 + 0.9i wählt (diese Wahl ist nicht bindend).
Im Falle ∆Pn ≤ 0 wird k um eins erhöht und (3.26) erneut berechnet. Analog wird so
lange fortgefahren, bis erstmals für ein k = s gilt ∆Pn > 0. Dann wird x(ν+1) ersetzt
durch xs+1 und damit nach Newton weiter iteriert.
3.3 Bestimmung von Lösungen algebraischer Gleichungen 111
Beispiel 3.10.
Das Verfahren von Jenkins und Traub ([JENK1970], [TRAU1966]) ist ein Iterations-
verfahren zur Ermittlung der betragskleinsten Nullstelle eines Polynoms Pn mit kom-
plexen Koeffizienten. Es ist für alle Startwerte x(0) ∈ [−∞, |xi |min ] global konvergent
von mindestens zweiter Ordnung. Es behandelt auch den Fall von zwei oder mehr be-
tragsgleichen Nullstellen. Je nachdem, ob die betragskleinste Nullstelle einfach, zweifach
oder mehr als zweifach ist, wird der vom Computer auszuführende Algorithmus automa-
tisch durch entsprechend eingebaute logische Entscheidungen modifiziert. Nachdem die
betragskleinste(n) Nullstelle(n) näherungsweise ermittelt ist (sind), wird durch Abdivi-
dieren der Nullstelle(n) das Restpolynom bestimmt. Hiervon liefert das gleiche Verfahren
eine Näherung für die nächste(n) Nullstelle(n).
112 3. Verfahren zur Lösung algebraischer Gleichungen
3.4 Anwendungsbeispiel
Beispiel 3.11.
Gegeben: Für die im Beispiel 43 (Turbinenbeispiel) hergeleitete kubische Gleichung er-
gibt sich mit den Daten für h, A, P in einem konkreten Fall
P3 (V̇ ) = V̇ 3 − 11144.16 · V̇ + 44233.6 = 0 .
Gesucht: Der optimale Volumenstrom V̇ , der mit Hilfe der Muller-Iteration berechnet
werden soll. Für die Abbruchbedingung wähle man die Schranke = 0.00001,
damit nur wenige, mit dem Taschenrechner nachvollziehbare Schritte, durch-
geführt werden müssen.
Lösung:
1. Schritt. (Bestimmung der Startwerte): Man setzt a0 = 44233.6, a1 = −11144.16 und
a2 = 0 und verwendet die künstlichen Funktionswerte (3.19), d. h.
f0 := a0 − a1 = 55377.76
f1 := a0 + a1 = 33089.44
f2 := a0 = 44233.6
Die kleinste Nullstelle und somit der optimale Volumenstrom ist erreicht bei
V̇ (5) = 3.974853244.
3.5 Entscheidungshilfen 113
3.5 Entscheidungshilfen
Die Entscheidungshilfen für die Auswahl eines geeigneten Verfahrens zur Nullstellenbe-
stimmung bei algebraischen Polynomen sind mit Vorbemerkungen und einem Überblick
in Abschnitt 3.3.1 zusammengefasst.
Direkte Verfahren
zur Lösung linearer Gleichungssysteme
Gegeben sei ein System von m linearen Gleichungen mit n Unbekannten xi der Form
⎧
⎪
⎪ a11 x1 + a12 x2 + . . . + a1n xn = a1 ,
⎪
⎨ a21 x1 + a22 x2 + . . . + a2n xn = a2 ,
.. .. (4.1)
⎪
⎪ . .
⎪
⎩
am1 x1 + am2 x2 + . . . + amn xn = am ,
wobei die Koeffizienten aik ∈ IR und die rechten Seiten ai ∈ IR, i = 1(1)m,
k = 1(1)n, vorgegebene Zahlen sind. In Matrizenschreibweise lautet (4.1)
mit ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
a11 a12 ··· a1n x1 a1
⎜ a21 a22 ··· a2n ⎟ ⎜ x2 ⎟ ⎜ a2 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
A = (aik ) = ⎜ .. .. ⎟, x =⎜ .. ⎟, a =⎜ .. ⎟.
⎝ . . ⎠ ⎝ . ⎠ ⎝ . ⎠
am1 am2 ··· amn xn am
116 4. Lösung linearer Gleichungssysteme
Ein Vektor x , dessen Komponenten xi , i = 1(1)n, jede der m Gleichungen des Systems
(4.1) zu einer Identität machen, heißt Lösungsvektor oder kurz Lösung von (4.1) bzw.
(4.2). In Abschnitt 4.13 werden überbestimmte Systeme Ax = a mit (m, n)-Matrix
A, a ∈ IRm , x ∈ IRn , m > n behandelt, sonst nur Systeme aus n Gleichungen mit n
Unbekannten, d. h. der Fall Ax = a mit (n, n)-Matrix A, x ∈ IRn , a ∈ IRn .
Theoretisch ist das Problem der Auflösung linearer Gleichungssysteme vollständig ge-
klärt. So lassen sich die Lösungen z. B. mit der Cramerschen Regel explizit berechnen;
dieser Weg erfordert jedoch einen immensen Rechenaufwand. Zur Lösung eines Systems
von n Gleichungen mit n Unbekannten erfordert die Cramersche Regel (n2 − 1) n! + n
Punktoperationen; dies sind 359 251 210 Punktoperationen für n = 10 und bereits 1021
für n = 20. Daher ist die Cramersche Regel (schon ab n = 3) für praktische Rechnungen
total unbrauchbar. Man muss also nach anderen, effektiveren Lösungsmöglichkeiten su-
chen. Ein solches Verfahren ist z. B. der Gaußsche Algorithmus, der für die Lösung von
n Gleichungen mit n Unbekannten nur n3 (n2 + 3n − 1) Punktoperationen benötigt, d. h.
im Fall n = 10 nur 430, im Fall n = 20 nur 3060 (gegenüber 1021 bei der Cramerschen
Regel).
Rechenzeit
(3 ns im Durchschnitt pro Multiplikation/Division)
Es gibt also brauchbare numerische Methoden. Trotz dieser Tatsache ist die praktische
Bestimmung der Lösung für sehr große n eine problematische numerische Aufgabe; die
Gründe hierfür sind:
1. Arbeitsaufwand (Rechenzeit) und Speicherplatzbedarf bei sehr großen Systemen,
siehe Beispiele
2. Verfälschung der Ergebnisse durch Rundungsfehler (numerische Instabilität)
3. Schlecht konditionierte Probleme (mathematische Instabilität)
Im Folgenden sollen nicht nur Lösungsmethoden angegeben, sondern auch auf numeri-
sche Effekte, anfallenden Arbeitsaufwand u.a.m. aufmerksam gemacht werden. Von einem
Lösungsverfahren wird gefordert, dass es Attribute wie Bandstruktur und Symmetrie
berücksichtigt, rechenzeitoptimal ist, Rundungsfehler unter Kontrolle gehalten werden
u.a.m.
Abb. 4.3.
3. Beispiel des Aerodynamischen Instituts der RWTH Aachen. Numerische Simulation
einer ablösenden Strömung um Tragflügelprofile, gerechnet mit den Navier-Stokes-
Gleichungen:
Wenn 3-dimensional gerechnet wird und ein (31 × 31 × 51)-Gitter mit je 4 Glei-
chungen verwendet wird, so erhält man nichtlineare Systeme aus 196 044 Glei-
chungen mit 196 044 Unbekannten, die iterativ (etwa mit 5 Iterationen) gelöst
werden. Rechnet man bis zum Wirbelablösen 10 000 Zeitschritte, so ergeben sich
5 × 10 000 = 50 000 lineare Gleichungssysteme aus rund ca. 200 000 Gleichungen,
die zu lösen sind.
Abb. 4.4.
4.1 Aufgabenstellung und Motivation 119
4. Beispiel des Instituts für Bildsame Formgebung der RWTH Aachen zum Freiform-
schmieden: Ein 8-kantiger Schmiedeblock aus Stahl wird mit 4 Hämmern bearbei-
tet, dann um 45◦ gedreht und wieder bearbeitet usw. Ziel ist die Berechnung des
Temperaturfeldes. Hier entstehen bei 1400 Knoten Gleichungssysteme von 4200
Gleichungen. Bei 50 Zeitschritten und 5 Iterationen pro Zeitschritt hat man also
250 lineare Gleichungssysteme mit je 4200 Gleichungen zu lösen.
Abb. 4.5.
Betrachtet man einen Rundblock mit 200 mm Durchmesser und 250 mm Länge, so
entstehen bei 3000 Knoten 9000 nichtlineare Gleichungen, d. h. bei 100 Zeitschrit-
ten mit 5 Iterationen pro Zeitschritt 500 lineare Systeme mit je 9000 Gleichungen.
Abb. 4.6.
120 4. Lösung linearer Gleichungssysteme
Beispiel 4.2.
⎛ ⎞
1 −2 3
A = ⎝ 4 5 6 ⎠ vom Typ (3, 3)
7 8 −9
Beispiel 4.4.
Die Matrix A aus dem Beispiel 4.2 ist regulär und zugleich streng regulär, während
⎛ ⎞
1 2 3
A = ⎝ 4 5 6 ⎠ mit det(A) = 0
7 8 9
Beispiel 4.6. ⎛ ⎞ ⎛ ⎞
1 0 0 5 7 −8
L1 = ⎝ 3 2 0 ⎠ R1 = ⎝ 0 1 3 ⎠
−3 5 7 0 0 4
untere Dreiecksmatrix obere Dreiecksmatrix
⎛ ⎞ ⎛ ⎞
1 0 0 1 7 −8
L2 = ⎝ 3 1 0 ⎠ R2 = ⎝ 0 1 3 ⎠
−3 5 1 0 0 1
normierte untere Dreiecksmatrix normierte obere Dreiecksmatrix
Beispiel 4.8.
⎛ ⎞ ⎛ ⎞
1 0 0 0 1 0
⎝ 0 aus E durch Vertauschung
E = 1 0 ⎠, P = ⎝ 1 0 0 ⎠
der 1. und 2. Zeile
0 0 1 0 0 1
⎛ ⎞⎛ ⎞ ⎛ ⎞
0 1 0 1 −2 3 4 5 6
PA = ⎝ 1 0 0 ⎠⎝ 4 5 6 ⎠=⎝ 1 −2 3 ⎠
0 0 1 7 8 −9 7 8 −9
122 4. Lösung linearer Gleichungssysteme
Satz 4.9.
Jede (n, n)-Matrix A mit det (Ak ) = 0 für k = 1(1)n−1 kann ohne Zeilenvertau-
schungen eindeutig in das Produkt LR aus einer normierten unteren Dreiecksmatrix
L und einer oberen Dreiecksmatrix R zerlegt werden:
⎛ ⎞ ⎛ ⎞
1 r11 r12 · · · r1n
⎜ l21 1 ⎟ ⎜ r22 · · · r2n ⎟
⎜ ⎟ ⎜ ⎟
L=⎜ . . ⎟ , R=⎜ .. .. ⎟.
⎝ .. .. ⎠ ⎝ . . ⎠
ln1 ln2 ··· 1 rnn
Beispiel 4.10.
⎛ ⎞
1 −1 −1
Gegeben: A = ⎝ −2 6 3 ⎠ vom Typ (3, 3).
−1 13 6
LR-Zerlegung von A:
⎛ ⎞ ⎛ ⎞⎛ ⎞
1 −1 −1 1 0 0 1 −1 −1
A = ⎝ −2 6 3 ⎠ = ⎝ −2 1 0 ⎠⎝ 0 4 1 ⎠ = LR
−1 13 6 −1 3 1 0 0 2
Die Multiplikation der beiden Matrizen L und R rechts und ein Koeffizientenvergleich
mit der Matrix A links ergibt zeilenweise die zu berechnenden Koeffizienten der Drei-
ecksmatrizen L und R.
!
1. Zeile: a11 = 1 = 1 · r11 ⇒ r11 = 1
a12 = −1 = 1 · r12 ⇒ r12 = −1
a13 = −1 = 1 · r13 ⇒ r13 = −1
2. Zeile: a21 = −2 = l21 · r11 + 1 · 0 ⇒ l21 = −2
a22 = 6 = l21 · r12 + 1 · r22 ⇒ r22 = 4
a23 = 3 = l21 · r13 + 1 · r23 + 0 · r33 ⇒ r23 = 1
3. Zeile: a31 = −1 = l31 · r11 ⇒ l31 = −1
a32 = 13 = l31 · r12 + l32 · r22 ⇒ l32 = 4 (13 − 1) = 3
1
4.2 Definitionen und Sätze 123
Algorithmus 4.11.
Die Elemente der Matrizen L und R gemäß Dreieckszerlegung nach Satz 4.9 werden
wie folgt berechnet:
1. Für k = 1(1)n: r1k := a1k
2. Für jedes i = 2(1)n
ai1
2.1 Für k = 1: li1 :=
a11
k−1
2.2 Für k = 2(1)i−1: lik := (aik − lij rjk )/rkk
j=1
i−1
2.3 Für k = i(1)n: rik := aik − lij rjk
j=1
Satz 4.12.
Für eine (n, n)-Matrix A mit det A = 0 gilt mit einer (n, n)-Permutationsmatrix P
die Zerlegung
PA = LR,
wobei L und R durch P und A eindeutig bestimmt sind.
In PA sind die Zeilen von A permutiert. Es gilt mit det(P) = (−1)k , k = Anzahl der
Zeilenvertauschungen,
Jede reguläre Matrix A kann somit durch Linksmultiplikation mit einer Permutationsma-
trix P (also durch Zeilenvertauschungen) in eine streng reguläre Matrix PA transformiert
werden.
Beispiel 4.13.
Für die Matrix
⎛ ⎞
1 −1 −1
A = ⎝ −1 1 −1 ⎠ , det(A2 ) = 0
−1 −1 1
existiert keine LR-Zerlegung ohne Zeilenvertauschung nach Satz 4.9. Nach Vertauschung
124 4. Lösung linearer Gleichungssysteme
der 2. und der 3. Zeile von A lässt sich eine eindeutige LR-Zerlegung angeben, denn mit
⎛ ⎞⎛ ⎞ ⎛ ⎞
1 0 0 1 −1 −1 1 −1 −1
PA = ⎝ 0 0 1 ⎠ ⎝ −1 1 −1 ⎠ = ⎝ −1 −1 1 ⎠ ,
0 1 0 −1 −1 1 −1 1 −1
gilt ⎛ ⎞⎛ ⎞
1 0 0 1 −1 −1
PA = LR = ⎝ −1 1 0 ⎠⎝ 0 −2 0 ⎠
−1 0 1 0 0 −2
Beispiel 4.15.
⎛ ⎞
1 2
1 3 4
A = ⎝ 3 4 ⎠ vom Typ (3, 2) ⇒ T
A = vom Typ (2, 3)
2 4 6
4 6
Beispiel 4.17.
⎛ ⎞
1 2 −3
A=⎝ 2 4 5 ⎠ = AT ⇒ A ist symmetrisch
−3 5 6
Beispiel 4.19.
⎛ ⎞ ⎛ ⎞
cos α 0 sin α cos α 0 − sin α
Q1 = ⎝ 0 1 0 ⎠ , Q1T = ⎝ 0 1 0 ⎠
− sin α 0 cos α sin α 0 cos α
1 2 −1 1 2 1
Q2 = √ , =√ Q2T
,
5 1 2 5 −1 2
1 2 −1 2 1 1 5 0 1 0
⇒ Q2 · Q2T = = = =E
5 1 2 −1 2 5 0 5 0 1
Satz 4.20.
Sei v ∈ IRn ein Vektor und E die (n, n)-Einheitsmatrix. Dann ist
2
H := E − vv T
v 2
Beispiel 4.21.
2 2 T 2
v = , v = v v = (2, 1) = 5
1 1
2 4 2
v vT = (2, 1) =
1 2 1
8 4
1 0 4 2 1 0
E − 2 2 v vT =
5 5
H := − 25 = −
v 0 1 2 1 0 1 4 2
5 5
1− 5
8
−5 4
−5 −5
3 4
= =
− 54 1 − 25 − 45 3
5
− 35 − 54 − 53 − 45 1 0
H TH = = ,
− 45 3
5 − 54 3
5 0 1
Diagonalmatrizen mit m = mr = 0,
bidiagonale Matrizen mit m = 1, mr = 0 oder m = 0, mr = 1,
tridiagonale Matrizen mit m = mr = 1,
fünfdiagonale Matrizen mit m = mr = 2.
Beispiel 4.23.
⎛ ⎞ ⎛ ⎞
1 0 0 0 1 0 0 0
⎜ 0 2 0 0 ⎟ ⎜ 2 1 0 0 ⎟
A = ⎜
⎝
⎟ B = ⎜ ⎟
0 0 5 0 ⎠ ⎝ 0 −1 2 0 ⎠
0 0 0 −6 0 0 3 4
Diagonalmatrix untere Bidiagonalmatrix
⎛ ⎞ ⎛ ⎞
1 2 0 0 1 2 0 0
⎜ 0 1 3 0 ⎟ ⎜ −1 2 3 0 ⎟
C = ⎜
⎝
⎟ D = ⎜ ⎟
0 0 4 5 ⎠ ⎝ 0 4 5 6 ⎠
0 0 0 1 0 0 1 3
obere Bidiagonalmatrix tridiagonale Matrix
⎛ ⎞
1 2 1 0 0
⎜ 3 4 5 2 0 ⎟
⎜ ⎟
F = ⎜ ⎜ 3 1 −2 4 3 ⎟
⎟ fünfdiagonale Matrix
⎝ 0 5 3 4 5 ⎠
0 0 7 1 2
Dies bedeutet, dass bei einer zyklisch tridiagonalen Matrix die Hauptdiagonale und die
beiden Nebendiagonalen mit Nichtnull-Elementen besetzt sein können und mindestens ei-
nes der beiden Elemente links unten (an1 ) und rechts oben (a1n ) von Null verschieden ist.
4.2 Definitionen und Sätze 127
Beispiel 4.25.
⎛ ⎞
⎛ ⎞ 1 2 0 0 4
1 2 0 1 ⎜ ⎟
⎜ −3 4 −3 4 1 0 0
1 0 ⎟ ⎜ ⎟
A=⎜
⎝ 0 1
⎟, B =⎜
⎜ 0 −1 5 −1 0 ⎟
⎟
2 1 ⎠ ⎝ ⎠
0 0 −3 6 3
−1 0 2 3
−5 0 0 2 1
sind zyklisch tridiagonale Matrizen.
für mindestens ein i muss das Zeichen echt größer“ gelten (d. h. der Betrag des Dia-
”
gonalelements ist größer oder gleich der Summe der Beträge der anderen Elemente in
derselben Zeile; für mindestens eine Zeile muß dabei echt größer“ gelten).
”
Die Matrix heißt stark diagonaldominant, wenn gilt
n
|aii | > |aik | für i = 1(1)n .
k=1
k=i
Das heißt, dass für alle Zeilen der Betrag des Diagonalelements echt größer als die
Summe der Beträge der übrigen Elemente der Zeile ist.
Beispiel 4.27.
⎛ ⎞
10 −1 2
A = ⎝ 5 8 1 ⎠ ist diagonaldominant
−1 2 3m
⎛ ⎞
10 −1 2
A = ⎝ 5 8 1 ⎠ ist stark diagonaldominant
−1 2 5m
Beispiel 4.29. ⎛ ⎞
5 −1 −2
A = ⎝ −1 6 −1 ⎠ = AT symmetrisch
−2 −1 5
ist positiv definit, da sämtliche Hauptabschnittsdeterminanten positiv sind
2. Die Zerlegung A = LR mit A = AT gemäß Satz 4.9 führt auf ein R mit rii > 0
für alle i.
Beispiel 4.30.
⎛ ⎞ ⎛ ⎞⎛ ⎞
4 −2 0 1 0 0 4m−2 0
⎜ ⎟ ⎜ 1 ⎟⎜ ⎟
A = ⎝ −2 5 −2 ⎠ = ⎝ − 2 1 0 ⎠⎝ 0 4m−2 ⎠ = LR
0 −2 5 0 − 12 1 0 0 4m
Beispiel 4.31.
⎛ ⎞ ⎛ ⎞
16 0 4 3 −1 −2
zu 1. A = ⎝ 0 9 0 ⎠, zu 2. B = ⎝ −1 6 −3 ⎠
4 0 17 −2 −3 8
⎛ ⎞
10 1 0 0
⎜ 1 10 −2 0 ⎟
zu 3. C = ⎜ ⎟
⎝ 0 −2 5 −3 ⎠
0 0 −3 6
Satz 4.32.
Eine stark diagonaldominante Matrix A = (aik ), i, k = 1(1)n, ist streng regulär.
Beispiel 4.33.
⎛ ⎞
10 −1 2
A=⎝ 5 8 1 ⎠ streng regulär, da det(A1 ) = |10| = 10 = 0 ,
−1 2 5
10 −1
det(A2 ) = = 85 = 0 , det(A3 ) = det(A) = 442 = 0
5 8
Satz 4.34.
Jede symmetrische, streng reguläre Matrix A = (aik ), i, k = 1(1)n, kann eindeutig in
das Produkt RT DR mit einer normierten oberen Dreiecksmatrix R, ihrer Transpo-
nierten RT und einer Diagonalmatrix D zerlegt werden.
⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
a11 a12 a13 1 0 0 d11 0 0 1 r12 r13
⎝ a12 a22 a23 ⎠ = ⎝ r12 1 0 ⎠ ⎝ 0 d22 0 ⎠ ⎝ 0 1 r23 ⎠
a13 a23 a33 r13 r23 1 0 0 d33 0 0 1
⎛ ⎞ ⎛ ⎞
1 0 0 d11 d11 r12 d11 r13
= ⎝ r12 1 0 ⎠ ⎝ 0 d22 d22 r23 ⎠
r13 r23 1 0 0 d33
⎛ ⎞
d11 d11 r12 d11 r13
= ⎝ r12 d11 2
d11 r12 + d22 d11 r13 r12 + d22 r23 ⎠ .
2 2
r13 d11 d11 r12 r13 + d22 r23 d11 r13 + d22 r23 + d33
130 4. Lösung linearer Gleichungssysteme
a22 = 2
d11 r12 + d22 ⇒ d22 = a22 − d11 r12
2
a23 = d11 r13 r12 + d22 r23 ⇒ r23 = (a23 − d11 r13 r12 ) / d22
a33 = 2
d11 r13 2
+ d22 r23 + d33 ⇒ d33 = a33 − d11 r13
2
− d22 r23
2
Algorithmus 4.35.
Gegeben: A = (aik ), i, k = 1(1)n, symmetrisch, streng regulär
Gesucht: Elemente von D und R bei der Zerlegung A = RT DR
1. d11 = a11
a1k
r1k = d11 , k = 2(1)n
n−1
3. dnn = ann − 2
djj rjn
j=1
⎛ A ⎞ = ⎛ RT ⎞ ⎛ D ⎞ ⎛ R ⎞
4 −2 0 1 0 0 4 0 0 1 − 12 0
⎝ −2 5 −2 ⎠ = ⎝ −1 1 0 ⎠ ⎝ 0 4 0 ⎠ ⎝ 0 1 − 21 ⎠
2
0 −2 5 0 − 12 1 0 0 4 0 0 1
Satz 4.37.
Jede symmetrische, positiv definite Matrix A = (aik ), i, k = 1(1)n, ist in das Produkt
RT DR gemäß Satz 4.34 zerlegbar, wobei alle Diagonalelemente von D (Pivotelemen-
te) positiv sind und A streng regulär ist.
4.2 Definitionen und Sätze 131
Satz 4.38.
Jede symmetrische, positiv definite Matrix A = (aik ), i, k = 1(1)n, kann eindeutig in
das Produkt RT R mit der oberen Dreiecksmatrix R = (rik ), rii > 0 für alle i, und
ihrer Transponierten RT zerlegt werden. Die Zerlegung heißt Cholesky-Zerlegung.
Cholesky-Zerlegung für n = 3: A = RT R
⎛ ⎞ ⎛ ⎞⎛ ⎞
a11 a12 a13 r11 0 0 r11 r12 r13
⎝ a12 a22 a23 ⎠ = ⎝ r12 r22 0 ⎠ ⎝ 0 r22 r23 ⎠
a13 a23 a33 r13 r23 r33 0 0 r33
⎛ 2 ⎞
r11 r11 r12 r11 r13
= ⎝ r11 r12 r12
2 2
+ r22 r12 r13 + r22 r23 ⎠
2 2 2
r11 r13 r12 r13 + r22 r23 r13 + r23 + r33
a23 = r12 r13 + r22 r23 ⇒ r23 = (a23 − r12 r13 )/r22
"
a33 = 2
r13 2
+ r23 2
+ r33 ⇒ r33 = a33 − r13 2 − r2
23
Algorithmus 4.39.
Beispiel 4.40.
⎛ ⎞
4 −2 0
Gegeben: Die Matrix A = ⎝ −2 5 −2 ⎠.
0 −2 5
A = RT R
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
4 −2 0 2 0 0 2 −1 0
⎝ −2 5 −2 ⎠ = ⎝ −1 2 0 ⎠ ⎝ 0 2 −1 ⎠
0 −2 5 0 −1 2 0 0 2
4.3 Lösbarkeitsbedingungen
für ein lineares Gleichungssystem
Ein lineares Gleichungssystem Ax = a mit einer (m, n)-Matrix A, x ∈ IRm , Rg(A) = r,
m ≥ n, heißt homogen, wenn a = 0 ist, andernfalls heißt es inhomogen.
Satz 4.41.
Ein inhomogenes Gleichungssystem Ax = a = 0 ist genau dann auflösbar, wenn der
Rang der erweiterten Matrix (A, a) gleich dem Rang der Matrix A ist: Rg(A, a) =
Rg(A). Die Gesamtheit der Lösungen setzt sich aus der Lösung xh des homogenen
Systems und einer speziellen Lösung des inhomogenen Systems zusammen.
4.4 Prinzip der direkten Methoden zur Lösung linearer Gleichungssysteme 133
Für Systeme (4.2) mit m = n gilt in der Formulierung über die Determinante
a) det A = 0: Es existiert genau eine Lösung, sie lautet x = A−1 a.
(Hier gilt: Rg(A, a) = Rg(A) = n)
b) det A = 0: Ist das System auflösbar, d. h. Rg(A, a) = Rg(A) < n, so ist die Lösung
nicht eindeutig bestimmt. Sie ergibt sich als Summe aus einer Linearkombination
der n−r linear unabhängigen Lösungen des homogenen Systems und einer speziellen
Lösung des inhomogenen Systems.
Aus den Sätzen 4.9 bzw. 4.12 ergeben sich die folgenden Algorithmen für die Lösung
linearer Systeme durch Dreieckszerlegung (Faktorisierung).
Beweis. Aus A x = a folgt nach Satz 4.9 mit A = LR die Beziehung LR x = a. Ziel
ist die Herstellung der oberen Dreiecksform R x = r , so dass gilt
Rx =a
L =⇒ Lr =a .
:=r
Das heißt man muss nach der Faktorisierung A = LR zuerst L r = a durch Vorwärtse-
limination lösen, bevor man aus R x = r rekursiv die Lösung x ermitteln kann. Daraus
ergibt sich die Reihenfolge der Schritte im Algorithmus.
Beispiel 4.43.
Gegeben: A x = a mit
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
4 −1 −1 2 x1
A=⎝ 8 0 −1 ⎠ a =⎝ 7 ⎠ x = ⎝ x2 ⎠
4 1 4 9 x3
Lösung:
1. Schritt: Dreieckszerlegung A = LR =⇒ L, R
⎛ ⎞ ⎛ ⎞⎛ ⎞
4 −1 −1 1 0 0 4 −1 −1
A=⎝ 8 0 −1 ⎠ = ⎝ 2 1 0 ⎠⎝ 0 2 1 ⎠
4 1 4 1 1 1 0 0 4
2. Schritt: Vorwärtselimination L r = a
⎛ ⎞⎛ ⎞ ⎛ ⎞
1 0 0 r1 2
⎝ 2 1 0 ⎠ ⎝ r2 ⎠ = ⎝ 7 ⎠
1 1 1 r3 9
⎫ ⎛ ⎞
r1 = 2 ⎬ 2
2r1 + r2 = 7 ⇒ r2 = 7−4 = 3 =⇒ r =⎝3⎠.
⎭
r1 + r2 + r3 = 9 ⇒ r3 = 9−3−2 = 4 4
3. Schritt: Rückwärtselimination
⎛ ⎞⎛ ⎞ ⎛ ⎞
−1 −1
4 x1 2
⎝ 02 1 ⎠ ⎝ x2 ⎠ = ⎝ 3 ⎠
00 4 x3 4
⎫ ⎛ ⎞
4x3 = 4 ⇒ x3 = 1 ⎬ 1
2x2 + x3 = 3 ⇒ x2 = (3 − 1)/2 = 1 =⇒ x =⎝1⎠
⎭
4x1 − x2 − x3 = 2 ⇒ x1 = (2 + 1 + 1)/4 = 1 1
4.4 Prinzip der direkten Methoden zur Lösung linearer Gleichungssysteme 135
Aus Satz 4.12 ergibt sich für A x = a durch Linksmultiplikation mit der Permutations-
matrix P
PA x = P a
und mit PA = LR
Rx = P a ,
L
Beispiel 4.45.
Gegeben: A x = a mit ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
3 3 6 −3 x1
A=⎝ 2 2 3 ⎠ a = ⎝ −1 ⎠ x = ⎝ x2 ⎠
1 0 1 0 x3
Lösung:
1. Schritt: PA = LR
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎛ ⎞
1 0 0 3 3 6 3 3 6 1 0 0 3 3 6
⎜ ⎟⎜ ⎟
PA = ⎝ 0 0 1 ⎠ ⎝ 2 2 3 ⎠ = ⎝ 1 0 1 ⎠ = LR = ⎝ 13 1 0 ⎠ ⎝ 0 −1 −1 ⎠
0 1 0 1 0 1 2 2 3 2
0 −1
3 0 1 0
2. Schritt: P a = L r
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎛ ⎞ ⎛ ⎞
1 0 0 −3 −3 1 0 0 r1 −3
⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 1 ⎟⎜ ⎟ ⎜ ⎟
⎝ 0 0 1 ⎠ ⎝ −1 ⎠ = ⎝ 0 ⎠ = ⎝ 3 1 0 ⎠ ⎝ r2 ⎠ =⇒ r =⎝ 1 ⎠
0 1 0 0 −1 2
3 0 1 r3 1
136 4. Lösung linearer Gleichungssysteme
Aus PA = LR folgt
Die folgenden direkten Eliminationsverfahren arbeiten nach den angegebenen Algorith-
men. Sie unterscheiden sich lediglich dadurch, dass spezielle Eigenschaften der Matrix
A in Ax = a ausgenutzt werden, wodurch eine zum Teil erhebliche Ersparnis an
Rechenaufwand erreicht werden kann. Im Wesentlichen wird hier nur mit Systemen aus n
Gleichungen für n Unbekannte gearbeitet, lediglich in Abschnitt 4.13 mit überbestimm-
ten Systemen. Bevor aber spezielle Lösungsverfahren für lineare Systeme mit speziellen
Matrizen entwickelt werden, soll der Gauß-Algorithmus in der Form des klassischen Gauß-
Schemas hergeleitet werden, und danach wird der Zusammenhang zur Dreieckszerlegung
hergestellt.
Das Prinzip des Gaußschen Algorithmus ist die Überführung eines Gleichungssystems
der Form (4.1) mit m = n in ein gestaffeltes System in oberer Dreiecksform
⎧
⎪
⎪ r11 x1 + r12 x2 + · · · + r1n xn = r1 ,
⎪
⎨ r22 x2 + · · · + r2n xn = r2 ,
.. .. .. (4.3)
⎪
⎪ . . .
⎪
⎩
rnn xn = rn ,
aus dem sich die xi , i = 1(1)n, durch Rückwärtselimination berechnen lassen, falls
r11 r22 . . . rnn = 0 ist.
Zunächst wird der Algorithmus hergeleitet und anschließend wird der Zusammenhang
mit den Algorithmen in Abschnitt 4.4 hergestellt.
dass das betragsgrößte Element der ersten Spalte von A in die erste Zeile kommt (Spal-
tenpivotsuche vgl. Abschnitt 4.5.2). Die durch die Umordnung entstandene Matrix heiße
(0) (0)
A(0) , ihre Elemente aik und die Komponenten der rechten Seite ai , so dass (4.1) in
das äquivalente System
n
(0) (0)
aik xk = ai , i = 1(1)n, (4.4)
k=1
übergeht. Ist det A = 0, so gilt für das betragsgrößte Element (Pivotelement) der er-
(0)
sten Spalte a11 = 0. Zur Elimination von x1 aus den Gleichungen i = 2(1)n multipli-
(0) (0)
ziert man die 1. Gleichung von (4.4) mit −ai1 /a11 und addiert sie jeweils zur i-ten
Gleichung, so dass sich für i = 2(1)n zusammen mit der unveränderten 1. Zeile ergibt
(1. Eliminationsschritt):
.. .. (4.5)
⎪
⎪
⎪
⎪ . .
⎩ (1) (1) (1)
an2 x2 ann xn =
+ . . . + an ,
mit
⎧
⎪
⎪ 0 für k = 1, i = 2(1)n,
⎨ (0)
(1) (1) (0) (0) ai1
aik = (0)
ai = ai − a1 , i = 2(1)n.
⎪
⎪ (0) (0) a (0)
⎩ aik − a1k i1 (0) sonst, a11
a11
Das System (4.5) besteht also aus einer Gleichung mit den n Unbekannten x1 , x2 , . . . , xn
und n−1 Gleichungen mit den n−1 Unbekannten x2 , . . . , xn .
Auf die n−1 Gleichungen i = 2(1)n von (4.5) wendet man das Eliminationsverfahren er-
neut an. Dazu muss man zunächst wieder eine Zeilenvertauschung durchführen, so dass
(1)
das betragsgrößte Element der
ai2 für i = 2(1)n in der 2. Gleichung erscheint; nach der
(1)
Zeilenvertauschung werden die Elemente der neu entstandenen Zeilen 2 bis n mit aik
(1)
bzw. ai bezeichnet:
⎧
⎪
⎪
(0)
a11 x1 +
(0)
a12 x2
(0)
+ . . . +a1n xn =
(0)
a1 ,
⎪
⎪
⎪
⎨ (1) (1) (1)
a22 x2 + . . . +a2n xn = a2 ,
(4.6)
⎪
⎪ .. ..
⎪
⎪ . .
⎪
⎩ (1) (1) (1)
an2 x2 + . . . +ann xn = an ,
(1)
wobei wegen det A = 0 gelten muss a22 = 0.
Verfährt man nun analog mit der 2. bis n-ten Gleichung von (4.6), so sind für jeden
138 4. Lösung linearer Gleichungssysteme
(j−1)
(j) (j−1) (j−1) aij
ai = ai − aj (j−1) , i = (j + 1)(1)n,
ajj
(n−1)
an aj
(j−1)
n (j−1)
ajk
xn = (n−1)
, xj = (j−1)
− (j−1)
xk , j = n−1, n−2, . . . , 1 . (4.8)
ann ajj k=j+1 ajj
(j−1)
Im Fall det A = 0 darf keines der Diagonalelemente ajj verschwinden. Ist es nach ir-
(j−1)
gendeinem Eliminationsschritt nicht mehr möglich, ein Element ajj = 0 zu finden, so
bedeutet dies, dass det A = 0 ist. Ob dann überhaupt eine Lösung existiert und wenn ja,
wieviele Parameter sie besitzt, folgt automatisch aus der Rechnung. Für die Determinan-
te von A gilt det A = (−1)k r11 r22 . . . rnn , wobei k die Anzahl der Zeilenvertauschungen
ist.
Da der Rang r von A gleich der Anzahl der nicht verschwindenden Diagonalelemente
(j−1)
rjj = ajj der Superdiagonalmatrix R (gegebenenfalls unter Spaltenvertauschungen)
ist, lässt sich die Anzahl n−r der Parameter nach Durchführung der n−1 Eliminations-
schritte sofort angeben.
4.5 Der Gauß-Algorithmus 139
Die Zeilen 1(0) , 2(1) , 3(2) bilden das gesuchte gestaffelte System (4.7), aus dem die Lösun-
gen xi rekursiv gemäß (4.8) bestimmt werden. Die Zeilenvertauschung der Zeilen 2(1) ,
3(1) erübrigt sich, falls |
(1)
a22 | ≥ |
(1)
a32 | ist; dann ist
(1) (1)
a2i = a2i und
(1) (1)
a3i = a3i für i = 2, 3
zu setzen.
Beispiel 4.47.
Gegeben: Das Gleichungssystem Ax = a mit
⎛ ⎞ ⎛ ⎞ ⎛
⎞
3 3 6 x1 −3
A=⎝ 2 2 3 ⎠, x = ⎝ x2 ⎠ , a =⎝ 1 ⎠
1 0 1 x3 0
Lösung:
Bezeichnung A a erfolgte Operationen
der Zeilen
1(0) 3m 3 6 −3
(0)
2 2 2 3 1 —
(0)
3 1 0 1 0
2
2(1) 0 0 −1 3 − · 1(0) + 2(0)
3m
1
3(1) 0 −1 −1 1 − · 1(0) + 3(0)
3m
Bemerkung: Würde man hier ohne Zeilenvertauschung von 2(1) und 3(1) arbeiten, so
käme es beim 2. Eliminationsschritt zur Division durch Null; dies wird durch die Zeilen-
vertauschung vermieden.
Beispiel 4.48.
Gegeben: Ein homogenes System Ax = 0 .
Gesucht: Die Lösungen des homogenen Systems mit dem Gaußschen Algorithmus.
Lösung:
1(0) 2 −2 4 0
2(0) −1 2 3 0 —
3(0) 1 −1 2 0
2(1) 0 1 5 0 − − 21 · 1(0) + 2(0)
3(1) 0 0 0 0 − 12 · 1(0) + 3(0)
Wegen det A = 0 existiert eine nichttriviale Lösung. Aus 3(1) folgt: x3 = t, t ∈ IR, aus
2(1) : x2 = −5t und aus 1(0) : x1 = −7t, so dass sich der vom Parameter t abhängige
Lösungsvektor ergibt ⎛ ⎞ ⎛ ⎞
x1 −7
⎝ x2 ⎠ = t ⎝ −5 ⎠
x3 1
d. h. alle Punkte einer Gerade durch den Nullpunkt sind Lösung des homogenen Sy-
stems.
4.5.2 Spaltenpivotsuche
Wenn die Koeffizienten gerundete Zahlen sind oder im Verlaufe der Rechnung gerundet
werden muss, sind Zeilenvertauschungen unerlässlich, um Verfälschungen des Ergebnis-
ses durch Rundungsfehler möglichst zu dämpfen. Man bezeichnet diese Strategie als
(j−1)
Spaltenpivotsuche oder teilweise Pivotsuche und die Diagonalelemente rjj = ajj als
Pivotelemente. Unter Verwendung der Spaltenpivotsuche wird der Gauß-Algorithmus in
den meisten Fällen stabil. Vollkommen stabil wird er, wenn man als Pivotelement je-
weils das betragsgrößte Element der gesamten Restmatrix verwendet, man spricht dann
von vollständiger Pivotsuche. Hierfür ist der Aufwand sehr groß; für die Praxis ist die
Spaltenpivotsuche bzw. die weiter unten erläuterte skalierte Pivotsuche in vielen Anwen-
dungsfällen ausreichend.
142 4. Lösung linearer Gleichungssysteme
Beispiel 4.49.
Gegeben: Das lineare Gleichungssystem
(1) 0.2420 · 10−3 x1 + 0.6004 · 10−2 x2 = 0.1743 · 10−2 ,
(2) 0.4000 · 100 x1 + 0.9824 · 101 x2 = 0.2856 · 101 .
Gesucht: Zu bestimmen sind x1 und x2 nach dem Gaußschen Algorithmus; die exakten
Lösungen sind x1 = 1, x2 = 0.25.
Lösung:
1. Fall: Die Gleichungen werden in der obigen Reihenfolge benutzt. Zur Elimination
von x1 aus Gleichung (2) wird (1) mit
0.4000
− = −0.1653 · 104 = c1
0.2420 · 10−3
multipliziert und c1 · (1) + (2) gebildet. Man erhält das gestaffelte System
0.2420 · 10−3 x1 + 0.6004 · 10−2 x2 = 0.1743 · 10−2
−0.1010 · 10 x2 = −0.2500 · 10−1 .
0
Als Lösungen ergeben sich daraus: x1 = 0.1062 · 101; x2 = 0.2475, sie weichen
stark von den exakten Lösungen ab.
2. Fall: Die Reihenfolge der Gleichungen (1) und (2) wird vertauscht, so dass jetzt
das betragsgrößte Element der ersten Spalte in der ersten Gleichung steht
1(0) 0.4000 · 100 x1 + 0.9824 · 101 x2 = 0.2856 · 101 ,
2(0) 0.2420 · 10−3 x1 + 0.6004 · 10−2 x2 = 0.1743 · 10−2 .
0.2420 · 10−3
− = −0.6050 · 10−3 = c2
0.4000
multipliziert und c2 · 1(0) + 2(0) gebildet. Man erhält das gestaffelte System
0.4000 · 100 x1 + 0.9824 · 101 x2 = 0.2856 · 101 ,
0.6000 · 10−4 x2 = 0.15000 · 10−4 .
Hieraus ergeben sich als Lösungen x1 = 0.1000 · 101 und x2 = 0.2500, sie
stimmen mit den exakten Lösungen überein.
Die Abweichungen der Lösungen im 1. Fall sind dadurch entstanden, dass mit dem sehr
großen Faktor c1 multipliziert wurde und die Rundungsfehler dadurch entsprechend an-
gewachsen sind.
4.5 Der Gauß-Algorithmus 143
Beispiel 4.50.
Gegeben: Das lineare Gleichungssystem
⎛ ⎞⎛ ⎞ ⎛ ⎞
7 8 9 x1 24
⎝ 8 9 10 ⎠ ⎝ x2 ⎠ = ⎝ 27 ⎠ .
9 10 8 x3 27
Gesucht: Eine näherungsweise Lösung mit dem Gaußschen Algorithmus unter Verwen-
dung der Gleitpunktarithmetik mit dreistelliger Mantisse
(i) ohne Pivotisierung
(ii) mit Spaltenpivotsuche.
Die exakte Lösung ist xex = (1, 1, 1)T .
Lösung:
−0.300
3(2) 0.00 0.00 −2.85 −3.00 − −0.120 · 2(1) + 3(1)
Bezeichnung A a Operationen
der Zeilen
1(0) 7.00 8.00 9.00 24.0
2(0) 8.00 9.00 10.0 27.0 —
3(0) 9.00 10.0 8.00 27.0
− 8.00
2(1) 9.00 · 1
(0)
0.00 0.110 2.89 3.00 + 2(0)
− 7.00
3(1) 9.00 · 1
(0)
0.00 0.220 2.78 3.00 + 3(0)
Bemerkung. (Pivotsuche ist nur entscheidend wirkungsvoll, wenn alle Zeilen- und Spal-
tenbetragssummen annähernd gleich groß sind)
n
n
zi := |aij | ≈ sk := |ajk | für i, k = 1(1)n. (4.9)
j=1 j=1
Matrizen, für die (4.9) gilt, heißen äquilibriert (vgl. [MAES1985], 2.2.2). Da sich aber
der Rechenaufwand bei einer Äquilibrierung (vgl. Abschnitt 4.14.3) beträchtlich erhöhen
würde, führt man bei nicht äquilibrierten Matrizen statt der Spaltenpivotsuche eine so-
genannte skalierte Spaltenpivotsuche durch. Man ersetzt bei der Elimination das Glied
(j−1) (j−1)
|ajj | = max | aij |
i=j(1)n
4.5 Der Gauß-Algorithmus 145
xn = rrnn
n ,
⎛ ⎞
n
xi = r1ii ⎝ri − rij xj ⎠ , i = n−1, n−2, . . . , 1.
j=i+1
Während des Eliminationsprozesses können die Elemente von A zeilenweise mit den
Elementen der Zerlegungsmatrizen überspeichert werden; die lij für i > j stehen dann
unterhalb der Hauptdiagonale, die lii = 1 werden nicht abgespeichert und die rij (i ≤ j)
stehen in und über der Hauptdiagonale. Ebenso können die Elemente von a durch die
von r überspeichert werden.
Beispiel 4.51.
2. Vorwärtselimination Lr = a:
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
1 0 0 r1 4 r1 4 ⏐
⏐
⎝ 1 1 0 ⎠ ⎝ r2 ⎠ = ⎝ 5 ⎠ ⇒ r = ⎝ r2 ⎠ = ⎝ 1 ⎠ ⏐
-
2 −3 1 r3 7 r3 2
3. Rückwärtselimination Rx = r:
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
2 1 1 x1 4 x1 1 .
⏐
⎝ 0 2 −1 ⎠ ⎝ x2 ⎠ = ⎝ 1 ⎠ ⇒ x = ⎝ x2 ⎠ = ⎝ 1 ⎠ ⏐
⏐
0 0 2 x3 2 x3 1
⇒ 1(0) 2 1 1 4
2(0) 2 3 0 5 — —
3(0) 4 −4 7 7
2 2
⇒ 2(1) =
2(1) 0 2m −1 1 − · 1(0) + 2(0) l21 = 2 =1
2
4 4
3(1) =
3(1) 0 −6 5 −1 − · 1(0) + 3(0) l31 = 2 =2
2
−6 −6
⇒ 3 (2)
0 0 2 2 − · 2(1) + 3(1) l32 = = −3
2m
2
Aus den Zeilen 1(0) , 2(1) und 3(2) ergeben sich die obere Dreiecksmatrix R sowie die
rechte Seite r des Systems R x = r .
Die Faktoren der Operationen in den Klammern sind die lij . Man erhält für die Zerle-
gungsmatrizen L und R und den Vektor r :
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
1 0 0 2 1 1 4
L=⎝ 1 1 0 ⎠, R = ⎝ 0 2 −1 ⎠ , r =⎝ 1 ⎠
2 −3 1 0 0 2 2
so dass aus Rx = r rekursiv die Lösung x = (1, 1, 1)T folgt.
4.5 Der Gauß-Algorithmus 147
Gegeben: Ax = a, det A = 0.
Gesucht: Lösung x .
Ganz analog lassen sich die letzten beiden Algorithmen unter Verwendung der skalierten
Spaltenpivotsuche formulieren. Dann ist lediglich noch die in Abschnitt 4.5.2 angegebene
Skalierung gemäß Formel (4.10) zu beachten.
Beispiel 4.54.
⇒ 1(0) 3m 3 6 −3
(0)
2 2 2 3 1 — —
(0)
3 1 0 1 0
(1) 2
2 0 0 −1 3 − · 1(0) + 2(0)
3m
1
3(1) 0 −1 −1 1 − · 1(0) + 3(0) —
3m
1
⇒ 2(1) 0 −1 −1 1 Vertauschung von l21 = = 1
3 3
2
3(1) 0 0 −1 3 2. und 3. Zeile l31 = 3 = 2
3
⎛ ⎞
0
⇒ 3(2) 0 0 −1 3 − ⎝ ⎠ · 2(1) + 3(1) l32 = 0
−1 =0
−1
4.5 Der Gauß-Algorithmus 149
2
3 0 1 0 0 −1 3
so dass sich aus R x = r rekursiv die Lösung x T = (3, 2, −3) errechnen lässt. Wegen der
Zeilenvertauschung gilt auch die Zerlegung
⎛ ⎞⎛ ⎞ ⎛ ⎞
1 0 0 3 3 6 3 3 6
LR = ⎝ 3 1 0 ⎠ ⎝ 0 −1 −1 ⎠ = ⎝ 1 0 1 ⎠
1
2
0 1 0 0 −1 2 2 3
⎛3 ⎞⎛ ⎞
1 0 0 3 3 6
= ⎝ 0 0 1 ⎠⎝ 2 2 3 ⎠ = PA
0 1 0 1 0 1
mit einer Permutationsmatrix P, die aus der Einheitsmatrix E durch Vertauschung der
2. und 3. Zeile entsteht.
Liegen Systeme mit gleicher Matrix A und m rechten Seiten aj , j = 1(1)m, vor, so kann
man statt Axj = aj schreiben
wobei A = (aik ), i, k = 1(1)n, die gemeinsame Matrix der m Systeme ist, X die (n, m)-
Matrix, die spaltenweise aus den m Lösungsvektoren xj , j = 1(1)m, aufgebaut ist und
A∗ die (n, m)-Matrix, deren m Spalten die m rechten Seiten aj sind. Die Dreiecks-
zerlegung der Matrix A braucht also nur einmal für alle m Systeme durchgeführt zu
werden, während Vorwärts- und Rückwärtselimination m-mal zu machen sind. Es sind
n3 /3 − n/3 + mn2 Punktoperationen erforderlich. Zusammengefasst ergibt sich damit
folgender
150 4. Lösung linearer Gleichungssysteme
Mit Spaltenpivotsuche muss analog mit der Permutationsmatrix P (vgl. Satz 4.12 bzw.
Algorithmus 4.44) multipliziert werden.
Beispiel 4.56.
Gegeben sind A x1 = a1 und A x2 = a2 mit
⎛ ⎞ ⎛ ⎞ ⎞ ⎛
3 3 6 −3 −6
A=⎝ 2 2 3 ⎠, a1 = ⎝ 1 ⎠ , a2 = ⎝ −3 ⎠
1 0 1 0 −2
⇒ 1(0) 3 3 6 −3 −6
2(0) 2 2 3 1 −3 —
3(0) 1 0 1 0 −2
2
2(1) 0 0 −1 3 1 − · 1(0) + 2(0)
31
3(1) 0 −1 −1 1 0 − · 1(0) + 3(0)
3
1
⇒ 2(1) 0 −1 −1 1 0 − · 1(0) + 3(0)
32
3(1) 0 0 −1 3 1 − 3 · 1(0) + 2(0)
⇒ 3(2) 0 0 −1 3 1 − 0
−1 · 2(1) + 3(1)
Axi = ei , i = 1(1)n,
mit det A = 0, ei i-ter Einheitsvektor. Fasst man die n rechten Seiten ei zu der Einheits-
matrix E zusammen und die n Lösungsvektoren xi zu einer Matrix X , so lassen sich die
n Systeme kompakt in der Form AX = E schreiben.
Man sollte dieses Verfahren nur anwenden, wenn A−1 explizit gebraucht wird. Auf kei-
nen Fall sollte es zur Lösung von m Systemen Axi = yi , i = 1(1)m, durch xi = A−1 yi
verwendet werden, weil dann 4n3 /3 − n/3 + mn2 Punktoperationen benötigt werden im
Gegensatz zu n3 /3 − n/3 + mn2 bei Anwendung des in Abschnitt 4.5.4 angegebenen
Verfahrens.
Beispiel 4.59.
2 3
Gegeben: Die Matrix A = .
1 4
Lösung: Rechenschema
Gestaffeltes System
2 3 x11 x12 1 0 −1 1 4 −3
= ⇒A = .
0 52 x21 x22 − 12 1 5 −1 2
R A−1 R∗
Beispiel 4.60.
Gegeben: Die Matrix ⎛ ⎞
2 −1 0 0
⎜ −1 2 −1 0 ⎟
A=⎜
⎝ 0 −1
⎟.
2 −1 ⎠
0 0 −1 2
Lösung: Rechenschema
3(1) 0 −1 2 −1 0 0 1 0
4(1) 0 0 −1 2 0 0 0 1
−2
3(2) 0 0 4
3 −1 1
3
2
3 1 0 − 3 · 2(1) + 3(1)
4(2) 0 0 −1 2 0 0 0 1
−3
4(3) 0 0 0 5
4
1
4
2
4
3
4 1 − 4 · 3(2) + 4(2)
Gestaffeltes System
⎛ ⎞⎛ ⎞ ⎛ ⎞
2 −1 0 0 x11 x12 x13 x14 1 0 0 0
⎜ 3
−1 ⎟⎜ ⎟ ⎜ 1
0 0 ⎟
⎜ 0 0 ⎟⎜ x21 x22 x23 x24 ⎟ ⎜ 1 ⎟
⎜ 2
⎟⎜ ⎟=⎜ 2
⎟
⎝ 0 0 4
3 −1 ⎠⎝ x31 x32 x33 x34 ⎠ ⎝ 1
3
2
3 1 0 ⎠
5 1 2 3
0 0 0 x41 x42 x43 x44 1
4
4 4
4
R A−1 R∗
⎛ ⎞
4 3 2 1
1⎜ 3 6 4 2 ⎟
⇒ A−1 = ⎜ ⎟.
5⎝ 2 4 6 3 ⎠
1 2 3 4
4.37 und 4.38 aufbauen. Die Verfahren haben nur Sinn ohne Spaltenpivotsuche. Deshalb
ist strenge Regularität stets Voraussetzung.
3. Schritt: (Rückwärtselimination Rx = r )
rn
3.1 xn = rnn
3.2 Für jedes
i = n−1(−1)1
n
xi = rii ri −
1
rik xk
k=i+1
4.7 Verfahren für Systeme mit symmetrischen Matrizen 155
Die Durchführung des Verfahrens verläuft analog zu Algorithmus 4.42. Hier ergeben sich
auch negative Hauptdiagonalelemente der Matrix D (siehe auch Beispiel 4.43).
Mit der Zerlegung A = RT R, wo R = (rik ) eine obere Dreiecksmatrix mit rii > 0 ist,
wird das System Ax = a in ein äquivalentes System Rx = r überführt in folgenden
Schritten:
1. (Faktorisierung) A = RT R ⇒ R,
2. (Vorwärtselimination) a = RT r ⇒ r ,
3. (Rückwärtselimination) Rx = r ⇒ x .
Die Elemente der Matrizen A, R, RT und des Vektors r werden wie folgt bezeichnet
⎛ ⎞
a11 a12 · · · a1n
⎜ a21 a22 · · · a2n ⎟
A=⎜ .
⎝ .. . .. .
..
⎟ = AT , d. h. aik = aki
⎠
an1 · · · an,n−1 ann
⎛ ⎞
r11 r12 · · · r1n
⎜ 0 r22 · · · r2n ⎟
R=⎜ .
⎝ .. .. . ⎟ , rii > 0
⎠
. ..
0 ··· rnn
⎛ ⎞ ⎛ ⎞
r11 0 ··· 0 r1
⎜ r12 r22 0 ⎟ ⎜ r2 ⎟
RT = ⎜ ⎝ ... .. .. ⎟ , rii > 0
⎠ r =⎜ ⎟
⎝ ... ⎠ .
. .
r1n · · · r1,n−1 rnn rn
Aus dem Koeffizientenvergleich A = RT R ergeben sich zeilenweise die Formeln für die
rik . Aus dem Koeffizientenvergleich a = RT r ergeben sich die ri durch Vorwärtselimina-
tion und aus dem Koeffizientenvergleich Rx = r ergeben sich die xi durch Rückwärts-
elimination. Man erhält den
156 4. Lösung linearer Gleichungssysteme
Algorithmus 4.62.
Gegeben: Ax = a mit symmetrischer, positiv definiter Matrix A = (aik ),
i, k = 1(1)n, n ≥ 2, a = (ai ), i = 1(1)n.
Gesucht: x = (xi ), i = 1(1)n.
Dann sind nacheinander folgende Schritte auszuführen:
1. Schritt: (Zerlegung A = RT R)
√
1.1 r11 = a11
1.2 Für jedes j = 2(1)n−1
#
j−1
1.2.1 rjj = ajj − 2
rij
i=1
2. Schritt: (Vorwärtselimination a = RT r )
a1
2.1 r1 = r11
2.2 Für jedes
j j−1
= 2(1)n
rj = aj − rij ri r1jj
i=1
3. Schritt: (Rückwärtselimination Rx = r )
rn
3.1 xn = rnn
3.2 Für jedesi = n−1(−1)1
n
xi = r1ii ri − rik xk
k=i+1
Beispiel 4.63.
1. Weg: Lösung des Systems Ax = a mit positiv definiter Matrix A gemäß Alg. 4.62
1. Schritt: A = RT R (Zerlegung)
⎛ ⎞ ⎛ √ ⎞⎛ √ ⎞
0 −1 2 2 0 − √12
2 √0 0
⎜ √ ⎟
⎝ 0 2 −1 ⎠ = ⎝ 0 2 0 ⎠⎝ 0 2 − √12 ⎠ .
−1 −1 2 − √12 − √12 1 0 0 1
2. Schritt: a = RT r (Vorwärtselimination)
⎛ ⎞ ⎛ √ ⎞⎛ ⎞
1 2 0 − √12 √1
√
⎝ 1 ⎠=⎜ ⎟ 2
⎝ 0 2 − √12 ⎠ ⎝ √12 ⎠ .
0 0 0 1 1
3. Schritt: Rx = r (Rückwärtselimination)
⎛ √ √ ⎞⎛ ⎞ ⎛ √ ⎞
2 √0 − 12 √2 x1 1
2 √2
⎝ 0 2 − 12 2 ⎠ ⎝ x2 ⎠ = ⎝ 1 ⎠,
2 2
0 0 1 x3 1
woraus sich rekursiv die Lösung x3 = x2 = x1 = 1 ergibt.
Die Lösung des Systems Ax = a mit dem Cholesky-Verfahren gemäß Algorithmus 4.62
kann verkürzt wie folgt geschrieben werden
(A, a) = (RT R, RT r ) = RT (R, r ) .
Durch zeilenweisen Koeffizientenvergleich der Matrix (A, a) mit dem Matrizenprodukt
RT (R, r ) können die Matrix R und der Vektor r berechnet werden, danach aus Rx = r
die Lösung.
⎛ ⎞ ⎛ √ ⎞⎛ √ ⎞
2 0 −1 1 2 0 0 2 0 − √12 √1
⎜ ⎟ ⎜ √ √ 2
⎝ 0 2 −1 1 ⎠=⎝ 0 2 0 ⎟ ⎜
⎠⎝ 0 2 − √12 √1
⎟
⎠.
2
−1 −1 2 0 − √2 − √2 1
1 1
0 0 1 1
A a RT R r
Mit der Zerlegung A = RT DR, wobei D eine Diagonalmatrix und R eine normierte
obere Dreiecksmatrix ist, wird das System Ax = a gemäß Algorithmus 4.61 in ein äqui-
valentes System Rx = r übergeführt. Die Zerlegung wird so vorgenommen, dass die
Anzahl der Punktoperationen wie in der 1. Darstellung n3 /6 + O(n2 ) ist.
Algorithmus 4.64.
1. (Faktorisierung A = RT DR)
1.1 d1 = a11
1.2 Für jedes j = 1(1)n:
1.2.1 Für jedes i = 1(1)j−1
hi = rij di (Zwischenspeicher) für j > 1
j−1
1.2.2 dj = ajj − hi rij für j = 2(1)n
i=1
2. (Vorwärtselimination RT z = a, D r = z )
2.1 z1 = a1
j−1
2.2 zj = aj − rij zi für j = 2(1)n
i=1
2.3 rj = zj /dj für j = 1(1)n
4.7 Verfahren für Systeme mit symmetrischen Matrizen 159
3. (Rückwärtselimination Rx = r )
3.1 xn = rn
n
3.2 xj = rj − rji xi für j = n−1(−1)1
i=j+1
Algorithmus 4.65.
Gegeben: Ax = a mit symmetrischer, positiv definiter Matrix A.
Gesucht: x = (xi ), i = 1(1)n.
1. (Faktorisierung A = RT DR)
Für jedes j = 1(1)n
1.1 Für jedes i = 1(1)j−1
1.1.1 h = aij
1.1.2 rij = h/di
1.1.3 Für jedes k = i+1(1)j
akj wird durch akj − hrik ersetzt
1.2 dj = ajj
2. (Vorwärtselimination RT z = a , D r = z )
Für jedes j = 1(1)n
2.1 zj = aj
2.2 Für jedes i = 1(1)j−1
zj := zj − rij zi
2.3 rj = zj /dj
3. (Rückwärtselimination Rx = r )
Für jedes j = n(−1)1
3.1 xj = rj
3.2 Für jedes i = j+1(1)n
xj := xj − rji xi
Ist die Matrix A symmetrisch und positiv definit, so kann auch das Verfahren der konju-
gierten Gradienten (CG-Verfahren) angewandt werden. Anstelle des linearen Gleichungs-
systems Ax = a wird hier iterativ die äquivalente Minimierungsaufgabe gelöst. Es gilt
der
Satz 4.66.
Die folgenden Aufgaben sind äquivalent:
(i) löse Ax = a und (ii) minimiere F (x ) = 12 x T Ax − x T a.
Die Richtung des stärksten Abstiegs einer Funktion ist durch ihren negativen Gradienten
gegeben.
Für F gilt : grad F (x ) = Ax − a.
Es besteht also die Aufgabe, jenen Punkt zu suchen, in dem der Gradient verschwindet.
T
gk+1 Ad k
2.4 βk = dkT Ad k
Für k = n − 1 (d. h. im n-ten Schritt) erhält man Un (d0 , d1 , . . . , dn−1 ) = IRn , denn
dn = 0 mit dn ∈
/ Un = IRn kann nicht existieren!
Das CG-Verfahren kann sowohl den direkten Verfahren zugerechnet werden, weil es nach
genau n Schritten bis auf Rundungsfehler die exakte Lösung liefert, als auch den ite-
rativen Verfahren, weil es im Allgemeinen wegen der raschen Konvergenz bereits nach
wenigen Schritten eine ausreichend gute Lösung liefert.
• die rasche Konvergenz; ihre Geschwindigkeit hängt allerdings von der Kondition
der Matrix A ab; je besser die Kondition um so geringer die benötigte Schrittzahl
bei vorgegebener Genauigkeit.
Ein Nachteil des Verfahrens ist die große Empfindlichkeit gegen Rundungsfehler.
Dass das Verfahren im Gegensatz zum Gauß-Algorithmus leicht vektorisierbar und par-
allelisierbar ist, ergibt sich aus der Art der Operationen, die pro Iterationsschritt aus-
zuführen sind. Es sind pro Schritt drei Skalarprodukte zu berechnen und eine Matrix-
Vektor-Multiplikation (Adk ). Da diese Operationen bei herkömmlicher Verarbeitung ca.
97 % der Rechenzeit ausmachen, muss hier auf besonders effektive Berechnung geachtet
werden.
In der Literatur ist eine Variante des CG-Verfahrens, das sogenannte CG-Verfahren
mit Vorkonditionierung zu finden. Man erreicht damit eine Verringerung der Zahl der
Iterationsschritte bei erhöhtem Rechenaufwand pro Schritt. Das Verfahren ist z. B. in
[BUNS1995], S. 156 ff., [SCHW1991], [SCHW1988], [MAES1985], S. 132-133 beschrieben.
Beispiel 4.68.
Gegeben: Das Gleichungssystem Ax = b mit der symmetrischen und positiv definiten
Matrix A sowie der rechte Seite b
5.0000 −2.0000 2.0000
A= , b=
−2.0000 10.0000 2.0000
1. Es werden gesetzt
1.0000 −1.0000 1.0000
x0 = , d0 = , g0 =
1.0000 −6.0000 6.0000
4.7 Verfahren für Systeme mit symmetrischen Matrizen 163
Beispiel 4.69.
α0 = 0.113646 , β0 = 0.029001
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
−1.5002 −1.4111 0.7731
⎜ 0.5454 ⎟ ⎜ −1.0004 ⎟ ⎜ 0.8844 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
x1 = ⎜
⎜ 0.3181
⎟,
⎟ g1 = ⎜
⎜ −2.8644 ⎟
⎟ d1 = ⎜
⎜ 2.6904
⎟
⎟
⎝ −1.2729 ⎠ ⎝ 0.4528 ⎠ ⎝ −1.0328 ⎠
2.1365 −4.3177 4.6077
α1 = 0.169390 , β1 = 0.528165
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
−1.3693 −1.7809 2.1893
⎜ 0.6952 ⎟ ⎜ 2.9576 ⎟ ⎜ −2.4905 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
x2 = ⎜
⎜ 0.7739
⎟,
⎟ g2 = ⎜
⎜ −0.6360 ⎟
⎟ d2 = ⎜
⎜ 2.0570
⎟
⎟
⎝ −1.4479 ⎠ ⎝ 1.8127 ⎠ ⎝ −2.3582 ⎠
2.9170 0.5088 1.9248
α2 = 1.082218
⎛ ⎞ ⎛ ⎞
1.0000 −0.0000
⎜ −2.0000 ⎟ ⎜ 0.0000 ⎟
⎜ ⎟ ⎜ ⎟
x3 = ⎜
⎜ 3.0000
⎟,
⎟ g3 = ⎜
⎜ −0.0000 ⎟
⎟
⎝ −4.0000 ⎠ ⎝ −0.0000 ⎠
5.0000 0.0000
Damit erhält man als Lösung x des Gleichungssystems x3 = (1, −2, 3, −4, 5)T .
Eine Matrix A = (aik ) heißt tridiagonal, falls gilt aik = 0 für |i − k| > 1, i, k = 1(1)n,
n ≥ 3. Ein Gleichungssystem (4.1) bzw. (4.2) mit tridiagonaler Matrix hat die Gestalt
⎛ ⎞ ⎛ ⎞ ⎛ ⎞
a11 a12 x1 a1
⎜ a21 a22 a23 ⎟ ⎜ x2 ⎟ ⎜ a2 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ a32 a33 a34 ⎟ ⎜ x3 ⎟ ⎜ a3 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ .. .. .. ⎟ ⎜ .. ⎟=⎜ .. ⎟. (4.12)
⎜ . . . ⎟ ⎜ . ⎟ ⎜ . ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ an−1,n−2 an−1,n−1 an−1,n ⎠ ⎝ xn−1 ⎠ ⎝ an−1 ⎠
an,n−1 ann xn an
Das System Ax = a lässt sich mit der Zerlegung A = LR, wo L eine bidiagonale
untere Dreiecksmatrix und R eine normierte bidiagonale obere Dreiecksmatrix ist, in
ein äquivalentes System Rx = r überführen. Voraussetzung für die Zerlegbarkeit ohne
Zeilenvertauschung ist strenge Regularität von A, d. h. es muss gelten (siehe Satz 4.9)
Ist diese Voraussetzung verletzt, muss mit Spaltenpivotsuche gearbeitet werden, wodurch
sich jedoch im Allgemeinen die Bandbreite erhöht (vgl. Abschnitt 4.12). Die Überführung
von Ax = a in Rx = r erfolgt in den Schritten:
1. (Faktorisierung) A = LR ⇒ L,R ,
2. (Vorwärtselimination) a = Lr ⇒ r,
3. (Rückwärtselimination) Rx = r ⇒ x.
166 4. Lösung linearer Gleichungssysteme
1. A = LR =⇒
⎛ ⎞ ⎛ ⎞⎛ ⎞
d1 c1 α1 1 γ1
⎜ b2 d2 c2 ⎟ ⎜ β2 α2 ⎟ ⎜ 1 γ2 ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎜ b3 d3 c3 ⎟ ⎜ β3 α3 ⎟⎜ 1 γ ⎟
⎜ ⎟ ⎜ ⎟⎜ 3 ⎟
A=⎜ . . . ⎟ =⎜ . . ⎟ ⎜ . . ⎟
⎜ .. .. .. ⎟ ⎜ .. .. ⎟⎜ .. .. ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎝ bn−1 dn−1 cn−1 ⎠ ⎝ βn−1 αn−1 ⎠ ⎝ 1 γn−1 ⎠
bn dn βn αn 1
⎛ ⎞
α1 α1 γ1
⎜ β2 β2 γ1 + α2 α2 γ2 ⎟
⎜ ⎟
⎜ β β γ + α α γ ⎟
⎜ 3 3 2 3 3 3 ⎟
=⎜ . . . ⎟
⎜ .. .. .. ⎟
⎜ ⎟
⎝ βn−1 βn−1 γn−2 + αn−1 αn−1 γn−1 ⎠
βn βn γn−1 + αn
Koeffizientenvergleich liefert:
α1 = d1 ⇒ α1 = d1 am
c1 = α1 γ1 ⇒ γ1 = c1
α1
bm
βi = b i cm i = 2(1)n
αi γi = ci ⇒ γi = ci
αi
em i = 2(1)n−1
immer zeilenweise, also zuerst am, dann bm, dann dmund emfür i = 2(1)n−1
und anschließend dmfür i = n; cmist nur Umbenennung.
4.9 Gleichungssysteme mit tridiagonaler Matrix 167
2. Vorwärtselimination a = L r
⎛ ⎞ ⎛ ⎞⎛ ⎞
a1 α1 r1
⎜ a2 ⎟ ⎜ b 2 α2 ⎟⎜ r2 ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎜ a3 ⎟ ⎜ b3 α3 ⎟⎜ r3 ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎜ .. ⎟ = ⎜ .. .. ⎟⎜ .. ⎟
⎜ . ⎟ ⎜ . . ⎟⎜ . ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎝ an−1 ⎠ ⎝ bn−1 αn−1 ⎠ ⎝ rn−1 ⎠
an bn αn rn
3. Rückwärtselimination R x = r
⎛ ⎞⎛ ⎞ ⎛ ⎞
1 γ1 x1 r1
⎜ 1 γ2 ⎟⎜ x2 ⎟ ⎜ r2 ⎟
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎜ 1 γ3 ⎟⎜ x3 ⎟ ⎜ r3 ⎟
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎜ .. .. ⎟⎜ .. ⎟=⎜ .. ⎟
⎜ . . ⎟⎜ . ⎟ ⎜ . ⎟
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎝ 1 γn−1 ⎠ ⎝ xn−1 ⎠ ⎝ rn−1 ⎠
1 xn rn
Algorithmus 4.70.
Gegeben: Ax = a mit tridiagonaler Matrix A, det (Ak ) = 0 für k = 1(1)n−1, n ≥ 3.
Gesucht: x = (xi ), i = 1(1)n.
Dann sind nacheinander folgende Schritte auszuführen:
1. (Zerlegung A = LR)
1.1 α1 = d1
1.2 γ1 = c1 /α1
1.3 Für jedes i = 2(1)n−1 sind zu berechnen:
1.3.1 αi = di − bi γi−1
1.3.2 γi = ci /αi
1.4 αn = dn − bn γn−1
168 4. Lösung linearer Gleichungssysteme
2. (Vorwärtselimination a = Lr )
2.1 r1 = a1 /d1
2.2 Für jedes i = 2(1)n sind zu berechnen:
ri = (ai − bi ri−1 )/αi
3. (Rückwärtselimination Rx = r )
3.1 xn = rn
3.2 Für jedes i = n−1(−1)1 sind zu berechnen:
xi = ri − γi xi+1
n3 n
+ n2 −
3 3
Punktoperationen, d. h. zum Beispiel bei n = 1000 benötigt der Gauß-Algorithmus ∼
109 /3 Punktoperationen, während das Verfahren speziell für tridiagonale Matrizen nur
knapp 5 000 Punktoperationen erfordert.
Beispiel 4.71.
Lösung:
1. Zerlegung A = LR
⎛ ⎞ ⎛ ⎞⎛ −1
⎞
2 −1 0 0 2 0 0 0 1 2 0 0
⎜ −1 2 −1 0 ⎟ ⎜ −1 3
0 0 ⎟ ⎜ 0 1 −2
0 ⎟
⎜ ⎟=⎜ 2 ⎟⎜ 3 ⎟
⎝ 0 −1 2 −1 ⎠ ⎝ 0 −1 4
0 ⎠⎝ 0 0 1 −3 ⎠
3 4
0 0 −1 2 0 0 −1 5
4 0 0 0 1
4.9 Gleichungssysteme mit tridiagonaler Matrix 169
2. Vorwärtselimination Lr = a
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ −5 ⎞
2 0 0 0 r1 −5 r1 2
⎜ −1 3
0 0 ⎟ ⎜ r2 ⎟ ⎜ 1 ⎟ ⎜ r2 ⎟ ⎜ −1 ⎟
⎜ 2 ⎟⎜ ⎟=⎜ ⎟ =⇒ ⎜ ⎟=⎜ 9 ⎟
⎝ 0 −1 4
0 ⎠⎝ r3 ⎠ ⎝ 4 ⎠ ⎝ r3 ⎠ ⎝ ⎠
3 4
0 0 −1 54 r4 −1 r4 1
3. Rückwärtselimination Rx = r
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞
1 −12 0 0 x1 −5
2 −2
⎜ 0 1 −2 0 ⎟ ⎜ x2 ⎟ ⎜ −1 ⎟ ⎜ 1 ⎟
⎜ 3 ⎟ ⎜ ⎟=⎜ 9 ⎟ =⇒ x =⎜ ⎟
⎝ 0 0 1 −3 ⎠ ⎝ x3 ⎠ ⎝ ⎠ ⎝ 3 ⎠.
4 4
0 0 0 1 x4 1 1
1. Faktorisierung: A = RT DR ⇒ R und D ,
2. Vorwärtselimination: RT z = a ⇒ z,
Dr = z ⇒ r,
3. Rückwärtselimination: Rx = r ⇒ x.
170 4. Lösung linearer Gleichungssysteme
⎛ ⎞ ⎛ ⎞
d1 c1 1 γ1
⎜ c1 d2 c2 ⎟ ⎜ 1 γ2 ⎟
⎜ ⎟ ⎜ ⎟
⎜ .. .. .. ⎟ ⎜ .. .. ⎟
A=⎜ . . . ⎟, R = ⎜ . . ⎟,
⎜ ⎟ ⎜ ⎟
⎝ cn−2 dn−1 cn−1 ⎠ ⎝ 1 ⎠
γn−1
cn−1 dn 1
⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
α1 r1 z1 x1 a1
⎜ α2 ⎟ ⎜ r2 ⎟ ⎜ z2 ⎟ ⎜ x2 ⎟ ⎜ a2 ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟
D=⎜ .. ⎟, r = ⎜ .. ⎟, z = ⎜ .. ⎟, x = ⎜ .. ⎟ , a = ⎜ .. ⎟
⎝ . ⎠ ⎝ . ⎠ ⎝ . ⎠ ⎝ . ⎠ ⎝ . ⎠
αn rn zn xn an
Algorithmus 4.72.
Gegeben: Ax = a, A symmetrisch, tridiagonal, positiv definit.
Gesucht: x = (xi ), i = 1(1)n.
Dann sind nacheinander folgende Schritte auszuführen:
1. (Zerlegung A = RT DR)
1.1 α1 = d1
1.2 γ1 = c1 /α1
1.3 Für jedes i = 2(1)n−1 sind durchzuführen:
1.3.1 αi = di − ci−1 γi−1
1.3.2 γi = ci /αi
1.4 αn = dn − cn−1 γn−1
2. (Vorwärtselimination RT z = a, D r = z )
2.1 z1 = a1
2.2 Für jedes i = 2(1)n ist zu berechnen:
zi = ai − γi−1 zi−1
2.3 Für jedes i = 1(1)n ist zu berechnen:
ri = zi /αi
3. (Rückwärtselimination Rx = r )
3.1 xn = rn
3.2 Für jedes i = n−1(−1)1 ist zu berechnen:
xi = ri − γi xi+1
Für die Determinante von A gilt
det A = det(RT DR) = det(RT ) det D det R = det D = α1 α2 . . . αn .
4.9 Gleichungssysteme mit tridiagonaler Matrix 171
Beispiel 4.73.
Lösung:
1. Zerlegung A = RT DR
⎛ ⎞ ⎛ ⎞⎛ ⎞⎛ ⎞
2 −1 0 0 1 0 0 0 2 0 0 0 1 − 21 0 0
⎜ −1 2 −1 0 ⎟ ⎜ −1 1 0 0 ⎟ ⎜0 3
0 0 ⎟⎜ 0 1 − 2
0 ⎟
⎜ ⎟=⎜ 2 ⎟⎜ 2 ⎟⎜ 3 ⎟
⎝ 0 −1 2 −1 ⎠ ⎝ 0 − 23 1 0 ⎠ ⎝ 0 0 43 0 ⎠ ⎝ 0 0 1 − 43 ⎠
0 0 −1 2 0 0 − 34 1 0 0 0 45
0 0 0 1
2. Vorwärtselimination RT z = a, D r = z
⎛ ⎞⎛ ⎞ ⎞⎛ ⎛ ⎞
1 0 0 0 z1 −5 −5
⎜ −1 0 ⎟ ⎜ z2 ⎟ ⎜ 1 ⎟ ⎜ 3 ⎟
⎜ 2 1 0 ⎟⎜ ⎟=⎜ ⎟ ⇒ z = ⎜ −2 ⎟
⎝ 0 −2 1 0 ⎠ ⎝ z3 ⎠ ⎝ 4 ⎠ ⎝ 3 ⎠
3
0 0 − 34 1 z4 −1 5
4
⎛ ⎞⎛ ⎞ ⎞⎛ ⎛ 5 ⎞
2 0 0 0 r1 −5 −2
⎜ 0 3
0 0 ⎟ ⎜ r2 ⎟ ⎜ −3 ⎟ ⎜ −1 ⎟
⎜ 2 ⎟⎜ ⎟=⎜ 2 ⎟⇒r =⎜ 9 ⎟
⎝ 0 0 4
0 ⎠ ⎝ r3 ⎠ ⎝ 3 ⎠ ⎝ ⎠
3 4
5 5
0 0 0 4 r4 4 1
3. Rückwärtselimination Rx = r
⎛ ⎞⎛ ⎞ ⎛ 5 ⎞ ⎛⎞
1 − 12 0 0 x1 −2 −2
⎜ 0 1 − 32 0 ⎟ ⎜ x2 ⎟ ⎜ −1 ⎟ ⎜ ⎟
⎜ ⎟⎜ ⎟=⎜ ⎟ ⇒ x = ⎜ 1 ⎟.
⎝ 0 0 3 ⎠⎝
1 −4 x3 ⎠ ⎝ 94 ⎠ ⎝ 3 ⎠
0 0 0 1 x4 1 1
172 4. Lösung linearer Gleichungssysteme
4.10 Gleichungssysteme
mit zyklisch tridiagonaler Matrix
Analog zur Vorgehensweise bei den tridiagonalen Systemen (Algorithmus 4.70) ergeben
sich auch hier durch Koeffizientenvergleich die Elemente von L, R, r und schließlich von
der Lösung x .
4.10 Gleichungssysteme mit zyklisch tridiagonaler Matrix 173
Algorithmus 4.74.
Gegeben: Ax = a mit zyklisch tridiagonaler Matrix A und det (Ak ) = 0 für
k = 1(1)n−1, n ≥ 4.
Gesucht: x = (xi ), i = 1(1)n.
Dann sind nacheinander folgende Schritte auszuführen:
1. (Faktorisierung A = LR)
1.1 α1 = d1
1.2 γ1 = c1 /α1
1.3 δ1 = e1 /α1
1.4 Für jedes i = 2(1)n−2 sind durchzuführen:
1.4.1 αi = di − bi γi−1
1.4.2 γi = ci /αi
1.4.3 βi = b i
1.4.4 δi = −βi δi−1 /αi
1.5 αn−1 = dn−1 − bn−1 γn−2
1.6 βn−1 = bn−1
1.7 ε3 = cn
1.8 Für jedes i = 4(1)n ist zu berechnen:
εi = −εi−1 γi−3
1.9 γn−1 = (cn−1 − βn−1 δn−2 )/αn−1
1.10 βn = bn − εn γn−2
n
1.11 αn = dn − εi δi−2 − βn γn−1
i=3
2. (Vorwärtselimination Lr = a)
2.1 r1 = a1 /α1
2.2 Für jedes i = 2(1)n−1 ist zu berechnen:
ri = (ai − ri−1 βi )/αi
n
2.3 rn = (an − εi ri−2 − βn rn−1 )/αn
i=3
3. (Rückwärtselimination Rx = r )
3.1 xn = rn
3.2 xn−1 = rn−1 − γn−1 xn
3.3 Für jedes i = n−2(−1)1 ist zu berechnen:
xi = ri − γi xi+1 − δi xn
Für die Determinante von A gilt:
Beispiel 4.75.
Gegeben: Das Gleichungssystem Ax = a mit zyklisch tridiagonaler Matrix A und
det(Ak ) =0 für k = 1(1)4
⎛ ⎞⎛ ⎞ ⎛ ⎞
2 −1 0 0 1 x1 5
⎜ −1 2 −1 0 ⎟ ⎜ ⎟ ⎜ ⎟
⎜ 0 ⎟ ⎜ x2 ⎟ ⎜ −8 ⎟
⎜ 0 −1 2 −1 0 ⎟ ⎜ x3 ⎟ = ⎜ 9 ⎟ .
⎜ ⎟⎜ ⎟ ⎜ ⎟
⎝ 0 0 −1 2 −1 ⎠ ⎝ x4 ⎠ ⎝ −6 ⎠
−1 0 0 −1 2 x5 2
Lösung:
1. Faktorisierung A = LR:
⎛ ⎞ ⎛ ⎞⎛ ⎞
2 −1 0 0 1 2 0 0 0 0 1 − 12 0 0 1
2
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎜ −1 2 −1 0 0 ⎟ ⎜ −1 3
0 0 0 ⎟⎜ 0 1 −3
2
0 1
⎟
⎜ ⎟ ⎜ 2 ⎟⎜ 3 ⎟
⎜ 0 −1 2 −1 ⎟ ⎜
0 ⎟ = ⎜ 0 −1 4
0 0 ⎟⎜ 0 0 1 − 34 1 ⎟
⎜ 3 ⎟⎜ 4 ⎟
⎜ ⎟ ⎜ ⎟⎜ ⎟
⎝ 0 0 −1 2 −1 ⎠ ⎝ 0 0 −1 5
4 0 ⎠⎝ 0 0 0 1 − 35 ⎠
−1 0 0 −1 2 −1 − 2 − 13 − 54
1
2 0 0 0 0 1
2. Vorwärtselimination Lr = a
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ 5
⎞
2 0 0 0 0 r1 5 2
⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ −1 3
0 0 0 ⎟⎜ r2 ⎟ ⎜ −8 ⎟ ⎜ −3
11
⎟
⎜ 2 ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ 0 −1 4
0 0 ⎟⎜ r3 ⎟=⎜ 9 ⎟ ⇒ r =⎜ 4 ⎟
⎜ 3 ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ 0 0 −1 5
4 0 ⎠⎝ r4 ⎠ ⎝ −6 ⎠ ⎝ −5
8
⎠
−1 − 2 − 31 − 54
1
2 r5 4 2
3. Rückwartselimination Rx = r
⎛ ⎞⎛ ⎞ ⎛ ⎞ ⎛ ⎞
1 − 21 0 0 1
2 x1 5
2 1
⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ 0 1 − 23 0 1
⎟⎜ x2 ⎟ ⎜ − 11 ⎟ ⎜ −2 ⎟
⎜ 3 ⎟⎜ ⎟ ⎜ 3 ⎟ ⎜ ⎟
⎜ 0 0 1 − 3 1 ⎟⎜ x3 ⎟=⎜ 4 ⎟ ⇒ x =⎜ 3 ⎟
⎜ 4 4 ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ 0 0 0 1 − 53 ⎠⎝ x4 ⎠ ⎝ −5
8
⎠ ⎝ −1 ⎠
0 0 0 0 1 x5 1 1
4.10 Gleichungssysteme mit zyklisch tridiagonaler Matrix 175
Die Matrix A = (aik ), i, k = 1(1)n, sei zyklisch tridiagonal, symmetrisch und positiv
definit.
Durch Koeffizientenvergleich in den drei Schritten ergibt sich der folgende Algorithmus.
176 4. Lösung linearer Gleichungssysteme
Algorithmus 4.76.
Gegeben: Ax = a mit symmetrischer, zyklisch tridiagonaler, positiv definiter
Matrix A.
Gesucht: x = (xi ), i = 1(1)n, n ≥ 4.
Dann sind nacheinander folgende Schritte auszuführen:
1. (Faktorisierung A = RT DR)
1.1 α1 = d1
1.2 γ1 = c1 /α1
1.3 δ1 = cn /α1
1.4 Für jedes i = 2(1)n−2 sind zu berechnen:
1.4.1 αi = di − ci−1 γi−1
1.4.2 γi = ci /αi
1.4.3 δi = −δi−1 ci−1 /αi
1.5 αn−1 = dn−1 − cn−2 γn−2
1.6 γn−1 = (cn−1 − cn−2 δn−2 )/αn−1
n−2
1.7 αn = dn − αi δi2 − cn−1 γn−1
i=1
2. (Vorwärtselimination RT z = a, D r = z )
2.1 z1 = a1
2.2 Für jedes i = 2(1)n−1 ist zu berechnen:
zi = ai − zi−1 γi−1
n−2
2.3 zn = an − δi zi − γn−1 zn−1
i=1
2.4 Für jedes i = 1(1)n ist zu berechnen:
ri = zi /αi
3. (Rückwärtselimination Rx = r )
3.1 xn = rn
3.2 xn−1 = rn−1 − γn−1 xn
3.3 Für jedes i = n−2(−1)1 ist zu berechnen:
xi = ri − γi xi+1 − δi xn
Für die Determinante von A gilt:
Eine Matrix A = (aik ), i, k = 1(1)n, heißt fünfdiagonal, falls gilt aik = 0 für |i − k| ≥ 2,
n ≥ 4, i, k = 1(1)n.