einf¨uhrung in die numerik vorlesungsskript - math.ovgu.de · kapitel 1 fehleranalyse 1.1...
Post on 30-Aug-2019
9 Views
Preview:
TRANSCRIPT
Einfuhrung in die NumerikVorlesungsskript
Prof. Dr. Lutz Tobiska
Sommersemester 2013
2
Einleitung
Aufgabenstellung der numerischen Mathematik ist die Entwicklung von Metho-den, mit deren Hilfe Losungen mathematischer Problemstellungen effektiv be-rechnet bzw. moglichst mit Fehlerangabe approximiert werden konnen. Numeri-sche Methoden sind der Schlussel zur Simulation komplexer Naherungsvorgangeauf Rechneranlagen. Man mochte teure Experimente wie z.B. Windkanalversu-che bei der Flugzeugkonstruktion oder Festigkeitstests bei Betonkonstruktionendurch beliebig oft und schnell wiederholbare Modellrechnungen ersetzen oder zu-mindest reduzieren. Dabei sind die verwendeten Verfahren aus einfachen Bau-steinen zusammengesetzt (z.B. Integralbestimmung, Losung von Gleichungssys-temen, Nullstellenberechnung, usw.). Wir werden uns in der Vorlesung Numerikvor allem mit diesen einfachen Bausteinen befassen.
Zur numerischen Losung eines Problems aus der Praxis gehort auch eine In-formation uber den gemachten Fehler, um das Resultat richtig einschatzen zukonnen. Der Gesamtfehler setzt sich zusammen aus den
Modellfehlern
• Idealisierungsfehler: Zur Beschreibung eines physikalischen Sachverhalteswird ein mathematisches Modell gebildet. Bei der Formulierung werden da-bei Vereinfachungen angenommen, etwa Vernachlassigung kapillarer Krafteoder linearisierte Materialgesetzte.
• Datenfehler: Die Daten eines mathematischen Modells (z.B. Koeffizienteneiner Differentialgleichung) sind aufgrund ungenauer Kenntnis der Material-eigenschaften notwendig mit Fehlern behaftet.
Numerischer Fehler
• Diskretisierungsfehler: Kontinuierliche Prozesse werden durch endliche Pro-zesse ersetzt (z.B. Approximation der Differentialgleichung durch eine Dif-ferenzengleichung)
• Abbruchfehler: Unendliche Algorithmen werden nach endlich vielen Schrit-ten abgebrochen (z.B. Fixpunktberechnung xn+1 = P (xn))
• Rundungsfehler: Auf einer Rechneranlage mussen alle Rechnungen auf ei-nem endlichen Zahlenbereich durchgefuhrt werden (z.B. 1
3≈ 0.3333).
3
4
Inhaltsverzeichnis
1 Fehleranalyse 71.1 Zahldarstellung und Rundungsfehler . . . . . . . . . . . . . . . . 71.2 Weitere Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Lineare Gleichungssysteme I 112.1 Fehlerabschatzungen . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Direkte Losungsverfahren . . . . . . . . . . . . . . . . . . . . . . . 212.3 Spezielle Gleichungssysteme . . . . . . . . . . . . . . . . . . . . . 282.4 Nicht regulare Systeme . . . . . . . . . . . . . . . . . . . . . . . . 34
3 Interpolation 373.1 Polynominterpolation . . . . . . . . . . . . . . . . . . . . . . . . . 383.2 Interpolationsfehler . . . . . . . . . . . . . . . . . . . . . . . . . . 423.3 Hermite-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 463.4 Spline-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4 Numerische Integration 594.1 Beispiele interpolatorischer Quadraturformeln . . . . . . . . . . . 594.2 Newton-Cotes-Formeln . . . . . . . . . . . . . . . . . . . . . . . . 634.3 Gaußsche Quadraturformeln . . . . . . . . . . . . . . . . . . . . . 65
5 Approximation 735.1 Gauß-Approximation . . . . . . . . . . . . . . . . . . . . . . . . . 735.2 Tschebyscheff-Approximation . . . . . . . . . . . . . . . . . . . . 78
6 Nichtlineare Gleichungen 856.1 Nullstellen reeller Funktionen . . . . . . . . . . . . . . . . . . . . 856.2 Konvergenzverhalten iterativer Verfahren . . . . . . . . . . . . . . 906.3 Interpolationsverfahren . . . . . . . . . . . . . . . . . . . . . . . . 916.4 Newton-Verfahren im Rd . . . . . . . . . . . . . . . . . . . . . . . 95
7 Lineare Gleichungssysteme II 977.1 Einzelschritt- und Gesamtschrittverfahren . . . . . . . . . . . . . 977.2 Abstiegsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5
6 INHALTSVERZEICHNIS
Kapitel 1
Fehleranalyse
1.1 Zahldarstellung und Rundungsfehler
Die Verarbeitung numerischer Algorithmen auf digitalen Rechenanlagen fuhrtzwangsweise zu Fehlern, die durch die Endlichkeit des Bereiches der auf einersolchen Maschine darstellbaren Zahlen bedingt ist. Sei b ∈ N, b ≥ 2 eine festeBasis. Dann besitzt jede reelle Zahl x bezuglich der Basis b eine Entwicklung derForm:
x = ±bE∞∑
i=1
mib−i , E =
s−1∑
j=0
ejbj
mit Koeffizienten mi, ej ∈ 0, 1, . . . , b−1 und einem Exponenten E ∈ Z. Die Ein-deutigkeit der Zahlendarstellung folgt durch die Festlegung m1 6= 0 undmi < b−1fur unendlich viele i und gilt fur alle x ∈ R\0.
Zur Approximation reeller Zahlen werden auf Rechneranlagen sogenannteGleitkommazahlen und Gleitkommaoperationen verwendet. Eine normalisierteGleitkommazahl (zur Basis b) ist eine reelle Zahl a in der Form
a = ±M · bE
mit der Mantisse M = 0.m1....mr und dem Exponenten E ∈ Z, wobei mi ∈0, · · · , b− 1. Fur a 6= 0 ist die Darstellung durch die Normierung m1 6= 0 ein-deutig. Fur a = 0 setzt man M = 0 und E beliebig.
Zur Darstellung solcher normalisierten Gleitkommazahlen brauchen wir also
r Ziffern + 1 Vorzeichen fur die Mantisse
s Ziffern + 1 Vorzeichen fur den Exponenten.
Gebrauchliche Basen sind b = 2 (Dualsystem), b = 10 (Dezimalsystem) und b =16 (Hexadezimalsystem). Die auf dem Rechner darstellbaren Gleitkommazahlen
7
8 KAPITEL 1. FEHLERANALYSE
heißen auch Maschinenzahlen, wegen der Endlichkeit gibt es eine großte (kleinstedarstellbare Zahl)
±(b− 1)b−1 + b−2 + · · ·+ b−r · b(b−1)(bs−1+···+b0) = ±(1− b−r)bbs−1
sowie eine kleinste positive und eine großte negative Zahl
aposmin = b−1 · b−(bs−1) = b−bs
, anegmax = −b−bs
MATLAB verwendet intern b = 2, r = 53, der Exponent variiert zwischen −1021und +1024. Betrachten wir das IEEE–Format, das 64 Bits verwendet:
a = ±M · 2c−1022
Bits Verwendung1 Vorzeichen
52 Mantisse M = 2−1 +m22−2 + . . .+m532
−53
(die erste Mantissenstelle ist immer 1 aus Normierungsgrunden, mitAusnahme der Null)
11 Charakteristik C = c020 + c12
1 + . . .+ c10210 ∈ (0, 2047)
Die ausgenommenen Werte C = 0 und C = 2047 der Charakteristik werdenzur Darstellung der Null (m2 = . . . = m53 = 0, c0 = . . . = c10 = 0) sowie derSondergroße NaN = Not a Number verwendet. Zahlen außerhalb des zulassigenBereiches
D := [amin, anegmin] ∪ 0 ∪ [aposmin, amax]
werden auf NaN abgebildet bzw. als overflow registriert. Fur x ∈ D wird eineRundungsoperation durchgefuhrt,
rd : D → (Menge der normierten Gleitkommazahlen) = F
|rd(x)− x| = minf∈F
|x− f | ∀x ∈ D.
Im IEEE–Format bedeutet dies
rd(x) = sgn(x)
0.m1 · · ·m53 · 2E falls m54 = 0,
(0.m1 . . .m53 + 2−53)2E falls m54 = 1.
Damit ist der absolute Fehler
|rd(x)− x| ≤ b−r
2bE = 2−54 · 2E
vom Exponenten E abhangig, der relative Fehler
∣∣∣∣rd(x)− x
x
∣∣∣∣ ≤1
2
b−rbE
|M |bE ≤1
2
b−rbE
b−1bE=
1
2b1−r
1.2. WEITERE BEISPIELE 9
ist beschrankt fur x ∈ D, x 6= 0 durch die Maschinengenauigkeit
eps =1
2b−r+1.
Fur x ∈ D ist dann
rd(x) = x · (1 + ε) mit |ε| ≤ eps.
Das IEEE–Format liefert
eps ≤ 1
22−52 ≈ 10−16.
Arithmetische Grundoperationen +,−, ·, \ werden auf der Rechenanlage durchentsprechende Maschinenoperationen ⊕,⊖,⊙,⊘ ersetzt, welche Maschinenzahlenin Maschinenzahlen uberfuhren. Dazu werden die Operationen maschinenintern(oft unter einer erhohten Stellenzahl fur die Mantisse) ausgefuhrt, in die normali-sierte Form uberfuhrt und dann gerundet. Liegt das Ergebnis nicht in D wird eineFehlermeldung ausgegeben. Man beachte, dass Assoziativ- bzw. Distributivgesetzim allgemeinen fur Maschinenoperationen nicht gelten.
1.2 Weitere Beispiele
a) Betrachten wir das lineare Gleichungssystem
0.780x1 + 0.563x2 = 0.2170.913x1 + 0.659x2 = 0.254
mit der eindeutigen Losung (x1, x2) = (1,−1). Wir nehmen die Basis b = 10und eine Mantissenlange von r = 3 an. Die Cramersche Regel ist nicht an-wendbar, denn det(A) = 0.
b) Wir wollen ln 2 berechnen und erinnern uns, dass
ln(1 + x) =
∞∑
n=1
(−1)n−1xn
n∀x ∈ (−1, 1].
Somit ist
ln 2 =
∞∑
n=1
(−1)n−1
n=
m∑
n=1
(−1)n−1
n+Rm.
10 KAPITEL 1. FEHLERANALYSE
Nach dem Leibnitzschen Konvergenzkriterium bzw. der Fehlerabschatzungfur alternierende Reihen ist |Rm| ≤ 1/(m+ 1) < 1/m. Somit folgt
∣∣∣ ln 2 − rd
(m∑
n=1
(−1)n−1
n
) ∣∣∣
≤ |Rm|︸︷︷︸Abbruchfehler
+
∣∣∣∣∣
m∑
n=1
(−1)n−1
n− rd
(m∑
n=1
(−1)n−1
n
)∣∣∣∣∣︸ ︷︷ ︸
Rundungsfehler
≤ 1
m+m · eps
Wir beobachten eine fur die Numerik typische Situation; mit zunehmendenm nimmt der Abbruchfehler ab, der Rundungsfehler aber zu. Das Minimumwird fur m2 = eps−1 angenommen. Der Gesamtfehler bleibt damit nur auf
1
m+m · eps ≈ 2
√eps
beschrankt (kann also nicht auf eine beliebig kleine Schranke gedruckt wer-den). Fur die numerische Bestimmung von ln 2 ist die folgende Darstellungbesser geeignet. Aus
ln(1 + x) =
∞∑
n=1
(−1)n−1xn
n, ln(1− x) = −
∞∑
n=1
xn
nx ∈ (−1,+1)
folgt durch Subtraktion
ln1 + x
1− x =
∞∑
n=1
[(−1)n−1 + 1]xn
n= 2
∞∑
k=1
x2k−1
2k − 1.
Fur x = 13
erhalten wir
ln 2 = 2
∞∑
k=1
1
2k − 1
(1
3
)2k−1
= 2
m∑
k=1
1
2k − 1
(1
3
)2k−1
+ Rm
mit der Fehlerschranke
|Rm| ≤(
2
2m+ 1
)(1
3
)2m+1 ∞∑
n=0
(1
3
)2n
=
(2
2m+ 1
)(1
3
)2m+11
1− 1/9
=2
3· 98· 1
2m+ 1· 9−m ≤ 10−7 fur m ≥ 6.
Wir sehen, dass die neue Reihendarstellung wesentlich schneller gegen ln 2konvergiert.
Kapitel 2
Losung linearerGleichungssysteme I
Seien A eine Matrix und b ein Vektor
A = (aj k) j = 1, · · · , m
k = 1, · · · , n
=
a1 1 · · · a1 n...
...am 1 · · · am n
, b = (bj)j=1,··· ,m =
b1...bm
.
Gesucht ist ein Vektor x = (xk)k=1,··· ,n mit der Eigenschaft
a11x1 + a12x2 + · · ·+ a1nxn = b1a21x1 + a22x2 + · · ·+ a2nxn = b2
...am1x1 + am2x2 + · · ·+ amnxn = bm
oder abgekurzt geschrieben: Ax = b. Das lineare Gleichungssystem Ax = b heißtunterbestimmt im Fall m < n, quadratisch im Fall m = n und uberbestimmtim Fall m > n. Es ist genau dann losbar, wenn Rang(A)=Rang[A, b], mit derzusammengesetzten Matrix
[A, b] =
a1 1 · · · a1 n b1...
......
am 1 · · · am n bm
.
Im quadratischen Fall sind die folgenden Aussagen aquivalent:
i) Ax = b ist fur jedes b eindeutig losbar.
ii) Rang(A) = n.
iii) det(A) 6= 0
11
12 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
iv) Alle Eigenwerte von A sind ungleich Null.
Wir beschaftigen uns im folgenden hauptsachlich mit der Losung von quadrati-schen Gleichungssystemen. Die dazu verwendeten Verfahren lassen sich grob inzwei Klassen einteilen: Ein direktes Verfahren zur Losung des GleichungssystemsAx = b ist ein Algorithmus, der (bei Vernachlassigung von Rundungsfehlern) inendlich vielen Schritten die Losung x liefert. Im Gegensatz dazu erzeugen die ite-rativen Verfahren sukzessive eine Folge von Vektoren x(t)t=1,2,···, die im Limesfur t→∞ immer bessere Approximationen zur Losung x sind.
2.1 Fehlerabschatzungen
Wir beschaftigen uns zunachst mit dem Problem der Konditionierung von qua-dratischen linearen Gleichungssystemen. Bei der Losung eines GleichungssystemsAx = b treten zwei Fehlereinflusse ein:
a) Fehler in der theoretischen Losung aufgrund von Eingangsfehlern in denElementen von A und b,
b) Fehler in der numerischen Losung aufgrund des Rundungsfehlers im Ver-laufe des Losungsprozesses.
Zur Erfassung dieser Fehler benotigen wir ein Maß fur die Große von Vektorenund Matrizen. Dazu dienen ublicherweise Normen auf dem n–dimensionalen Zah-lenraum Kn, K = R oder K = C. (Im Hinblick auf spatere Anwendungen lassenwir im folgenden auch komplexe Vektoren bzw. Matrizen zu.) Eine Abbildung‖ · ‖ : Kn → R+ heißt Norm, wenn sie folgende Eigenschaften besitzt:
(N1) ‖x‖ > 0, x ∈ Kn\0 (Definitheit),
(N2) ‖αx‖ = |α| · ‖x‖, x ∈ Kn, α ∈ K (positive Homogenitat),
(N3) ‖x+ y‖ ≤ ‖x‖+ ‖y‖, x, y ∈ Kn (Subadditivitat).
Beispiele:
‖x‖2 :=
(n∑
i=1
|xi|2) 1
2
euklische Norm (l2–Norm)
‖x‖∞ := maxi=1,··· ,n
|xi| Maximumnorm (l∞–Norm)
‖x‖1 :=
n∑
i=1
|xi| Manhattennorm (l1–Norm)
2.1. FEHLERABSCHATZUNGEN 13
Mit Hilfe einer Norm ‖ · ‖ auf Kn lasst sich die Konvergenz einer Folge vonVektoren gegen einen Vektor erklaren durch
x(t) → x (t→∞) : ⇔ ‖x(t) − x‖ → 0 (t→∞).
Aus der Dreiecksungleichung (N3) folgt uber die Beziehung ‖x‖ = ‖x − y + y‖die wichtige Ungleichung
‖x− y‖ ≥∣∣∣‖x‖ − ‖y‖
∣∣∣, x, y ∈ Kn (2.1.1)
welche u.a. die Stetigkeit von Normen als Funktionen von Kn in R impliziert.
Theorem 2.1.1 Auf dem endlichdimensionalen Vektorraum Kn sind alle Nor-men aquivalent, d.h.: Zu je zwei Normen ‖ · ‖, ‖ · ‖′ gibt es positive Konstantenm, M , mit denen gilt:
m‖x‖ ≤ ‖x‖′ ≤M‖x‖, x ∈ Kn (2.1.2)
Beweis: Es genugt, die Behauptung fur den Fall zu zeigen, dass eine der beidenNormen die Maximumnorm ‖ · ‖∞ ist. Sei ‖ · ‖ irgendeine zweite Norm. Bezuglichder kartesischen Einheitsvektoren e1, · · · , en hat jeder Vektor x ∈ Kn die Darstel-lung
x =n∑
i=1
xiei.
Folglich gilt
‖x‖ ≤ γ‖x‖∞, γ :=n∑
i=1
‖ei‖.
Die Norm ‖ · ‖ ist also auch stetig bezuglich der komponentenweisen Konvergenzvon Vektoren. Die Punktmenge
S ≡ x ∈ Kn, ‖x‖∞ = 1 ⊂ K
n
ist beschrankt und abgeschlossen. Die Norm ‖ · ‖ nimmt also auf S ihr Minimiumund Maximum an. Es existieren also x0, x1 ∈ S, so dass
0 < ‖x0‖ ≤ ‖x‖ ≤ ‖x1‖ <∞ ∀x ∈ S.
Fur beliebiges y ∈ Kn\0 ist y/‖y‖∞ ∈ S und folglich
‖x0‖ ≤‖y‖‖y‖∞
≤ ‖x1‖.
14 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Mit m ≡ ‖x0‖ und M ≡ ‖x1‖ gilt daher
m‖y‖∞ ≤ ‖y‖ ≤M‖y‖∞ ∀y ∈ Kn,
die Norm ‖ · ‖ ist also zur Maximumnorm aquivalent.
Die Beziehung 2.1.2 impliziert, dass die durch eine beliebige Norm induzierteKonvergenz von Vektoren stets aquivalent zur komponentenweisen Konvergenzist.
Wir betrachten nun den Vektorraum der n×n–Matrizen A ∈ Kn×n. Offenbarkann dieser mit dem Vektorraum der n ∗ n–Vektoren identifiziert werden. Somitubertragen sich alle Aussagen fur Vektornormen auf Normen fur Matrizen. Ins-besondere sind alle Normen fur n× n–Matrizen aquivalent, und die Konvergenzvon Matrizen ist die komponentenweise Konvergenz:
A(t) → A (t→∞) ⇔ a(t)jk → ajk (t→∞), j, k = 1, · · · , n.
Eine Norm ‖ · ‖ auf Kn×n heißt vertraglich mit einer Vektornorm ‖ · ‖ auf Kn,wenn gilt:
‖Ax‖ ≤ ‖A‖ · ‖x‖, x ∈ Kn, A ∈ K
n×n.
Sie heißt Matrizennorm, wenn sie submultiplikativ ist:
‖AB‖ ≤ ‖A‖ · ‖B‖, A,B ∈ Kn×n.
Beispielsweise ist die Quadratsummennorm (Frobeniusnorm)
‖A‖F :=
(n∑
j,k=1
|ajk|2) 1
2
eine mit der euklidischen Vektornorm vertragliche Matrizennorm. Fur eine belie-bige Vektornorm ‖ · ‖ auf Kn wird durch
‖A‖ := supx∈Kn\0
‖Ax‖‖x‖ = sup
x ∈ Kn
‖x‖ = 1
‖Ax‖
eine mit ‖ · ‖ vertragliche Matrizennorm erklart. Diese heißt die von ‖ · ‖ erzeugtenaturliche Matrizennorm. Fur naturliche Matrizennormen gilt ‖I‖ = 1.
Theorem 2.1.2 Die naturlichen Matrizennormen zu ‖ · ‖∞ und ‖ · ‖1 sind diemaximale Zeilensumme
‖A‖∞ := max1≤j≤n
n∑
k=1
|ajk|
2.1. FEHLERABSCHATZUNGEN 15
bzw. die maximale Spaltensumme
‖A‖1 := max1≤k≤n
n∑
j=1
|ajk|.
Beweis: Wir fuhren den Beweis nur fur ‖ · ‖∞. Offenbar ist die maximaleZeilensumme ‖ · ‖∞ eine Matrizennorm. Wegen
‖Ax‖∞ = max1≤j≤n
∣∣∣n∑
k=1
ajkxk
∣∣∣ ≤ max1≤j≤n
n∑
k=1
|ajk| max1≤k≤n
|xk| = ‖A‖∞‖x‖∞
ist sie vertraglich mit ‖·‖∞. Im Falle ‖A‖∞=0 ist A = 0, d.h. es gilt trivialerweise
‖A‖∞ = sup‖x‖∞=1
‖Ax‖∞.
Sei also ‖A‖∞ > 0 und m ∈ 1, . . . , n ein Index mit der Eigenschaft
‖A‖∞ = max1≤j≤n
n∑
k=1
|ajk| =n∑
k=1
|amk|.
Wir setzen fur k = 1, . . . , n: zk ≡ |amk|/amk fur amk 6= 0 und zk ≡ 0 sonst, somitgilt:z = (zk)k ∈ Kn, ‖z‖∞ = 1. Fur v := Az gilt dann
vm =n∑
k=1
amkzk =n∑
k=1
|amk| = ‖A‖∞.
Folglich ist
‖A‖∞ = vm ≤ ‖v‖∞ = ‖Az‖∞ ≤ sup‖y‖∞=1
‖Ay‖∞.
Die Eigenwerte λ ∈ K einer Matrix A ∈ Kn×n sind die Nullstellen ihrescharakteristischen Polynoms p(λ) = det(A − λ I). Folglich existieren genau n(ihrer Vielfachheit als Nullstelle entsprechend oft gezahlte) Eigenwerte λ, und zujedem λ existiert mindestens ein Eigenvektor w ∈ Kn\0 : Aw = λw. Sei nun‖ · ‖ eine beliebige Vektornorm und ‖ · ‖ eine damit vertragliche Matrizennorm,(wobei die beiden Normen der Einfachheit halber gleich bezeichnet werden). Miteinem normierten Eigenvektor zum Eigenwert λ gilt
|λ| = |λ|‖w‖ = ‖λw‖ = ‖Aw‖ ≤ ‖A‖ ‖w‖ = ‖A‖,
16 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
d.h. alle Eigenwerte von A liegen in einer Kreisscheibe in C mit Mittelpunkt Nullund Radius ‖A‖. Speziell mit ‖A‖∞ erhalt man die Abschatzung
max |λ| ≤ ‖A‖∞ = max1≤j≤n
n∑
k=1
|ajk|.
Eine Matrix A ∈ Kn×n heißt hermitesch, wenn gilt:
A = AT
bzw. ajk = akj, j, k = 1, . . . , n.
Reelle hermitesche Matrizen werden symmetrisch genannt. Der Begriff der Sym-metrie ist eng verknupft mit dem des Skalarprodukts. Eine Abbildung (·, ·):Kn × Kn → K wird ein Skalarprodukt genannt, wenn sie folgende Eigenschaf-ten hat:
(S1) (x, y) = (y, x), x, y ∈ Kn (Symmetrie),
(S2) (αx+ βy, z) = α(x, y) + β(y, z), x, y, z ∈ Kn, α, β ∈ K (Linearitat),
(S3) (x, x) > 0, x ∈ Kn\0 (Definitheit).
Jedes Skalarprodukt auf Kn ×Kn erzeugt durch
‖x‖ := (x, x)12 , x ∈ K
n,
eine zugehorige Vektornorm. Im folgenden wird fast ausschließlich das euklidischeSkalarprodukt verwendet:
(x, y)2 =n∑
j=1
xjyj, (x, x)2 = ‖x‖22.
Mit Hilfe des euklidischen Skalarprodukts laßt sich die Eigenschaft einer Matrix,hermitesch zu sein, aquivalent ausdrucken durch:
A = AT ⇔ (Ax, y)2 = (x,A y)2, x, y ∈ K
n.
Die von der euklidischen Vektornorm erzeugte naturliche Matrizennorm heißt dieSpektralnorm und wird mit ‖ · ‖2 bezeichnet.
Theorem 2.1.3 Fur hermitesche Matrizen gilt
‖A‖2 = max|λ|, λ Eigenwert von A. (2.1.3)
2.1. FEHLERABSCHATZUNGEN 17
Beweis: Bekanntlich besitzt eine hermitesche Matrix A ∈ Kn×n nur reelle Eigen-werte und zwar genau n (ihrer Vielfachheit entsprechend oft gezahlt), λ1, . . . , λn ∈R. Ferner existiert ein zugehoriges Orthonormalsystem von Eigenvektoren
w1, . . . , wn ⊂ Kn : Awi = λwi, (wi, wj)2 = δij, i, j = 1, . . . , n.
Jedes x ∈ Kn besitzt eine Darstellung der Form
x =
n∑
i=1
αiwi, αi = (x, wi)2,
und es gilt
‖x‖22 = (x, x)2 =
n∑
i,j=1
αiαj(wi, wj)2 =
n∑
i=1
|αi|2,
‖Ax‖22 = (Ax,Ax)2 =
n∑
i,j=1
λiαiλjαj(wi, wj)2 =
n∑
i=1
λ2i |αi|2.
Hiermit folgt
‖A‖2 ≤ max1≤i≤n
|λi|.
Wegen der allgemeinen Eigenwertschranke max |λ| ≤ ‖A‖ fur beliebige vertragli-che Matrizennormen folgt damit die Behauptung.
Fur allgemeine Matrizen A ∈ Kn×n gilt
‖A‖2 = max|λ| 12 , λ Eigenwert von ATA.
Lemma 2.1.4 Eine Matrix A ∈ Cn×n ist genau dann hermitisch, wenn (Ax, x)fur alle x ∈ Cn reell ist.
Beweis: Wir nutzen die Tatsache, dass fur jede Matrix A ∈ Cn×n
(Ax, y) = (x,ATy) ∀x, y ∈ C
n gilt.
a) Sei nun A hermitisch, d.h. AT
= A. Setze y = x und nutze die Eigenschaftdes Skalarproduktes (x, y) = (y, x). Dann folgt
(Ax, x) = (x,Ax) = (Ax, x).
b) Sei nun (Ax, x) reell, d.h.
(Ax, x) = (Ax, x) = (x,ATx) = (A
Tx, x)
18 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
bzw.((A − AT
) x, x)
= 0 fur alle x ∈ Cn. Wir nutzen ein Hilfsresultat furbeliebige B ∈ Cn×n
(B x, x) = 0 ∀x ∈ Cn ⇒ B = 0
woraus A = AT
folgen wurde. Wahlt man zunachst
x = (0, · · · , 0, 1︸︷︷︸k-te Stelle
, 0, · · · , 0) ,
so folgt das Verschwinden der Diagonale von B:
(B x, x) =∑
i,j=1
bij xjxi = bkk = 0 ∀k.
Die Wahl
x = (0, · · · , 0, 1︸︷︷︸k-te
, 0, · · · , 0, i︸︷︷︸l-te Stelle
, 0, · · · , 0)← imaginare Einheit
ergibt(B x, x) = bkk + bkl i− blk i− i2 bll = (bkl − blk) i = 0.
Schließlich fuhrt die Wahl
x = (0, · · · , 0, 1︸︷︷︸k-te
, 0, · · · , 0, 1︸︷︷︸l-te Stelle
, 0, · · · , 0)
auf(B x, x) = bkk + bkl + blk + bll = bkl + blk = 0.
Die Losung des Gleichungssystems fur bkl und blk ergibt bkl = blk = 0 furalle l 6= k, l, k ∈ 1, · · · , n.
Eine Matrix A ∈ Kn×n heißt positiv definit, wenn gilt:
(Ax, x)2 ∈ R, (Ax, x)2 > 0 ∀x ∈ Kn\0.
Eine hermitesche Matrix ist genau dann positiv definit, wenn alle ihre (reellen)Eigenwerte positiv sind. Wir werden spater sehen, daß lineare Gleichungssystememit positiv definiten Koeffizientenmatrizen besonders gunstige Losbarkeitseigen-schaften besitzen.
Wir kommen nun zur Fehleranalyse fur lineare Gleichungssysteme
Ax = b
2.1. FEHLERABSCHATZUNGEN 19
mit regularer Koeffizientenmatrix A ∈ Kn×n . Die Matrix A und der Vektor bseien mit Fehlern δA bzw. δb behaftet, so dass ein gestortes System
A x = b,
mit A = A + δA, b = b + δb und x = x + δx gelost wird. Wir wollen den Fehlerδx in Abhangigkeit von δA und δb abschatzen. Dazu sei im folgenden ‖ · ‖ einebeliebige Vektornorm und entsprechend ‖ · ‖ die zugehorige Matrizennorm.
Lemma 2.1.5 Die Matrix B ∈ Kn×n habe die Norm ‖B‖ < 1. Dann ist dieMatrix I +B regular, und es gilt
‖(I +B)−1‖ ≤ 1
1− ‖B‖ . (2.1.4)
Beweis: Fur alle x ∈ Kn gilt
‖(I +B)x‖ ≥ ‖x‖ − ‖B x‖ ≥ (1− ‖B‖)‖x‖.
Wegen 1−‖B‖ > 0 ist also I+B injektiv und folglich regular. Mit der Abschatzung
1 = ‖I‖ = ‖(I +B)(I +B)−1‖ = ‖(I +B)−1 +B(I +B)−1‖≥ ‖(I +B)−1‖ − ‖B‖‖(I +B)−1‖ = ‖(I +B)−1‖(1− ‖B‖) > 0.
erhalt man die behauptete Ungleichung.
Nach diesen Vorbereitungen konnen wir den folgenden allgemeinen Storungs-satz fur lineare Gleichungssysteme beweisen:
Theorem 2.1.6 Die Matrix A ∈ Kn×n sei regular, und es sei
‖δA‖ < ‖A−1‖−1. (2.1.5)
Dann ist die gestorte Martix A = A+ δA ebenfalls regular, und fur den relativenFehler der Losung gilt:
‖δx‖‖x‖ ≤
cond(A)
1− cond(A)‖δA‖\‖A‖
‖δb‖‖b‖ +
‖δA‖‖A‖
, (2.1.6)
mit der Konditionszahl cond(A) von A,
cond(A) := ‖A‖‖A−1‖.
Beweis: Aufgrund der Voraussetzungen ist
‖A−1δA‖ ≤ ‖A−1‖‖δA‖ < 1,
20 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
so dass auch A+ δA = A[I + A−1δA] nach Lemma 2.1.5 regular ist. Aus
(A+ δA)x = b+ δb, (A+ δA)x = b+ δAx
folgt fur δx = x− x zunachst
(A+ δA)δx = δb− δAx.
Somit haben wir
‖δx‖ ≤ ‖(A+ δA)−1‖‖δb‖+ ‖δA‖‖x‖≤ ‖(I + A−1δA)−1‖‖A−1‖‖δb‖+ ‖δA‖‖x‖
≤ ‖A−1‖1− ‖A−1‖‖δA‖‖δb‖+ ‖δA‖‖x‖
≤ ‖A−1‖‖A‖‖x‖1− ‖A−1‖‖A‖‖δA‖\‖A‖
‖δb‖‖A‖‖x‖ +
‖δA‖‖A‖
.
Wegen ‖b‖ = ‖Ax‖ ≤ ‖A‖‖x‖ folgt schließlich
‖δx‖ ≤ cond(A)
1− cond(A)‖δA‖\‖A‖
‖δb‖‖b‖ +
‖δA‖‖A‖
‖x‖
und damit die Behauptung des Satzes.
Die Konditionszahl cond(A) hangt offenbar von der bei ihrer Definition zu-grundegelegten Vektornorm ab. Meistens verwendet man die Maximumnorm ‖·‖∞oder die euklidische Norm ‖ · ‖2. Im ersten Fall ist
cond∞(A) := ‖A‖∞‖A−1‖∞
mit der maximalen Zeilensumme ‖ · ‖∞. Speziell fur hermitesche Matrizen giltnach Lemma 2.1.2
cond2(A) := ‖A‖2‖A−1‖2 =|λmax||λmin|
mit den betragsmaßig großten bzw. kleinsten Eigenwerten λmax und λmin von A;die Große cond2(A) wird auch die Spektralkonditionszahl von A genannt.
Ist cond(A)‖δA‖‖A‖−1 ≪ 1, so wird in Theorem 2.1.6
‖δx‖‖x‖ . cond(A)
‖δb‖‖b‖ +
‖δA‖‖A‖
,
d.h. die Konditionszahl cond(A) ist gerade der Verstarkungsfaktor, mit dem sichdie relativen Fehler in A und b auf den relativen Fehler in x auswirken. Diese
2.2. DIREKTE LOSUNGSVERFAHREN 21
Fehlerabschatzung erlaubt folgenden Schluss:
Regel: Die Kondition von A sei cond(A)∼ 10s. Sind die Elemente von A und bmit einem relativen Fehler der Art
‖δA‖\‖A‖ ∼ 10−k, ‖δb‖\‖b‖ ∼ 10−k(k > s)
behaftet, so muss mit einem relativen Fehler im Ergebnis der Großenordnung
‖δx‖\‖x‖ ∼ 10s−k
gerechnet werden, d.h. man verliert im ungunstigsten Fall s Stellen an Genauig-keit.
Beispiel:
A =
[1.2969 0.86480.2161 0.1441
], A−1 = 108 ·
[0.1441 −0.8648−0.2161 1.2969
]
‖A‖∞ = 2.1617, ‖A−1‖∞ = 1.513 · 108 ⇒ cond(A) ≈ 3.3 · 108.
Bei der Losung des Gleichungssystems Ax = b gehen also im ungunstigen Fall 8wesentliche Stellen an der Genauigkeit, mit der die Elemente ajk und bj gegebensind, verloren. Dieses System ist extrem schlecht konditioniert.
Wir demonstrieren anhand der Spektralkondition, dass die Abschatzung inTheorem 2.1.6 im wesentlichen scharf ist. Sei A eine positiv definite n×n–Matrixmit kleinstem und großtem Eigenwert λ1 bzw. λn sowie zugehorigen normiertenEigenvektoren w1 bzw. wn. Wir wahlen δA ≡ 0, b ≡ wn, δb ≡ ǫw1(ǫ 6= 0).Dann haben die Gleichungen Ax = b und A x = b+ δb die Losungen
x =1
λnwn, x =
1
λnwn + ǫ
1
λ1w1.
Folglich ist fur δx = x− x
‖δx‖2‖x‖2
= ǫλn
λ1
‖w1‖2‖wn‖2
= cond2(A)‖δb‖2‖b‖2
.
2.2 Direkte Losungsverfahren
Im folgenden diskutieren wir direkte Losungsmethoden fur (reelle) quadratischelineare Gleichungssysteme der Form
Ax = b. (2.2.7)
22 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Besonders leicht losbar sind gestaffelte Systeme, z.B. solche mit einer oberenDreicksmatrix A = (ajk) als Koeffizientenmatrix
a11x1 + a12x2 + . . . + a1nxn = b1a22x2 + . . . + a2nxn = b2
...annxn = bn
Im Falle ajj 6= 0, j = 1, . . . , n, erhalt man die Losung durch sog. sukzessivesRuckwartseinsetzen
xn =bnann
, j = n− 1, . . . , 1 : xj =1
ajj
(bj −
n∑
k=j+1
ajkxk
).
Dazu sind offensichtlich n2/2+O(n) arithmetische Operationen erforderlich (1 arith-metische Operation:=1 Multiplikation (+1 Addition) oder 1 Division).
Das klassische direkte Verfahren zur Losung (regularer) Gleichungssysteme istdas Gaußsche Eliminationsverfahren. Dabei wird das gegebene System Ax = bschrittweise in ein oberes Dreiecksystem Rx = c umgeformt, welches dieselbeLosung x besitzt und dann durch Ruckwartseinsetzem gelost wird. Dazu stehendie folgenden elementaren Umformungen zur Verfugung:
(i) Vertauschung zweier Gleichungen,
(ii) Addition des Vielfachen einer Gleichung zu einer anderen.
(Die Vertauschung zweier Spalten von A ist ebenfalls zulassig, wenn die Unbe-kannten xi entsprechend umnumeriert werden.)
In der praktischen Durchfuhrung des Gaußschen Eliminationsverfahrens wen-det man die elementaren Umformungen auf die zusammengesetzte Matrix [A, b]an. Im folgenden wird A als regular angenommen.
Zunachst setzt man A(0) ≡ A, b(0) ≡ b. Bestimme a(0)r1 6= 0, r ∈ 1, . . . , n.
(Solch ein Element existiert, da A sonst singular ware). Vertausche die 1–te unddie r–te Zeile. Das Resultat sei [A(0), b(0)]. Subtrahiere fur j = 2, . . . , n das qj1–fache,
qj1 ≡ a(0)j1 \a
(0)11 (= a
(0)r1 \a(0)
rr ),
der 1–ten Zeile von der j–ten Zeile, das Resultat ist
[A(1), b(1)] =
a(0)11 a
(0)12 . . . a
(0)1n b
(0)1
0 a(1)22 . . . a
(1)2n b
(1)2
......
......
0 a(1)n2 . . . a
(1)nn b
(1)n
.
2.2. DIREKTE LOSUNGSVERFAHREN 23
Ein Zeilen- oder Spaltentausch wird durch Multiplikation mit einer Permuta-tionsmatrix
Pkl =
I0 . . . 1... I
...1 . . . 0
I
beschrieben, wobei I Einheitsmatrizen der Dimensionen k−1, l−1−k, und n− lbezeichnen. Linksmultiplikation PA vertauscht die Zeilen k und l, Rechtsmultipli-kation AP die Spalten k und l. Eine Permutationsmatrix ist eine durch Zeilenper-mutation (Spaltenpermutation) aus der Einheitsmatrix entstehende Matrix. Dieobigen Matrix Pkl ist der Permutation (1, . . . , k−1, l, k+1, . . . , l−1, k, l+1, . . . , n)zugeordnet. Die Determinante einer Permutation ist det P = ±1, je nachdem obdie Permutation gerade (+1) oder ungerade (−1) ist. Ferner gilt P−1 = P . DieElimination der Unbekannten xk in den Zeilen k+1, . . . , n mittels der k-ten Zeilekann als Linksmultiplikation mit der Frobeniusmatrix
Gk =
1. . .
1−qk+1,k 1
.... . .
−qn,k 1
, det Gk = 1
beschrieben werden. Man rechnet nach, dass fur die inverse Matrix
G−1k =
1. . .
1qk+1,k 1
.... . .
qn,k 1
= 2I −Gk
gilt. Den oben beschriebenen Ubergang [A(0), b(0)] → [A(1), b(1)] kann man nunmit Hilfe von Matrizenmultiplikationen beschreiben
[A(0), b(0)] = P1[A(0), b(0)], [A(1), b(1)] = G1[A
(0), b(0)],
wobei P1 eine Permutationsmatrix und G1 eine Frobenius–Matrix sind.
Die Gleichungssysteme Ax = b und A(1)x = b(1) haben offenbar dieselbeLosung:
Ax = b⇔ A(1)x = G1P1Ax = G1P1b = b(1)
24 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Das Element ar1 = a(0)11 heißt ”Pivotelement” und der ganze Teilschritt seiner
Bestimmung ”Pivotsuche”. Aus Grunden der numerischen Stabilitat trifft mangewohnlich die Wahl
|ar1| = max1≤j≤n
|aj1|. (2.2.8)
Der ganze Prozeß inkl. Zeilenvertauschung wird dann ”Spaltenpivotierung” ge-nannt. Sind die Elemente der Matrix A von sehr unterschiedlicher Großenord-nung, so empfielt es sich, totale Pivotierung vorzunehmen. Diese besteht aus derWahl
|ars| = max1≤j,k≤n
|ajk| (2.2.9)
und anschließender Vertauschung der ersten mit der r–ten Zeile und der ers-ten mit der s–ten Spalte. Entsprechend der Spaltenvertauschung mussen danndie Unbekannten xk umnumeriert werden. Bei großen Gleichungssystemen ist dietotalte Pivotisierung meist zu aufwendig, so dass man sich mit der Spaltenpivo-tierung begnugt.
Die im ersten Schritt erzeugte Matrix A(1) ist wieder regular. Dasselbe giltauch fur die um die erste Zeile und Spalte reduzierte Teilmatrix, so dass aufsie der Eliminationsprozeß analog zu Schritt 1 angewendet werden kann. DurchWeiterfuhrung dieses Eliminationsprozesses erhalt man in n − 1 Schritten eineKette von Matrizen
[A, b]→ [A(1), b(1)]→ . . .→ [A(n−1), b(n−1)] =: [R, c],
wobei
[A(i), b(i)] = GiPi[A(i−1), b(i−1)], [A(0), b(0)] := [A, b],
mit Permutationsmatrizen Pi und (regularen) Frobenius–Matrizen Gi sind.
Das Endresultat
[R, c] = Gn−1Pn−1 · · ·G1P1[A, b] (2.2.10)
entspricht einem oberen Dreickssystem Rx = c, welches dieselbe Losung wie dasAusgangssystem Ax = b besitzt.
Im i–ten Eliminationschritt [A(i−1), b(i−1)] → [A(i), b(i)] werden in der i–tenSpalte die Elemente unterhalb der Diagonalen annuliert. Den frei werdenen Platzbenutzt man zur Abspeicherung der wesentlichen Elemente qi+1,i, . . . , qn,i derFrobenius–Matrizen G−1
i (i = 1, . . . , n − 1). Da im i–ten Eliminationschritt dievorausgehenden Zeilen 1 bis i nicht verandert werden, arbeitet man also mit
2.2. DIREKTE LOSUNGSVERFAHREN 25
Matrizen der Form
r11 r12 . . . r1 i r1,i+1 . . . r1 n c1λ21 r22 . . . r2 i r2,i+1 . . . r2 n c2λ31 λ32 r3 i r3,i+1 . . . r3 n c3...
.... . .
......
......
λi 1 λi 2 rii ri,i+1 . . . ri n ciλi+1,1 λi+1,2 λi+1,i a
(i)i+1,i+1 . . . a
(i)i+1,n b
(i)i+1
......
......
......
λn,1 λn,2 . . . λn,i a(i)n,i+1 . . . a
(i)n,n b
(i)n
Dabei sind die Subdiagonalelemente λk+1,k, . . . , λnk der k–ten Spalte gewissePermutationen der wesentlichen Elemente qk+1,k, . . . , qnk von G−1
k , da die Zeilen-vertauschungen (nur diese!) an der gesamten Matrix vorgenommen werden. AlsEndresultat erhalt man eine Matrix
r11 . . . r1 n c1l21 r22 r2 n c2...
. . .. . .
......
ln 1 . . . ln,n−1 rnn cn
.
Theorem 2.2.1 Die Matrizen
L =
1 0l21 1...
. . .. . .
ln 1 . . . ln,n−1 1
und R =
r11 r12 . . . r1 n
r22 . . . r2 n
. . ....
0 rnn
bilden eine sog. LR-Zerlegung der Matrix PA:
PA = LR, P = Pn−1 · · ·P1. (2.2.11)
Beweis: Sei P eine Permutationsmatrix, die nur Zeilen mit einem Index großeroder gleich k + 1 tausche. Dann gilt
Gk = PGkP−1 =
1. . .
1−qk+1,k 1
.... . .
−qn,k 1
,
26 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
dies bedeutet, dass beim Ubergang von Gk auf Gk nur die entsprechenden Zei-leneintrage der qij geandert werden. Wir schreiben nun
R = Gn−1Pn−1Gn−2P−1n−1Pn−1Pn−2Gn−3 (Pn−1Pn−2)
−1 (Pn−1Pn−2Pn−3) · · ·Die Permutationsmatrizen Pn−1, Pn−1Pn−2, Pn−1Pn−2Pn−3 tauschen hochstens dieletzten 2, 3, bzw. 4 Zeilen, somit gilt
R = Gn−1Gn−2 · · · G1Pn−1Pn−2 · · ·P1A.
Nun sind das Produkt Gn−1Gn−2 · · · G1 eine untere Dreiecksmatrix mit Diago-naleintragen Null und P = Pn−1Pn−2 · · ·P1 eine geeignete Permutationsmatrix.Setzen wir
L =(Gn−1Gn−2 · · · G1
)−1
= G−11 G−1
2 · · · G−1n−1,
so erhalten wir die behauptete Darstellung LR = PA.
Theorem 2.2.2 Die LR-Zerlegung (nach Gauß) einer regularen Marix ist, wennsie existiert, eindeutig bestimmt.
Beweis: Sei A = L1R1 = L2R2. Dann ist
L−12 L1 = R2R
−11 = I, L−1
1 L2 = R1R−12 = I,
da L−12 L1, L
−11 L2 untere Dreiecksmatrizen mit Einsen auf der Hauptdiagona-
len und R2R−11 , R1R
−12 obere Dreiecksmatrizen sind. Folglich ist L1 = L2 und
R1 = R2.
Das Gaußsche Verfahren liefert Rx = c aus Ax = b. Aquivalent hierzu istdie Losung zweier Dreieckssysteme, falls man die Dreieckszerlegung PA = LRbereits hat:
PAx = LRx = Pb ⇔ Ly = Pb und Rx = y.
Dies ist vor allem dann sinnvoll, wenn man mehrere Gleichungssysteme mit ver-schiedenen rechten Seiten und gleicher Koeffizientenmatrix A hat.
Die Losung eines n × n Gleichungssystems Ax = b mit Hilfe des GaußschenVerfahrens erfordert
n3
3+ O(n2) arithm. Operationen.
Dasselbe gilt fur die Bestimmung der Dreieckszerlegung PA = LR.
Hinweis: MATLAB-Routinen
[L,R, P ] = lu(A) liefert L, R und P einer gegebenen Matrix A[L∗, R] = lu(A) liefert L∗, R gemaß L∗R = P−1LR = A
2.2. DIREKTE LOSUNGSVERFAHREN 27
Varianten und Anwendungen der Gaußelimination
Nachiteration. Rundungsfehler liefern in der Regel LR 6= PA und damit nureine Naherung x0 mit dem Defekt d0 := b − Ax0 6= 0. Unter Verwendung derbereits erstellten Dreieckszerlegung LR ∼ PA lost man nun (naherungsweise)die sog. Defektgleichung
Ax = d0
und erhalt daraus eine Korrektur x fur x0:
x1 := x0 + x.
Bei exakter Losung der Defektgleichung ware x1 tatsachlich die exakte Losungdes Gleichungssystems, im allgemeinen kann man auch bei nur naherungsweiserLosung der Defektgleichung eine bessere Naherung als x0 erwarten. Dazu mussallerdings der Defekt mit erhohter Genauigkeit berechnet werden. Diese Uberle-gungen stutzen sich auf die folgende Fehleranalyse.
Wir nehmen an, dass sich die relative Storung der Matrix A durch eine kleineZahl ε abschatzen lasst. Nach dem Storungstheorem 2.1.6 gilt
‖x0 − x‖‖x‖ ≤ cond(A)
1− cond(A)‖A− LR‖‖A‖
‖A− LR‖‖A‖︸ ︷︷ ︸∼ε
Der Verlust an Stellen entspricht der Große von cond(A). Zusatzlich auftretendeRundungsfehler werden vernachlassigt. Der exakte Defekt werde durch b − Ax0
ersetzt, wobei A eine genauere Approximation fur A ist,
‖A− A‖‖A‖ ≤ ε≪ ε.
Nach Konstruktion gilt dann
x1 = x0 + x = x0 + (LR)−1[b− Ax0]
= x0 + (LR)−1[Ax− Ax0 + (A− A)x0].
Nun haben wir
x1 − x = x0 − x+ (LR)−1A(x− x0) + (LR)−1(A− A)x0
= (LR)−1(LR− A)(x0 − x) + (LR)−1(A− A)x0.
Wegen
(LR)−1 = (A− A+ LR)−1 =[A(I − A−1(A− LR)
]−1
=(I −A−1(A− LR)
)−1
A−1
28 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
folgt die Abschatzung
‖(LR)−1‖ ≤ ‖A−1‖ 1
1− ‖A−1‖ ‖A− LR‖= ‖A−1‖
(1− cond(A)
‖A− LR‖‖A‖
)−1
.
Dies fuhrt letztlich zu
‖x1 − x‖‖x‖ ≤ cond(A)
‖A− LR‖‖A‖︸ ︷︷ ︸∼ε
‖x0 − x‖‖x‖︸ ︷︷ ︸
∼cond(A)ε
+‖A− A‖‖A‖︸ ︷︷ ︸∼ε
‖x0‖‖x‖
.
In der Praxis reichen oft schon 2-3 Nachiterationen aus, um den Fehler in x aufdie Großenordnung der Genauigkeit der Defektauswertung zu drucken.
Determinantenbestimmung. Ausgehend von der LR-Zerlegung PA = LRerhalten wir mit detP = 1
detA = detP detA = detL detR =n∏
i=1
rii.
Rangbestimmung. Ist der Gaußalgorithmus mit Spaltenpivotierung durchfuhrbar,d.h. lassen sich immer ein nicht verschwindenes Pivotelement finden und ist auchdas letzte Diagonalelement a
(n−1)nn 6= 0, so ist Rang(A) = n. Gilt jedoch im i-ten
Eliminationsschritt fur alle Elemente in der i-ten Spalte
a(i−1)ji = 0, j = i, . . . , n,
so ist A singular. In diesem Fall wird Totalpivotierung vorgenommen (Spaltenund Zeilentausch andern den Rang einer Matrix nicht!). Gilt nun nach dem i-tenIterationschritt
a(i)jk = 0, j, k = i+ 1, . . . , n,
so ist Rang(A) = i.
2.3 Spezielle Gleichungssysteme
Die Losung sehr großer Gleichungssysteme mit dem Gaußschen Eliminationsver-fahren ist mit Schwierigkeiten verbunden, insbesondere wenn der auf dem Com-puter verfugbare Hauptspeicher zur Speicherung nicht ausreicht. Die teilweiseAuslagerung und Verwendung externer Speicher treibt die Rechenzeit aufgrunddes erforderlichen Datentransfers in die Hohe und ist deshalb keine echte Alterna-tive. Viele in der Praxis vorkommende Matrizen besitzen jedoch eine besondereStruktur, die es erlaubt, bei der Durchfuhrung des Gaußschen Verfahrens Spei-cherplatz zu sparen.
2.3. SPEZIELLE GLEICHUNGSSYSTEME 29
Bandmatrizen
Definition 2.3.1 Eine Matrix A ∈ Rn×n heißt Bandmatrix vom Bandtyp (ml, mr),0 ≤ ml, mr ≤ l − 1, wenn
ajk = 0 fur k < j −ml oder k > j +mr, j, k = 1, . . . , n.
Die Elemente von A sind also bis auf die Hauptdiagonale und hochstens ml +mr
Nebendiagonalen gleich Null. Die Große m = ml +mr + 1 heißt Bandbreite.
Untere Dreiecksmatrizen sind vom Typ (n − 1, 0), obere Dreiecksmatrizenvom Type (0, n− 1) und Tridiagonalmatrizen vom Typ (1, 1). Ein Beispiel einer(16× 16)-Matrix vom Bandtyp (4,4) ist:
A =
B II B I
I B II B
mit B =
4 −1−1 4 −1
−1 4 −1−1 4
Die Speicherung einer Bandmatrix A ∈ Rn×n vom Bandtyp (ml, mr) erfolgt ubli-cherweise in der Form einer (ml +mr + 1)× n-Matrix A durch die Zuordnung
ai−j+mr+1,j = ai,j.
Die obige Matrix A in kompakter Speicherung also als
A =
0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 −1 −1 −1 0 −1 −1 −1 0 −1 −1 −1 0 −1 −1 −14 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4−1 −1 −1 0 −1 −1 −1 0 −1 −1 −1 0 −1 −1 −1 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 01 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0
Die Bander werden also in der oberen linken und unteren rechten Ecke mit Nullenaufgefullt und dann “vertikal zusammengeschoben”.
Die Inverse einer Bandmatrix ist im allgemeinen voll besetzt. So gilt fur dieobige (4× 4)-Matrix B
B−1 =1
209
56 15 4 115 60 16 44 16 60 151 4 15 56
.
30 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Fur die LR-Zerlegung benotigt man jedoch keinen zusatzlichen Speicherplatzausserhalb des Bandes, wir haben zum Beispiel B = LR mit
L =
1 0 0 0−0.2500 1 0 0
0 −0.2667 1 00 0 −0.2679 1
R =
4 −1 0 00 3.7500 −1 00 0 3.7333 −10 0 0 3.7321
Theorem 2.3.1 Ist A ∈ Rn×n eine Bandmatrix vom Typ (ml, mr), fur die dasGaußsche Eliminationsverfahren ohne Zeilenvertauschung durchfuhrbar ist, dannsind auch alle reduzierten Matrizen Bandmatrizen desselben Typs, und die Fakto-ren der Dreieckszerlegung von A sind Bandmatrizen vom Typ (ml, 0) bzw. (0, mr).
Beweis: Man uberlegt sich leicht, dass die Frobeniusmatrizen Bandmatrizen vomTyp (ml, 0) sind.
Zur Durchfuhrung der Gaußschen Elimination genugt es also das “Band” zuspeichern, bei Großenordnungen von n ∼ 10000 und m ∼ 100 bedeutet diesbereits eine betrachtliche Reduktion des Speicherplatzbedarfes. Eine extreme Er-sparnis ergibt sich bei Tridiagonalmatrizen
a1 b1
c2. . .
. . .. . .
. . . bn−1
cn an
.
Hier lassen sich die Elemente der LR-Zerlegung
L =
1
γ2. . .. . . 1
γn 1
R =
α1 β1
. . .. . .
αn−1 βn−1
αn
durch einfache, rekursive Beziehungen bestimmen:
α1 = a1, β1 = b1
i = 2, . . . , n− 1 : γi = ci/αi−1, αi = ai − γiβi−1, βi = bi
γn = cn/αn−1, αn = an − γnβn−1
Hierzu sind 3n− 2 Speicherplatze und 2n− 2 Operationen erforderlich.Wesentlich fur Theorem 2.3.1 war, dass das Gaußsche Verfahren ohne Zeilen-
vertauschungen durchgefuhrt werden kann, da andernfalls die Bandbreite anwachst.Wir betrachten im folgenden zwei Klassen von Matrizen, bei denen dies der Fallist.
2.3. SPEZIELLE GLEICHUNGSSYSTEME 31
Diagonaldominante Matrizen
Definition 2.3.2 Eine Matrix A ∈ Rn×n heißt diagonaldominant, wenn
n∑
j=1,j 6=i
|aij| ≤ |aii|, i = 1, . . . , n.
Theorem 2.3.2 Die Matrix A ∈ Rn×n sei regular und diagonaldominant. Dannexistiert eine LR-Zerlegung A = LR, die mit Gaußscher Elimination ohne Pivo-tierung bestimmt werden kann.
Beweis: Wir stellen zunachst fest, dass a11 6= 0 ist, da andernfalls aus der Dia-gonaldominanz
n∑
j=2
|a1j | ≤ |a11| = 0
folgt, also die erste Zeile nur aus Nullelementen besteht und die Matrix nichtregular ware. Der erste Eliminationsschritt A→ A(1) kann somit ohne Pivotierungdurchgefuhrt werden. Die Elemente a
(1)ij werden aus bestimmt aus:
j = 1, . . . , n : a(1)1j = a1j
i = 2, . . . , n, j = 1, . . . , n : a(1)ij = aij − qi,1a1j , qi,1 =
ai,1
a11.
Fur die neuen Zeilen i = 2, . . . , n gilt:
n∑
j=2,j 6=i
|a(1)ij | ≤
n∑
j=2,j 6=i
|aij|+ |qi,1|n∑
j=2,j 6=i
|a1j |
≤n∑
j=1,j 6=i
|aij| − |ai1|+ |qi,1|n∑
j=2
|a1j| − |qi,1| |a1i|
≤ |aii| − |qi,1a1i| ≤ |a(1)ii |
Die Matrix A(1) ist regular und wieder diagonaldominant, und folglich ist a(1)22 6= 0.
Die Eigenschaft bleibt also in jedem Schritt der Gaußschen Elimination erhalten.Wir konnen die Elimination folglich ohne Zeilenvertauschungen durchfuhren.
Positiv definite Matrizen
Theorem 2.3.3 Fur positiv definite Matrizen A ∈ Rn×n ist das Gaußsche Eli-
minationsverfahren ohne Zeilenvertauschung durchfuhrbar, und die dabei auftre-tenden Pivotelemente a
(i)ii sind positiv.
32 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Beweis: Da eine positiv definite Matrix positive Diagonalelemente hat, ist ins-besondere a11 > 0. Die Beziehung
a(1)ij = aij −
ai1a1j
a11
= aij −a1iaj1
a11
= a(1)ij , i, j = 2, . . . , n, (2.3.12)
zeigt, dass die im ersten Eliminationsschritt erzeugte (n − 1) × (n − 1)-Matrix
A(1) =(a
(1)ij
)i,j=2,...,n
symmetrisch ist. Wir wollen zeigen, dass sie auch positiv
definit ist. Dann folgt wieder a(1)22 > 0 und der Eliminationsprozeß kann mit
positivem Pivotelement fortgesetzt werden. Die positive Definitheit zeigen wirdurch Induktion. Dazu sei x = (x2, . . . , xn)T ∈ Rn−1\0 und x = (x1, x)
T ∈ Rn
mit
x1 = − 1
a11
n∑
k=2
a1kxk.
Nun ist
0 <
n∑
j,k=1
ajkxjxk =
n∑
j,k=2
ajkxjxk + 2x1
n∑
k=2
a1kxk + a11x21.
Die Nullerganzung
0 = − 1
a11
n∑
j,k=2
ak1a1jxkxj +1
a11
n∑
k=2
a1kxk
2
fuhrt zu
0 <n∑
j,k=1
ajkxjxk =n∑
j,k=2
ajk −
ak1a1j
a11
xjxk + a11
x1 +
1
a11
n∑
k=2
a1kxk
2
=n∑
j,k=2
a(1)jk xjxk
und damit zu xT A(1)x > 0.
Fur positiv definite Matrizen existiert also stets eine LR-Zerlegung A = LRmit positiven Pivotelementen
rii = a(i)ii > 0, i = 1, . . . , n.
Wegen A = AT gilt aber auch
A = AT = (LR)T = (LDR)T = RTDLT
2.3. SPEZIELLE GLEICHUNGSSYSTEME 33
mit den Matrizen
R =
1 r12/r11 . . . r1n/r11. . .
...1 rn−1,n/rn−1,n−1
1
, D =
r11
. . .
rnn
.
Wegen der Eindeutgkeit der LR-Zerlegung folgt aus
A = LR = RTDLT
notwendig L = RT bzw. R = DLT . Positiv definite Matrizen gestatten also einesog. Cholesky-Zerlegung
A = LDLT = LLT
mit der Matrix L = LD1/2. Bei der Berechnung der Cholesky-Zerlegung genugtes, die Matrizen D und L zu bestimmen. Dies reduziert den Speicherplatzbedarfauf n(n + 1)/2 und die benotigten Operationen auf n3/6 + O(n2).
Der Algorithmus von Cholesky zur Berechnung der Zerlegungsmatrix
L =
l11...
. . .
ln1 . . . lnn
geht direkt von der Beziehung A = LLT aus, die als System von n(n + 1)/2Gleichungen fur die Großen ljk, k ≤ j, auffassen kann. Ausmultiplizieren von
l11...
. . .
ln1 . . . lnn
l11 . . . ln1
. . ....
lnn
=
a11 . . . a1n...
...an1 . . . ann
ergibt in der ersten Spalte von L:
l211 = a11, l21 l11 = a21, . . . , ln1 l11 = an1
woraus sichl11 =
√a11, j = 2, . . . , n : lj1 =
aj1
l11,
berechnen. Seien fur ein i ∈ 2, . . . , n die Elemente ljk, k = 1, . . . , i − 1, j =k, . . . , n, schon bekannt. Dann erhalt man aus
l2i1 + l2i2 + · · ·+ l2ii = aii, lii > 0
lj1li1 + lj2li2 + · · ·+ ljilii = aji
die nachsten Elemente lii und lji, j = i+ 1, . . . , n.
34 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
2.4 Nicht regulare Systeme
Wir betrachten nun allgemeinere Gleichungssysteme der Form Ax = b, bei denenA ∈ Rm×n eine gegebene, nicht notwendig quadratische Matrix und b ∈ Rm dierechte Seite bezeichnen. Wir lassen auch den Fall Rang(A) < Rang[A, b] zu, d..h.das System muß nicht im eigentlichen Sinne losbar sein. Wir verallgemeinern denLosungsbegriff in dem Sinne, dass ein Vektor x ∈ Rn gesucht wird, dessen Defektd ≡ b−Ax die kleinste euklidische Norm besitzt. Im Falle Rang(A) = Rang[A, b]fallt dieser Losungsbegriff mit dem ublichen zusammen.
Theorem 2.4.1 Es existiert eine verallgemeinerte Losung x ∈ Rn von Ax = bmit kleinsten Fehlerquadraten
‖Ax− b‖2 = minx∈Rn‖Ax− b‖2.
Dies ist aquivalent dazu, dass x Losung der Normalgleichung
ATAx = AT b
ist. Im Fall Rang(A) = n ist x eindeutig bestimmt, andernfalls ist jede Losungvon der Form x+ y, mit y ∈ Kern(A).
Beweis. Fur eine Minimallosung x gilt notwendig
0 =∂
∂xi‖Ax− b‖22
∣∣∣x=x
=∂
∂xi
n∑
j=1
∣∣∣∣∣
m∑
k=1
ajkxk − bj∣∣∣∣∣
2∣∣∣x=x
= 2
n∑
j=1
aji
(m∑
k=1
ajkxk − bj)
= 2(ATAx− AT b)i,
somit lost x die Normalgleichung. Sei umgekehrt x Losung der Normalgleichung.Fur beliebiges x ∈ Rn gilt dann
‖b− Ax‖22 = ‖b− Ax+ A(x− x)‖22= ‖b− Ax‖22 + 2(b− Ax, A(x− x)) + ‖A(x− x)‖22.
Das orthogonale Komplement von Bild(A) in Rm ist Kern(AT ). Nun liegt aberb−Ax im Kern von AT , falls x eine Losung der Normalgleichung ist. Folglich
‖b− Ax‖22 = ‖b− Ax‖22 + ‖A(x− x)‖22 ≥ ‖b−Ax‖22.
Bleibt die Losbarkeit des Normalgleichungssystems zu untersuchen. WegenRm = Bild(A)⊕Kern(AT ) kann b eindeutig in
b = s+ r, s ∈ Bild(A), r ∈ Kern(AT )
2.4. NICHT REGULARE SYSTEME 35
zerlegt werden. Fur ein x ∈ Rn mit Ax = s gilt dann
ATAx = AT s = AT s+ AT r = AT b,
d.h. x lost das Normalgleichungssystem. Im Fall Rang(A) = n ist Kern(A) = 0und Bild(A) = R
n. Aus ATAx = 0 folgt wegen Kern(AT ) ⊥ Bild(A) notwendigAx = 0 bzw. x = 0. Die Matrix ATA ∈ Rn×n ist also regular und x eindeutigbestimmt. Im Fall Rang(A) < n gilt fur jede weitere Losung x1 der Normalglei-chung
b = Ax1 + (b− Ax1) ∈ Bild(A) + Kern(AT ).
Aus der Eindeutigkeit dieser orthogonalen Zerlegung schliessen wir Ax1 = s = Axund x− x1 ∈ Kern(A).
Die klassische Anwendung des Satzes ist die sog. Gaußsche Ausgleichsrech-nung. Die Aufgabenstellung besteht dabei in folgendem: Zu gegebenen Funk-tionen ϕ1, . . . , ϕn und Punkten (xj , yj) ∈ R2, j = 1, . . . , m, m > n, ist eineLinearkombination
ϕ(x) =n∑
k=1
ckϕk(x)
derart zu bestimmen, dass die mittlere Abweichung(
m∑
j=1
|ϕ(xj)− yj|2)1/2
moglichst klein wird. Zur Losung dieser Aufgabe setzen wir
y := (y1, . . . , yn)T , c := (c1, . . . , cn)
T ,
ak := (ϕk(x1), . . . , ϕk(xm))T , k = 1, . . . , n, A := [a1, . . . , an]
Dann ist das FunktionalF (c) = ‖Ac− y‖2
bezuglich c ∈ Rn zu minimieren. Dies ist gleichbedeutend damit, fur das uber-bestimmte (m > n) Gleichungssystem Ac = y eine verallgemeinerte Losung mitkleinsten Fehlerquadraten zu bestimmen. Im Fall Rang(A) = n ist diese eindeutigals Losung der Normalgleichung
ATAc = ATy
gegeben. Fur den Spezialfall polynomialer Funktionen ϕk(x) = xk−1 ist die MatrixA gegeben durch
A =
1 x1 . . . xn−11
1 x2 . . . xn−12
......
1 xm . . . xn−1m
, m > n.
36 KAPITEL 2. LINEARE GLEICHUNGSSYSTEME I
Wegen der Regularitat der Vandermondeschen Determinante∣∣∣∣∣∣∣∣∣
1 x1 . . . xn−11
1 x2 . . . xn−12
......
1 xn . . . xn−1n
∣∣∣∣∣∣∣∣∣
=
n∏
j,k=1j<k
(xk − xj) 6= 0
fur paarweise verschiedene Stutzstellen xj ist dann Rang(A) = n und die Gauß-sche Ausgleichsaufgabe eindeutig losbar.
Beispiel. Zu den Meßdaten
xj -2 -1 0 1 2yj 0.5 0.5 2 3.5 3.5
soll mit Hilfe der Gaußschen Ausgleichsrechnung eine lineare Funktion y(x) =a + bx angepasst werden. Dies ist aquivalent zur Losung des uberbestimmtenGleichungssystems
1 −21 −11 01 11 2
[ab
]=
0.50.52
3.53.5
.
Das zugehorige Normalgleichungssystem lautet[5 00 10
] [ab
]=
[109
]
mit der Losung a = 2 und b = 0.9.
Basierend auf den folgenden Satz gelingt die Berechnung von x auch ohneexplizite Aufstellung der Normalengleichung.
Theorem 2.4.2 Sei A ∈ Rm×n eine rechteckige Matrix mitm ≥ n und Rang(A) =n. Dann existiert eine eindeutig bestimmte (orthogonale) Matrix Q ∈ Rm×n mitder Eigenschaft QTQ = I und eine eindeutig bestimmte obere DreiecksmatrixR ∈ Rn×n mit positiven Diagonalelementen rii > 0, i = 1, . . . , n, so dass A = QR.
Wenden wir die QR-Zerlegung auf die Normalgleichung an, so folgt
ATAx = RTQTQRx = RTRx = RTQT b
und wegen der Regularitat vonRT Rx = QT b.Dieses System kann durch Ruckwartseinsetzengelost werden. Wegen
ATA = RTR
ist mit R also die Cholesky-Zerlegung von ATA bestimmt, ohne ATA explizitberechnen zu mussen. Ein stabiles Verfahren zur Berechnung der QR-Zerlegungist das Householder-Verfahren.
Kapitel 3
Interpolation
Ein Grundproblem der Numerik ist die Darstellung und Auswertung von Funk-tionen. Dabei ergeben sich folgende Aufgabenstellungen:
• Eine Funktion f ist nur auf einer diskreten Menge von Argumenten x0, . . . , xn
bekannt und soll mit dieser Information rekonstruiert werden, um z.B. dieFunktion grafisch darzustellen oder um die Funktion an gewissen Zwischen-stellen auszuwerten.
• Eine analytisch gegebene Funktion f soll auf einer Rechenanlage so dar-gestellt werden, dass Funktionswerte f(x) zu beliebigen x leicht innerhalbeiner vorgegebene Toleranz berechnet werden konnen.
In beiden Fallen hat man ein System mit unendlich vielen Freiheitsgraden, namlichdie funktionale Abhangigkeit y = f(x), durch einen endlichen Datensatz zu simu-lieren. Hierzu bedient man sich gewisser Klassen P einfach strukturierter Funk-tionen; z.B.
Polynome: p(x) = c0 + c1x+ c2x2 + · · ·+ cnx
n
Rationale Funktionen: r(x) =c0 + c1x+ c2x
2 + · · ·+ cnxn
b0 + b1x+ b2x2 + · · ·+ bmxm
Trigonometrische Polynome: t(x) =a0
2+
n∑
k=1
ak cos(kx) + bk sin(kx)
Exponentialsummen: e(x) =
n∑
k=1
ak exp(bkx).
Geschieht die Zuordnung eines Elementes g ∈ P zur Funktion f durch Fixierenvon Funktionswerten, etwa
g(xi) = yi := f(xi), i = 0, . . . , n,
37
38 KAPITEL 3. INTERPOLATION
so spricht man von Interpolation. Ist g ∈ P als in einem gewissen Sinn besteDarstellung von f zu bestimmen, etwa als
maxa≤x≤b
|f(x)− g(x)| minimal fur g ∈ P,
so spricht man von Approximation. Die jeweilige Wahl der Konstruktion vong ∈ P hangt von der zu erfullenden Aufgabe ab. In diesem Abschnitt behan-deln wir Interpolationsaufgaben, im Abschnitt 5 wenden wir uns der Frage derApproximation zu.
3.1 Polynominterpolation
Wir bezeichnen mit Pn den Vektorraum der Polynome vom Grad kleiner odergleich n:
Pn := p(x) = c0 + c1x+ · · ·+ cnxn : ci ∈ R, i = 0, . . . , n.
Die Lagrange’sche Interpolationsaufgabe besteht darin, zu n + 1 paarweise ver-schiedenen Stutzstellen (Knoten) x0, x1, . . . , xn ∈ R ein Polynom p ∈ Pn zubestimmen, das in xi vorgegebene Werte yi annimmt, fur das also p(xi) = yi gilt.
Theorem 3.1.1 Die Lagrange’sche Interpolationsaufgabe besitzt eine eindeutigbestimmte Losung.
Beweis: Wir zeigen zunachst die Eindeutigkeit. Angenommen es gabe zwei Lo-sungen p1, p2 ∈ Pn, dann verschwindet das Polynom p := p1−p2 ∈ Pn in den n+1paarweise verschiedenen Punkten x0, . . . xn , und ist folglich Null. Zur Existenzbetrachten wir die Forderungen p(xi) = yi, i = 0, . . . , n, als n + 1 Gleichungenzur Bestimmung der n + 1 unbekannten Koeffizienten ci, i = 0, . . . , n. Da dashomogene System nur die Nulllosung hat (Eindeutigkeit!), hat das inhomogeneSystem stets eine Losung.
Zur expliziten Konstruktion des Interpolationspolynoms p ∈ Pn verwendetman die Lagrange’schen Basispolynome
L(n)i (x) =
n∏
k=0,k 6=i
x− xk
xi − xk
∈ Pn, i = 0, . . . , n.
Eine wichtige Eigenschaft der Lagrange’schen Basispolynome ist
L(n)i (xj) = δij .
Damit kann das gesuchte Polynom in der Form
p =
n∑
i=0
yiL(n)i ∈ Pn
3.1. POLYNOMINTERPOLATION 39
geschrieben werden. Das Polynom heißt Lagrange’sches Interpolationspolynomzu den Stutzpunkten (xi, yi), i = 0, . . . , n.
Tipp: In MATLAB konnen die n+1 Stutzpunkte (xi, yi) in zwei Vektoren x undy gespeichert werden. Dann liefert
c=polyfit(x,y,n)
im Vektor c = (c(1), . . . , c(n + 1)) die Koeffizienten des zugeordneten Interpola-tionspolynoms
p(x) =n+1∑
i=1
c(i)xn+1−i.
Benotigt man die Ableitung des Interpolationspolynoms, so kann man mit
d=polyder(c)
den Koeffizientenvektor der Ableitung bestimmen.
Das Lagrange’sche Interpolationspolynom hat den Nachteil, dass sich die ver-wendeten Basisfunktionen von Pn bei Hinzunahme einer weiteren Stutzstellevollig andern, vgl. Abb. 3.1. Dieser Effekt kann durch Ubergang zu einer an-
0 0.5 1 1.5 2 2.5 3−0.5
0
0.5
1
1.5
0 0.5 1 1.5 2 2.5 3 3.5 4−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Abbildung 3.1: Basispolynome L(3)i auf der Knotenmenge 0, 1, 2, 3 (links) und
L(4)i auf der Knotenmenge 0, 1, 2, 3, 4 (rechts).
deren (hierarchischen) Basis in Pn vermieden werden. Dazu verwendet man dieBasisfunktionen
N0(x) := 1, Ni+1(x) := (x− xi)Ni(x), i = 0, . . . , n− 1.
Die rekursive Definition fuhrt auf
Ni(x) =
i−1∏
j=0
(x− xj), i = 1, . . . , n.
40 KAPITEL 3. INTERPOLATION
Aus dem Ansatz
p(x) =
n∑
i=0
aiNi(x)
folgt dann das gestaffelte Gleichungssystem
y0 = p(x0) = a0
y1 = p(x1) = a0 + a1(x1 − x0)
...
yn = p(xn) = a0 + a1(xn − x0) + · · ·+ an(xn − x0) . . . (xn − xn−1),
aus dem sukzessive die ai ermittelt werden konnen. Insbesondere ist die Hin-zunahme eines weiteren Stutzpunktes (xn+1, yn+1) problemlos moglich, denn diebislang berechneten ai, i = 0, . . . , n, andern sich bei Hinzunahme der neuen Ba-sisfunktion Nn+1 nicht. Praktisch bestimmt man die Koeffizienten ai jedoch aufeine andere Weise, die im folgenden beschrieben wird.
Sei pi,i+k ∈ Pk das Polynom, dass die Stutzpunkte (xi, yi), . . . , (xi+k, yi+k)interpoliert. Zu den Punkten (xi, yi) definieren wir die dividierten Differenzeny[xi, . . . , xi+k] rekursiv durch
i = 0, . . . , n : y[xi] := yi
k = 1, . . . , n− i : y[xi, . . . , xi+k] :=y[xi+1, . . . , xi+k]− y[xi, . . . , xi+k−1]
xi+k − xi
Theorem 3.1.2 Fur i = 0, . . . , n und k = 0, . . . , n− i gilt
pi,i+k(x) = y[xi] + y[xi, xi+1](x− xi) + . . .
· · ·+ y[xi, . . . , xi+k](x− xi) . . . (x− xi+k−1).
Beweis: Induktion bezuglich der Indexdifferenz k = (i + k) − i. Fur k = 0ist pi,i = yi = y[xi], i = 0, . . . , n. Sei die Behauptung nun fur k − 1 richtig. Dapi,i+k ∈ Pk das Polynom bezeichnet, dass die Stutzpunkte (xi, yi), . . . , (xi+k, yi+k)interpoliert, gilt mit gewissem a ∈ R
pi,i+k(x) = pi,i+k−1(x) + a(x− xi) . . . (x− xi+k−1).
Zu zeigen ist also a = y[xi, . . . , xi+k]. Nach Induktionsannahme ist
pi,i+k−1(x) = · · ·+ y[xi, . . . , xi+k−1]xk−1,
pi+1,i+k(x) = · · ·+ y[xi+1, . . . , xi+k]xk−1,
wobei “. . . ” fur Polynomanteile vom Grad kleiner oder gleich k − 2 steht. DasPolynom
q(x) :=(x− xi) pi+1,i+k(x)− (x− xi+k) pi,i+k−1(x)
xi+k − xi
3.1. POLYNOMINTERPOLATION 41
interpoliert die k + 1 Stutzpunkte (xj , yj), j = i, . . . , i+ k. In der Tat haben wirfur j = i+ 1, . . . , i+ k − 1
q(xi) = pi,i+k−1(xi) = yi,
q(xj) =(xj − xi) pi+1,i+k(xj)− (xj − xi+k) pi,i+k−1(xj)
xi+k − xi
=(xj − xi) yj − (xj − xi+k) yj
xi+k − xi
= yj,
q(xi+k) = pi+1,i+k(xi+k) = yi+k.
Aus der Eindeutigkeit des Interpolationspolynoms folgt q = pi,i+k. Fur den Ko-effizienten bei xk in q beziehungsweise pi,i+k gilt somit
a =y[xi+1, . . . , xi+k]− y[xi, . . . , xi+k−1]
xi+k − xi= y[xi, . . . , xi+k],
was zu zeigen war.
Setzt man i = 0 und k = n, so erhalt man die Newton’sche Darstellung desInterpolationspolynoms zu den Stutzpunkten (x0, y0), . . . , (xn, yn)
p(x) =n∑
i=0
y[x0, . . . , xi] Ni(x).
Zur Berechnung der Koeffizienten ai = y[x0, . . . , xi] nutzt man folgendes Schema
xn − x0 . . . x2 − x0 x1 − x0 x0 y0 y[x0, x1] y[x0, x1, x2] . . . y[x0, . . . , xn]x3 − x1 x2 − x1 x1 y1 y[x1, x2] y[x1, x2, x3]
x3 − x2 x2 y2 y[x2, x3]...
...xn yn
Bei Hinzunahme eines weiteren Stutzpunktes (xn+1, yn+1) berechnet man denKoeffizienten y[x0, . . . , xn+1] im Newton’schen Interpolationspolynom einfach durchBerechnung einer weiteren Diagonalen im obigen Schema.
Die im Beweis des Theorem 3.1.2 verwendete Beziehung zwischen den Poly-nomen pi,i+k kann direkt zur rekursiven Berechnung des Interpolationspolynomsp = p0,n verwendet werden. Das durch
pi,i(x) = yi i = 0, . . . , n,
pi,i+k(x) = pi,i+k−1(x) + (x− xi)pi+1,i+k(x)− pi,i+k−1(x)
xi+k − xik = 0, . . . , n− i,
42 KAPITEL 3. INTERPOLATION
erzeugte Polynom p0,n ist die Neville’sche Darstellung des Interpolationspoly-noms zu den Stutzstellen (x0, y0), . . . , (xn, yn) (kurz: Neville’sches Interpolations-polynoms). Die praktische Berechnung erfolgt nach folgendem Schema:
x0 y0 p0,1(x) p0,2(x) p0,3(x) . . . p0,n−1(x) p0,n(x)x1 y1 p1,2(x) p1,3(x) p1,4(x) . . . p1,n(x)x2 y2 p2,3(x) p2,4(x) p2,5(x) . . ....
......
......
xn−1 yn−1 pn−1,n(x)xn yn
Auch hier ist die Hinzunahme eines weiteren Stutzpunktes (xn+1, yn+1) problem-los moglich. Die Neville’sche Darstellung des Interpolationspolynoms bietet einesehr effiziente Moglichkeit zur Berechnung einzelner Funktionswerte p(ξ) (ξ 6= xi)ohne vorherige Bestimmung der Koeffizienten in der Newtonschen Darstellung.Dazu setzt man im obigen Neville-Schema einfach x = ξ und verwendet zurBerechnung von pi,k := pi,k(ξ) die Rekursionsformeln
i = 0, . . . , n pi,i = yi
k = 1, . . . , n− i pi,i+k = pi,i+k−1 +pi,i+k−1 − pi+1,i+k
(ξ − xi+k)/(ξ − xi)− 1.
3.2 Interpolationsfehler
Wir wollen den Interpolationsfehler abschatzen, der bei Ersetzung einer gegebe-nen Funktion f durch ihr Interpolationspolynom pn mit den Knoten x0, x1, . . . , xn
entsteht. Dazu bezeichne [x0, . . . , xn] ⊂ [a, b] das kleinste Intervall, dass alle inder Klammer eingeschlossenen Punkte enthalt.
Theorem 3.2.1 Sei f ∈ Cn+1[a, b]. Dann gibt es zu jedem x ∈ [a, b] ein ξx ∈[x0, . . . , xn, x], so dass
f(x)− pn(x) =f (n+1)(ξx)
(n+ 1)!
n∏
j=0
(x− xj).
Beweis: Fur x ∈ x0, . . . , xn ist die Behauptung aufgrund der Interpolationsei-genschaft des Polynoms pn trivial. Sei nun x ∈ [a, b]\x0, . . . , xn. Wir setzen
l(t) :=n∏
j=0
(t− xj), c(x) :=f(x)− pn(x)
l(x).
Die Funktion t 7→ F (t) := f(t)− pn(t) − c(x)l(t) besitzt in [a, b] mindestens dien + 2 Nullstellen x0, x1, . . . , xn, x. Aus der wiederholten Anwendung des Satzes
3.2. INTERPOLATIONSFEHLER 43
von Rolle folgt, dass die Ableitung t 7→ F (n+1) eine Nullstelle ξx ∈ [x0, . . . , xn, x]hat. Wegen
0 = F (n+1)(ξx) = f (n+1)(ξx)− p(n+1)n (ξx)− c(x)l(n+1)(ξx)
= f (n+1)(ξx)− c(x)(n+ 1)!
folgt die Behauptung des Satzes.
Wir wollen den Fehler bei der Lagrange’schen Interpolation diskutieren. Furgroßes n wird 1/(n+1)! sehr klein. Das Produkt wird klein, wenn die Stutzstellenimmer starker zusammenrucken. Sind also die Ableitungen von f auf [a, b] be-schrankt, so gilt
maxa≤x≤b
|f(x)− pn(x)| → 0 n→∞.
Oft sind jedoch die Ableitung der zu interpolierenden Funktionen nicht beschrankt,z.B.
f(x) =1
1 + x2, |f (n)(x)| ≈ 2n n! O(|x|−(n+2)),
so dass gleichmaßige Konvergenz nicht zu erwarten ist. Der Weierstraß’sche Ap-proximationssatz besagt, dass jede auf [a, b] stetige Funktion beliebig genau durchPolynome approximiert werden kann. Die Vermutung, dass dies mit Lagrange-schen Interpolationspolynomen auf aquidistanten Stutzstellen geschehen kann,ist jedoch im allgemeinen falsch.
Beispiel: Seien f(x) = |x|, x ∈ [−1,+1], Stutzstellen xi = −1 + ih, i =0, . . . , 2m, h = 1/m; x 6∈ −1, 0, 1. In Abb. 3.2 sind die Interpolationspolynomeder Betragsfunktion fur m = 4, m = 8, m = 12 und m = 16 dargestellt. Manerkennt an den Intervallgrenzen deutlich einen Trend zum Uberschwingen. Einewesentlich Verbesserung kann durch Ubergang zu nicht aquidistanten Knotenerreicht werden. Abb. 3.3 zeigt das Interpolationspolynom vom Grade m = 16bei Verwendung der Tschebyscheff-Knoten xi = − cos(πi/16), i = 0, . . . , 16. DieTschebyscheff-Knoten auf dem Intervall [a, b] sind durch
xi =a+ b
2− b− a
2cos
πi
n, i = 0, . . . , n,
gegeben.
Fur die mit den Stutzpunkten (xi, f(xi)) gebildeten dividierten Differenzenschreiben wir f [xi, . . . , xi+k] = y[xi, . . . , xi+k].
Theorem 3.2.2 Sei f ∈ Cn+1[a, b]. Dann gilt fur x ∈ [a, b]\x0, . . . , xn dieDarstellung
f(x)− pn(x) = f [x0, . . . , xn, x]
n∏
j=0
(x− xj)
44 KAPITEL 3. INTERPOLATION
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
1.2
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
1.2
1.4
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−0.5
0
0.5
1
1.5
2
2.5
3
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−2
0
2
4
6
8
10
12
Abbildung 3.2: Interpolationspolynome der Betragsfunktion der Ordnung 4, 8,12 und 16 bei Verwendung aquidistanter Stutzstellen.
und es ist
f [x0, . . . , xn, x] =
∫ 1
0
∫ t1
0
· · ·∫ tn−1
0
∫ tn
0
f (n+1)(x0 + t1(x1 − x0) + · · ·
+ tn(xn − xn−1) + t(x− xn))dtdtn · · · dt2dt1.
Beweis: Wir fuhren den Beweis durch vollstandige Induktion nach der Anzahlder Stutzstellen (in der Reihung x0, . . . , xn). Fur n = 0 gilt
f(x)− p0(x) = f(x)− f(x0) = f [x0, x](x− x0)
= (x− x0)
∫ 1
0
f ′(x0 + t(x− x0)) dt.
Gelte die Behauptung nun fur n − 1 ≥ 0, wir zeigen die Richtigkeit fur n. Ausder Newton’schen Darstellung des Interpolationspolynoms und der Richtigkeit
3.2. INTERPOLATIONSFEHLER 45
−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
1.2
Abbildung 3.3: Interpolationspolynom der Betragsfunktion der Ordnung 16 beiVerwendung der Tschebyscheff-Knoten.
fur n− 1 folgt
f(x)− pn(x) = f(x)−n∑
i=0
f [x0, . . . , xi]
i−1∏
j=0
(x− xj)
= f(x)− pn−1(x)− f [x0, . . . , xn]n−1∏
j=0
(x− xj)
= f [x0, . . . , xn−1, x]n−1∏
j=0
(x− xj)− f [x0, . . . , xn]n−1∏
j=0
(x− xj).
Wegen f [x0, . . . , xn−1, x] = f [x, x0, . . . , xn−1] (Ubungsaufgabe) und der Definitiondividierter Differenzen erhalten wir
f(x)− pn(x) = f [x0, . . . , xn, x]n∏
j=0
(x− xj).
Ferner ist nach Induktionsvoraussetzung (Anderung der Notationen im erstenTerm n 7→ n− 1, t 7→ tn und im zweiten Term n 7→ n− 1, x = xn)
f [x0, . . . , xn−1, x]− f [x0, . . . , xn]
=
∫ 1
0
∫ t1
0
· · ·∫ tn−1
0
f (n)
(x0 + t1(x1 − x0) + · · ·+ tn(x− xn−1)
)
− f (n)(x0 + t1(x1 − x0) + · · ·+ tn(xn − xn−1)
)dtn · · · dt1
=
∫ 1
0
∫ t1
0
· · ·∫ tn
0
d
dtf (n)
(x0 + · · ·+ tn(xn − xn−1) + t(x− xn)
)dtdtn · · · dt1,
46 KAPITEL 3. INTERPOLATION
woraus wegen
d
dtf (n)
(x0+ · · ·+ tn(xn − xn−1) + t(x− xn)
)
= f (n+1)(x0 + · · ·+ tn(xn − xn−1) + t(x− xn)
)(x− xn)
und der Definition dividierter Differenzen die Behauptung folgt.
Die obige Darstellung dividierter Differenzen gestattet fur differenzierbareFunktionen die stetige Fortsetzung fur den Fall, dass einige der Stutzstellen zu-sammenfallen:
f [x0, · · · , xr, xr, · · · , xn] := limε→0
f [x0, · · · , xr, xr + ε, · · · , xn].
Im Extremfall fallen alle Stutzstellen zusammen, wir haben
f [x0, · · · , x0] =
∫ 1
0
∫ t1
0
· · ·∫ tn−1
0
f (n)(x0) dtn · · · dt1 =1
n!f (n)(x0)
und die Newton’sche Interpolationsformel geht in das Taylor-Polynom n-ten Gra-des uber
pn(x) =n∑
i=0
f [x0, · · · , xi]i−1∏
j=0
(x− xj) =n∑
i=0
1
i!f (i)(x0)(x− x0)
i.
3.3 Hermite-Interpolation
Die Langrange’sche Interpolation kann auf den Fall erweitert werden, in demneben Funktionswerten auch Werte der Ableitungen einer Funktion f in gewis-sen (oder allen) Knoten bekannt sind. Die Hermite-Interpolation kann wie folgtcharakterisiert werden:
Gegeben seien paarweise verschiedene Knoten xi, i = 0, . . . , n, derFunktionswert und die Ableitungen y
(k)i bis zur Ordnung k = mi im
Knoten xi, i = 0, . . . , n. Gesucht ist ein Polynom N -ten Grades,
N =
n∑
i=0
(1 +mi)− 1,
mit der Eigenschaft p(k)(xi) = y(k)i , k = 0, . . . , mi, i = 0, . . . , n. Die
Knoten xi werden auch als mi + 1-fache Stutzstellen bezeichnet.
Theorem 3.3.1 Die Hermite’sche Interpolationsaufgabe ist eindeutig losbar.
3.3. HERMITE-INTERPOLATION 47
Beweis: Fur die N + 1 unbekannten Koeffizienten ci des gesuchten Polynomsergeben sich aus den Interpolationsbedingungen N + 1 lineare Bestimmungsglei-chungen. Dieses Gleichungssystem ist genau dann eindeutig losbar, wenn daszugeordnete homogene System nur die Nulllosung besitzt. Sei also p(k)(xi) = 0fur k = 0, . . . , mi, i = 0, . . . , n. Da pN in x0, . . . , xn verschwindet, konnen wir dasgesuchte Polynom in der Form
pN(x) = r0(x)n∏
i=0
(x− xi)
darstellen, wobei r0 ein Polynom vom Grade kleiner oder gleich N − n − 1 ist.Die Ableitung
p′N(x) = r′0(x)n∏
i=0
(x− xi) + r0(x)n∑
j=0
∏
j 6=i
(x− xi)
verschwindet in allen Knoten xi fur die mi ≥ 1, folglich verschwindet dort auchr0. Somit haben wir fur ein gewisses Polynom r1
pN(x) = r1(x)n∏
i=0
2∏
j=1mi≥1
(x− xi), deg r1 ≤ N − n− 1−n∑
i=0mi≥1
1.
Sei m = maxi=0,...,n
mi. Nach m− 1 Schritten haben wir
pN(x) = rm−1(x)n∏
i=0
min(m,mi+1)∏
j=1
(x− xi),
deg rm−1 ≤ N − n− 1−n∑
i=0mi≥1
1− · · · −n∑
i=0mi≥m−1
1 =
n∑
i=0mi=m
1− 1.
Das Polynom rm−1 verschwindet in allen xi, fur die mi = m gilt, d.h. in mehrpaarweise verschiedenen Punkten als sein Grad angibt, kann somit nur das Null-polynom sein.
Wir suchen ahnlich wie bei der Lagrange’schen Interpolation eine Darstellungdes Hermite’schen Interpolationspolynoms in der Form
HN(x) =
n∑
i=0
mi∑
k=0
y(k)i Lik(x). (3.3.1)
Betrachten wir zunachst die durch
lij(x) =(x− xi)
j
j!
n∏
k=0k 6=i
(x− xk
xi − xk
)mk+1
, i = 0, . . . , n, j = 0, . . . , mi,
48 KAPITEL 3. INTERPOLATION
definierten N + 1 Polynome vom Grade kleiner oder gleich N . Die Polynomehaben die Eigenschaft
lij(xi) = l(1)ij (xi) = · · · = l
(j−1)ij (xi) = 0, l
(j)ij (xi) = 1,
lij(xk) = l(1)ij (xk) = · · · = l
(mk)ij (xk) = 0,
Die gesuchten Hermite’schen Basispolynome Lij werden nun rekursiv wie folgtdefiniert:
Limi(x) = limi
(x), i = 0, . . . , n
Lij(x) = lij(x)−mi∑
k=j+1
l(k)ij (xi)Lik(x), j = mi − 1, mi − 2, . . . , 0
Theorem 3.3.2 Die rekursiv definierten Hermite’schen Basispolynome Lij erfullendie Beziehungen
L(r)ij (xk) =
1 falls i = k und j = r0 andernfalls
Hieraus folgt unmittelbar die Darstellung (3.3.1).
Beweis: Ubungsaufgabe.
Theorem 3.3.3 Sei f ∈ CN+1[a, b]. Dann gibt es zu jedem x ∈ [a, b] einen Punktξx ∈ [x0, . . . , xn, x], so dass fur den Fehler des Hermite’schen Interpolationspoly-nom pN folgende Darstellung gilt:
f(x)− pN(x) =1
(N + 1)!f (N+1)(ξx)
n∏
i=0
(x− xi)mi+1
Beweis: Analog zum Fehler des Lagrange’schen Interpolationspolynoms.
Beispiel: Bestimmen Sie das Hermite’sche Interpolationspolynom zu den vorge-gebenen Werten f(0) = 2, f ′(0) = 1, f ′′(0) = 4 und f(1) = −1. Der Ansatz
H3(x) = Ax3 +Bx2 + Cx+D
fuhrt auf das Gleichungssystem
2 = H3(0) = D
1 = H ′3(0) = C
4 = H ′′3 (0) = 2B
−1 = H3(1) = A+B + C +D
3.4. SPLINE-INTERPOLATION 49
mit der Losung A = −6, B = 2, C = 1 und D = 2. Das gesuchte Interpolati-onspolynom ist damit H3(x) = −6x3 + 2x2 + x+ 2. Fur den oben beschriebenenallgemeinen Ansatz folgt
x0 = 0, m0 = 2, x1 = 1, m1 = 0, N = (m0 + 1) + (m1 + 1)− 1 = 3.
Mit den Hilfspolynomen
l00(x) = −(x− 1), l01(x) = −x(x − 1), l02(x) = −x2
2(x− 1), l10(x) = x3
erhalten wir die Hermite’schen Basispolynome Lij rekursiv, zunachst Limifur
i = 0, 1:
L02(x) = l02(x) = −x2
2(x− 1), L10(x) = l10(x) = x3.
Die Rekursionsbeziehung liefert
L01(x) = l01(x)− l01′′(0)L02(x) = −x(x− 1)− (−2)(−x2
2(x− 1)) = x(1− x2)
L00(x) = l00(x)− l00′(0)L01(x)− l00′′(0)L02(x) = 1− x3
Aus der allgemeinen Darstellung
HN(x) =
n∑
i=0
mi∑
k=0
y(k)i Lik(x)
folgt im vorliegenden Fall
H3(x) = 2(1− x3) + x(1− x2) + 2x2(1− x)− x3 = −6x3 + 2x2 + x+ 2.
Der Vorteil der Bestimmung der Hermite’schen Basispolynome besteht darin, dasswir die Hermite-Interpolation fur beliebige Daten f(0), f ′(0), f ′′(0), f(1) gelosthaben, es gilt namlich
H3(x) = f(0)(1− x3) + f ′(0)x(1− x2) +f ′′(0)
2x2(1− x) + f(1)x3.
3.4 Spline-Interpolation
Wie wir gesehen haben, eignen sich Lagrangesche Interpolationspolynome nichtbesonders gut zur Approximation von (nicht glatten) Funktionen, da sie bei
50 KAPITEL 3. INTERPOLATION
Vermehrung der Stutzstellenzahl dazu neigen, zwischen den Stutzstellen immergroßere Werte anzunehmen. Dies ist die Folge ihrer Steifheit bedingt durch dieForderung von C∞-Ubergangen in den Knoten. Zur Reduzierung dieser Steifheitsetzt man die interpolierende Funktion stuckweise polynomial bezuglich einerZerlegung a = x0 < x1 < · · · < xn = b an. In den Knoten werden dann geeigne-te Differenzierbarkeitseigenschaften vorausgesetzt. Wir bezeichnen die Lange desTeilintervalls Ii = [xi−1, xi] durch hi = xi − xi−1, die Große h = maxi=1,...,n hi
charakterisiert die Feinheit der Intervallzerlegung. Auf einer solchen Intervallzer-legung werden Vektorraume von stuckweise polynomialen Funktionen betrachtet
S(k,r)h [a, b] = p ∈ Cr[a, b] : p|Ii
∈ Pk(Ii), i = 1, . . . , n, k, r = 0, 1, 2, . . .
Zu einem Satz von Stutzwerten in Punkten aus dem Intervall [a, b] wird dann ei-
ne Interpolierende p ∈ S(k,r)h [a, b] mit Hilfe geeigneter Interpolationsbedingungen
bestimmt. Wir betrachten nun einige Beispiele.
Beispiel: Die stetige, stuckweise lineare Lagrange-Interpolation (Fall k = 1,r = 0) approximiert eine gegebene Funktion f auf [a, b] durch einen Polygonzugin den Stutzstellen xi, i = 0, . . . , n:
p ∈ S(1,0)h [a, b] = p ∈ C[a, b], p|Ii
linear , p(xi) = f(xi), i = 0, . . . , n.
Die Anwendung der Fehlerabschatzung fur die Lagrange-Interpolation separatauf jedem der Teilintervalle Ii
f(x)− p(x) =f ′′(ξx)
2(x− xi−1)(x− xi)
|f(x)− p(x)| ≤ 1
2maxx∈Ii
|f ′′(x)| h2i
4x ∈ Ii
ergibt die globale Fehlerabschatzung
maxx∈[a,b]
|f(x)− p(x)| ≤ h2
8maxx∈[a,b]
|f ′′(x)|.
Fur die Konstruktion der Interpolierenden verwendet man die KnotenfunktionaleNi(f) = f(xi), i = 0, . . . , n, die eine Knotenbasis ϕj ∈ S
(1,0)h [a, b], j = 0, . . . , n
eindeutig durch die Bedingung
Ni(ϕj) = δij, i, j = 0, . . . , n
festlegen. Die Interpolierende kann dann in der Form
p(x) =
n∑
i=0
Ni(f)ϕi(x)
3.4. SPLINE-INTERPOLATION 51
xxj−1 xj xj+1
ϕj
Abbildung 3.4: Lokale Basisfunktion ϕj zum Knotenfunktional Nj .
dargestellt werden.
Beispiel: Stetige, stuckweise kubische Lagrange-Interpolation (Fall k = 3, r =0). Zur Erzielung globaler Stetigkeit verwendet man die Knotenfunktionale
Ni(f) = f(xi), i = 0, . . . , n.
Diese werden in jedem Intervall Ii durch zwei weitere
Nij(f) = f(xij), xij ∈ (xi−1, xi), j = 0, 1
erganzt. Damit ist eindeutig eine global stetige Interpolierende p ∈ S(3,0)h [a, b]
festgelegt. Man erhalt die Fehlerabschatzung
maxx∈[a,b]
|f(x)− p(x)| ≤ h4
4!maxx∈[a,b]
|f (4)(x)|
intervallweise aus der Abschatzung fur kubische Lagrange-Interpolation.
Beispiel: Stetig differenzierbare, stuckweise kubische Hermite-Interpolation (Fallk = 3, r = 1). Zur Erzielung globaler stetiger Differenzierbarkeit verwendet mandie Knotenfunktionale
Ni(f) = f(xi), Nn+i+1(f) = f ′(xi), i = 0, . . . , n.
Damit ist eindeutig eine global stetig differenzierbare Interpolierende p ∈ S(3,1)h [a, b]
festgelegt. Man erhalt fur den Fehler
maxx∈[a,b]
|f(x)− p(x)| ≤ h4
4!maxx∈[a,b]
|f (4)(x)|,
diesmal aus der intervallweisen Anwendung der Abschatzung fur kubische Hermite-Interpolation.
Die Forderung nach hoherer Glattheit, etwa eine Interpolation in S(k,k−1)h [a, b]
fuhrt auf die Spline-Interpolation, die von großer praktischer Bedeutung, etwa in
52 KAPITEL 3. INTERPOLATION
der Computer-Graphik ist.
Unsere Aufgabe besteht nun darin, aus vorgegebenen Werten
(xi, yi), i = 0, . . . , n,
eine global zweimal stetige, stuckweise kubische Interpolation sn zu bestimmen.Ein solche Funktion wird interpolierender kubischer Spline genannt.
Theorem 3.4.1 Der interpolierende kubische Spline sn existiert und ist eindeu-tig bestimmt durch zusatzliche Vorgabe von einer der folgenden Randbedingungen
(a) s′′n(a) = s′′n(b) = 0 (naturlicher kubischer Spline)
(b) s′n(a) = s′n(b) und s′′n(a) = s′′n(b) (periodischer kubischer Spline)
(c) s′n(a) = y′(a) und s′n(b) = y′(b) (gebundener kubischer Spline)
Beweis. Jeder kubische Spline hat bezuglich der Zerlegung a = x0 < x1 <· · · < xn = b die Form s|Ii
= pi, pi ∈ P3(Ii), i = 1, . . . , n. Jedes der kubischenPolynome pi hat 4 unbestimmte Koeffizienten, dies ergibt 4n Freiheitsgrade. Zuihrer Bestimmung stehen folgende lineare Beziehungen zur Verfugung:
s(xi) = yi, i = 0, . . . , n : 2n Gleichungens′ ∈ C[a, b] : n-1 Gleichungens′′ ∈ C[a, b] : n-1 GleichungenZusatzbedingungen : 2 GleichungenInsgesamt : 4n Gleichungen
Zum Nachweis der Existenz einer Losung des linearen Gleichungssystems von4n Gleichungen mit 4n Unbekannten genugt es, wie ublich zu zeigen, dass daszugeordnete homogene System nur die Nulllosung besitzt.
Sei f ∈ C2[a, b] eine beliebige, die Daten (xi, yi), i = 0, . . . , n, interpolierendeFunktion mit
s′′n(x) (f ′(x)− s′n(x))∣∣bx=a
= 0.
Dann erhalten wir durch elementweise partielle Integration
∫ b
a
s′′n(x) (f ′′(x)− s′′n(x)) dx =n∑
i=1
∫ xi
xi−1
s′′n(x) (f ′′(x)− s′′n(x)) dx
=
n∑
i=1
s′′n(x)(f ′(x)− s′n(x))
∣∣xi
xi−1− s′′′n (x)(f(x)− sn(x))
∣∣xi
xi−1
+
∫ xi
xi−1
s(iv)(x)(f ′(x)− s′n(x)) dx.
3.4. SPLINE-INTERPOLATION 53
Nun sind s(iv)(x) = 0 fur x ∈ Ii, i = 1, . . . , n, und f(xi) − sn(xi) = 0 furi = 0, . . . , n. Die verbleibene Summe reduziert sich wegen der Stetigkeit von s′′nund f ′ − s′n auf die Differenz der Werte an den Intervallenden, es folgt
∫ b
a
s′′n(x) (f ′′(x)− s′′n(x)) dx = s′′n(b) (f ′(b)− s′n(b))− s′′n(a) (f ′(a)− s′n(a)) = 0.
Wir betrachten nun den Fall eines Null interpolierenden Splines, d.h. dasssn(xi) = 0 fur i = 0, . . . , n gilt. Fur den gebundenen Spline sei zusatzlichs′n(b) = s′n(a) = 0. Dann genugt die Nullfunktion f = 0 allen oben getroffe-nen Voraussetzungen, unabhangig davon, ob ein naturlicher, periodischer odergebundener Spline vorliegt. Wir bekommen
∫ b
a
|s′′n(x)|2 dx = 0,
somit sn linear auf [a, b]. Wegen sn(a) = sn(b) = 0 folgt sn(x) ≡ 0.
Die im Beweis beobachtete Orthogonalitatseigenschaft hat die interessanteKonsequenz, dass sich der interpolierende kubische Spline durch besonders gerin-ge Oszillation auszeichnet.
Theorem 3.4.2 Unter allen Funktionen f ∈ C2[a, b], die den Interpolationsbe-dingungen f(xi) = yi, i = 0, . . . , n, und einer der Randbedingungen
(a) f ′′(a) = f ′′(b) = 0
(b) f ′(a) = f ′(b) und f ′′(a) = f ′′(b)
(c) f ′(a) = y′(a) und f ′(b) = y′(b)
genugen, besitzt die interpolierende kubische Splinefunktion sn die kleinste Ge-samtkrummung in folgendem Sinne:
∫ b
a
|s′′n(x)|2 dx ≤∫ b
a
|f ′′n(x)|2 dx.
Beweis. Aus der Identitat f = sn + (f − sn) folgt mit der oben beobachtetenOrthogonalitat der zweiten Ableitungen
∫ b
a
|f ′′(x)|2 dx−∫ b
a
|s′′n(x)|2 dx
= 2
∫ b
a
s′′n(x)(f ′′(x)− s′′n(x)) dx+
∫ b
a
|f ′′(x)− s′′n(x)|2 dx ≥ 0
die Behauptung.
54 KAPITEL 3. INTERPOLATION
Bemerkung. Der Name “Spline” erklart sich durch die physikalische Inter-pretation der obigen Minimaleigenschaft. Beschreibt y = f(x) die Lage einerdunnen Holzlatte, so mißt
E =
∫ b
a
(y′′(x)
(1 + y′(x)2)3/2
)2
dx
die “Biegeenergie” der Latte. Aufgrund des Hamiltonschen Prinzips stellt sichdie Latte so ein, dass diese Energie minimiert wird. Fur kleine Auslenkungen giltnaherungsweise
E ≈∫ b
a
y′′(x)2 dx,
der interpolierende kubische Spline beschreibt also annahernd die Lage einerdunnen Holzlatte, die an den Knoten xi fixiert ist. Bei der gebundenen Spline-Interpolation haben wir die Latte an den Randknoten unter zusatzlicher Vorgabeder Steigungen an den Randknoten eingespannt. Die naturlichen Randbedingun-gen entsprechen der Situation, dass die Latte ausserhalb des Intervalls [a, b] ge-rade ist. Bei den periodischen Randbedingungen handelt es sich um ringformiggeschlossenen Latten. Derartige dunne Holzlatten wurde tatsachlich als Zeichen-werkzeug benutzt und tragen im Englischen den Namen “Spline”.
Zur expliziten Bestimmung des interpolierenden kubischen Spline sn benutzenwir die Bezeichnungen
yi = sn(xi), Mi = s′′n(xi), i = 0, . . . , n.
Da sn|Ii= pi ∈ P3(Ii), ist p′′i auf Ii linear, somit
p′′i (x) = Mi−1xi − xhi
+Mix− xi−1
hix ∈ Ii, i = 1, . . . , n.
Zweimalige Integration ergibt
pi(x) = Mi−1(xi − x)3
6hi+Mi
(x− xi−1)3
6hi+ Ai−1(x− xi−1) +Bi−1
mit noch zu bestimmenden Konstanten Ai−1 und Bi−1, i = 1, . . . , n. Die Inter-polationsbedingungen pi(xi−1) = sn(xi−1) = yi−1 und pi(xi) = sn(xi) = yi liefern
yi−1 = Mi−1h2
i
6+Bi−1, yi = Mi
h2i
6+ Ai−1hi +Bi−1
woraus
Bi−1 = yi−1 −Mi−1h2
i
6, Ai−1 =
yi − yi−1
hi− hi
6(Mi −Mi−1)
3.4. SPLINE-INTERPOLATION 55
folgt. Die Forderung der Stetigkeit der ersten Ableitungen in xi, i = 1, . . . , n− 1,fuhrt auf
Mihi
2+ Ai−1 = p′i(xi) = p′i+1(xi) = −Mi
hi+1
2+ Ai
woraus sich ein lineares Gleichungssystem von n−1 Gleichungen, i = 1, . . . , n−1:
hi
hi + hi+1
Mi−1 + 2Mi +hi+1
hi + hi+1
Mi+1 =6
hi + hi+1
(yi+1 − yi
hi+1
− yi − yi−1
hi
)
fur die n+ 1 gesuchten Großen Mi, i = 0, . . . , n, ergibt.
Im Fall der naturlichen Randbedingungen komplettieren die beiden Zusatzbe-dingungen M0 = s′′n(a) = Mn = s′′n(b) = 0 das Gleichungssystem, das dann nachElimination von M0 und Mn die Form
2h2
h1 + h20 . . . 0
h2
h2 + h3
2h3
h2 + h3
0
. . .. . .
. . .
0hn−2
hn−2 + hn−12
hn−1
hn−2 + hn−1
0 . . . 0hn−1
hn−1 + hn2
M1
M2
...
Mn−2
Mn−1
=
d1
d2
...
dn−2
dn−1
annimmt. Hierbei wurde zur Abkurzung
di =6
hi + hi+1
(yi+1 − yi
hi+1− yi − yi−1
hi
), i = 1, . . . , n− 1,
benutzt.
Im Fall gebundener Splines erganzt man die Gleichungen durch die Forderun-gen
y′(a) = p′1(x0) = −M0h1
2+ A0 = −M0
2h1
6−M1
h1
6+y1 − y0
h1
y′(b) = p′n(xn) = Mnhn
2+ An−1 = Mn
2hn
6+Mn−1
hn
6+yn − yn−1
hn
56 KAPITEL 3. INTERPOLATION
und man erhalt das System
2 1 0 . . . 0
h1
h1 + h22
h2
h1 + h2. . . 0
. . .. . .
. . .
0 . . .hn−1
hn−1 + hn2
hn
hn−1 + hn
0 . . . 0 1 2
M0
M1
...
Mn−1
Mn
=
d0
d1
...
dn−1
dn
mit
d0 =6
h1
(y1 − y0
h1
− y′(a))
dn =6
hn
(y′(b)− yn − yn−1
hn
).
Beide Gleichungssysteme (naturliche und gebundene Splines) sind strikt dia-gonal dominant (keine Pivotierung erforderlich) und konnen effizient mit demThomas-Verfahren (Tridiagonalmatrix) gelost werden.
Interpolierende Spline-Funktionen besitzen bessere Approximationseigenschaf-ten fur
h = maxi=1,...,n
hi → 0
als die Lagrangeschen Interpolationspolynome.
Theorem 3.4.3 Sei f ∈ C4[a, b]. Fur den interpolierenden kubischen Spline mits′′n(a) = f ′′(a) = s′′n(b) = f ′′(b) = 0 gilt
maxa≤x≤b
|f (j)(x)− s(j)n (x)| ≤ C h4−j max
a≤x≤b|f (4)(x)|, j = 0, 1, 2.
Beweis. Siehe: H. Werner, R. Schaback, Praktische Mathematik II, Springer-Verlag, 1972
Theorem 3.4.4 Sei f ∈ C4[a, b]. Fur den gebundenen interpolierenden kubi-schen Spline gilt
maxa≤x≤b
|f(x)− sn(x)| ≤ 5
384h4 max
a≤x≤b|f (4)(x)|.
Beweis. Siehe: Hall und Meyer, Optimal error bounds for cubic spline interpo-lation, J. Approx. Theory 16(1976), 105-122
Tipp: In MATLAB kann der gebundene, interpolierende kubische Spline mit Hil-fe des Komandos yy=spline(x,y,xx) berechnet werden. Dazu enthalt der Vektor x
3.4. SPLINE-INTERPOLATION 57
die Knotenwerte xi, der Vektor y die zu interpolierenden Funktionswerte yi sowieals erstes und letztes Element die Steigungen y′(a), y′(b), und xx die zu berechnen-den Abzissen. Die zum Abzissenvektor gehorigen Werte der Splineapproximationstehen dann im Ausgabevektor yy. In dieser Variante des spline-Befehls gilt alsosize(y) = size(x) + 2.Sind die Langen der Vektoren x und y gleich, so werden die zwei erforderlichenZusatzbedingungen fur die Eindeutigkeit des interpolierenden kubischen Splinesaus der Forderung bestimmt, dass die dritten Ableitungen von sn in den Knotenx1 und xn−1 stetig sind (“not a knot” Bedingungen oder “Einheitlichkeitsbedin-gungen”). Das beinhaltet, dass sn auf den Intervallen [x0, x2] und [xn−2, xn] nichtnur stuckweise kubisch sondern kubisch ist. Man sagt auch, die Knoten x1 undxn−1 sind keine aktiven Knoten. Die aus den Einheitlichkeitsbedingungen resul-tierenden Gleichungen sind
M1 −M0
h1=M2 −M1
h2
Mn −Mn−1
hn=Mn−1 −Mn−2
hn−1.
Durch diese zusatzlichen Gleichungen ist die Koeffizientenmatrix des Gleichungs-systems zur Bestimmung der M0, . . . , Mn, nicht mehr strikt diagonal dominant,eliminiert man jedoch M0 und Mn aus diesem System, erhalt man wieder ei-ne strikt diagonal dominante Matrix und die eindeutige Losbarkeit ist gesichert.Beispielsweise fuhrt die Elimination von
M0 =
(h1 + h2
h2
)M1 −
h1
h2M2
in der Gleichungh1
h1 + h2M0 + 2M1 +
h2
h1 + h2M2 = d1
auf (2 +
h1
h2
)M1 +
(1− h1
h2
)M2 = d1
und die strikte Diagonaldominanz ist erfullt.
Beispiel: In Abbildung 3.5 sind der kubische interpolierende Spline der Funk-tion y = 1/(1 + x2) auf dem Intervall [−5, 5] sowie das Lagrange’sche Interpo-lationspolynom vom Grade 12 dargestellt. Hierzu wurden die MATLAB-Befehlecs=spline(x,y,xx) sowie p=polyfit(x,y,12) und interpol=polyval(p,xx) verwendet.In den Vektoren x und y sind die 13 aquivalenten Stutzstellen xi = −5+5(i−1)/6und die zugeordneten Funktionswerte y(xi), i = 1, . . . , 13, gespeichert. Deutlichist die schlechte Approximation der Runge-Funktion durch das Interpolations-polynoms am Rande zu sehen. Die Abschwachung der C∞ Glattheit im Innerndes Intervalls [−5, 5] und die Hinzunahme zweier Zusatzbedingungen am Rande(Einheitlichkeitsbedingungen) fuhren zu drastischen Verbesserungen der Appro-ximationsgute.
58 KAPITEL 3. INTERPOLATION
−5 −4 −3 −2 −1 0 1 2 3 4 5−4
−3.5
−3
−2.5
−2
−1.5
−1
−0.5
0
0.5
1
Abbildung 3.5: Stuckweise kubische Spline-Approximation im Vergleich zum La-grangeschen Interpolationspolynom vom Grade 12.
Kapitel 4
Numerische Integration
Die Berechnung bestimmter Integrale kann in der Praxis meist nur naherungsweisemit Hilfe von “Quadraturformeln” erfolgen. Dazu macht man fur eine Funktionf ∈ C[a, b] den Ansatz
I(f) =
∫ b
a
f(x) dx ≈ In(f) = (b− a)n∑
i=0
ωif(xi)
mit Stutzstellen a ≤ x0 < x1 < · · · < xn ≤ b und Gewichten ωi ∈ R. Ein typischesBeispiel ist die Trapezregel, bei der n = 1, x0 = a, x1 = b, ω0 = ω1 = 1/2 giltund zu ∫ b
a
f(x) dx ≈ b− a2
(f(a) + f(b))
fuhrt.
4.1 Beispiele interpolatorischer Quadraturformeln
Ein naheliegender Weg zur Konstruktion von Quadraturformeln ist der uber diePolynominterpolation. Zu paarweise verschiedenen Stutzstellen a ≤ x0 < x1 <· · · < xn ≤ b wird das Lagrangesche Interpolationspolynom gebildet
pn(x) =n∑
i=0
f(xi)L(n)i (x)
und dann die Integration uber f durch die uber pn ersetzt, wir setzen also
In(f) :=
∫ b
a
pn(x) dx = (b− a)n∑
i=0
f(xi)1
b− a
∫ b
a
L(n)i (x) dx =
n∑
i=0
ωif(xi).
Die Gewichte ωi hangen offenbar nur von den Stutzstellen x0, . . . , xn, insbesonderenicht von f ab. Der Quadraturfehler einer interpolatorischen Quadraturformellaßt sich leicht angeben:
59
60 KAPITEL 4. NUMERISCHE INTEGRATION
Theorem 4.1.1 Fur interpolatorische Quadraturformeln gilt:
I(f)− In(f) =
∫ b
a
f [x0, x1, . . . , xn, x]n∏
j=0
(x− xj)dx.
Beweis. Folgt sofort aus der Fehlerdarstellung der Lagrange-Interpolation.
Aus der Fehlerdarstellung folgt, dass die interpolatorische QuadraturformelIn exakt fur Polynome vom Grade kleiner oder gleich n ist.
Mittelpunkts- oder Rechteckregel
Wir ersetzen f im Intervall [a, b] durch den Funktionswert im Mittelpunkt desIntervalls und erhalten
In(f) = (b− a)f(a + b
2
)
mit dem Gewicht ω0 = 1 und der Stutzstelle x0 = (a + b)/2. Fur f ∈ C2[a, b]genugt der Quadraturfehler der Beziehung
E0(f) =
∫ b
a
f(x) dx− (b− a)f(a + b
2
)=
(b− a)3
24f ′′(ξ), ξ ∈ (a, b).
Dies folgt aus der Taylorentwicklung
f(x) = f(x0) + f ′(x0)(x− x0) + f ′′(x0 + θ(x− x0))(x− x0)
2
2
durch Integration und Anwendung des Zwischenwertsatzes der Integralrechnung
∫ b
a
f ′′(x0 + θ(x− x0))(x− x0)
2
2dx = f ′′(ξ)
∫ b
a
(x− x0)2
2dx =
(b− a)3
24f ′′(ξ).
Die Mittelpunktsformel ist nach Konstruktion exakt fur Polynome vom GradeNull, die Fehlerabschatzung zeigt, dass sie fur Polynome vom Grade kleiner odergleich 1 exakt ist. Zur Erhohung der Genauigkeit, kann man das Intervall [a, b] inTeilintervalle der Breite H = (b − a)/m zerlegen und wendet die Mittelpunkts-formel auf jedem der Teilintervalle an. Die Mittelpunkte der Teilintervalle sindgegeben durch xk = a+(2k+1)H/2, k = 0, 1, . . . , m−1. Im Ergebnis erhalt mandie zusammengesetzte Mittelpunktsformel
I0,m(f) = H
m−1∑
k=0
f(xk), m ≥ 1.
4.1. BEISPIELE INTERPOLATORISCHER QUADRATURFORMELN 61
Sei wieder f ∈ C2[a, b]. Dann gilt fur den Quadraturfehler der zusammengesetztenMittelpunktsformel
E0,m(f) := I(f)− I0,m(f) =H3
24
m−1∑
k=0
f ′′(ξk)
=b− a24
H2 1
m
m−1∑
k=0
f ′′(ξk) =b− a24
H2f ′′(ξ).
Trapezregel
Ersetzt man f durch das Lagrangesche Interpolationspolynom bezuglich der Kno-ten x0 = a und x1 = b, so erhalt man die Trapezregel
I1(f) =b− a
2(f(a) + f(b))
mit den Gewichten ω0 = ω1 = 1/2 und den Stutzstellen x0 = a und x1 = b. Istf ∈ C2[a, b], so ist der Quadraturfehler gegeben durch
E1(f) = I(f)− I1(f) = −(b− a)3
12f ′′(ξ), ξ ∈ (a, b).
Tatsachlich erhalt man aus der Formel fur den Interpolationsfehler und dem Mit-telwertsatz der Integralrechnung
E1(f) =1
2
∫ b
a
f ′′(ξ(x))(x− a)(x− b) dx
= −f′′(ξ)
2
∫ b
a
(x− a)(b− x) dx = −(b− a)3
12f ′′(ξ).
Zur Erhohung der Genauigkeit, kann man das Intervall [a, b] in Teilintervalle derBreiteH = (b−a)/m zerlegen und wendet die Trapezregel auf jedem der Teilinter-valle getrennt an. Dazu seien xk = a+ kH , k = 0, 1, . . . , m, die Quadraturknotenund wir erhalten
I1,m(f) =H
2
m−1∑
k=0
(f(xk) + f(xk+1))
= H
(1
2f(x0) + f(x1) + · · ·+ f(xm−1) +
1
2f(xm)
).
Cavalieri-Simpson-Formel
Ersetzt man f durch das Lagrangesche Interpolationspolynom bezuglich der Kno-ten x0 = a, x1 = (a+b)/2 und x2 = b, so erhalt man die Cavalieri-Simpson-Formel
I2(f) =b− a
6
(f(a) + 4f
(a+ b
2
)+ f(b)
)
62 KAPITEL 4. NUMERISCHE INTEGRATION
mit den Gewichten ω0 = ω2 = 1/6 und ω1 = 2/3 sowie den Stutzstellen x0 = a,x1 = (a+ b)/2 und x2 = b. Im vorliegenden Fall gilt fur den Quadraturfehler
E2(f) = I(f)− I2(f) =1
3!
∫ b
a
f ′′′(ξ(x))(x− x0)(x− x1)(x− x2) dx,
aber der Faktor (x− x1) andert sein Vorzeichen im Integrationsintervall, so dasswir den Mittelwertsatz der Integralrechnung nicht direkt anwenden konnen. Wirschreiben unter Verwendung der Newtonschen Darstellung des Interpolationsfeh-lers den Quadraturfehler wie folgt um
E2(f) =
∫ b
a
f [x0, x1, x2, x](x− x0)(x− x1)(x− x2) dx
=
∫ b
a
f [x0, x1, x2, x]− f [x0, x1, x2, x1]
x− x1(x− x0)(x− x1)
2(x− x2) dx
+ f [x0, x1, x2, x1]
∫ b
a
(x− x0)(x− x1)(x− x2) dx.
Das im zweiten Summanden stehende Integral
Q =
∫ b
a
(x− x0)(x− x1)(x− x2) dx
verschwindet, denn die Transformation x = x0 +x2− t ergibt Q = −Q. Der unterdem Integralzeichen stehende Term (x−x0)(x−x1)
2(x−x2) wechselt in [a, b] dasVorzeichen nicht; Anwendung des Mittelwertsatzes der Integralrechnung und dieDarstellbarkeit der dividierten Differenzen in Integralform ergibt nun
E2(f) =
∫ b
a
f [x0, x1, x2, x, x1](x− x0)(x− x1)2(x− x2) dx
=
∫ b
a
f (4)
4!(ξ(x))(x− x0)(x− x1)
2(x− x2) dx
=f (4)
4!(ξ)
∫ b
a
(x− x0)(x− x1)2(x− x2) dx.
Ist f ∈ C4[a, b], so ist der Quadraturfehler gegeben durch
E2(f) = −(b− a)5
2880f (4)(ξ).
Man beachte, dass nach Konstruktion Polynome vom Grade kleiner oder gleich 2mit der Cavalieri-Simpson-Formel exakt integriert werden. Die Fehlerabschatzungzeigt nun, dass sogar Polynome vom Grade kleiner oder gleich 3 exakt integriertwerden.
4.2. NEWTON-COTES-FORMELN 63
Zur Erhohung der Genauigkeit, kann man das Intervall [a, b] wieder in Teil-intervalle der Breite H = (b− a)/m zerlegen und wendet die Cavalieri-Simpson-Formel auf jedem der Teilintervalle getrennt an. Dazu seien xk = a + kH/2,k = 0, 1, . . . , 2m, die Quadraturknoten und wir erhalten
I2,m(f) =H
6
(f(x0) + 2
m−1∑
r=1
f(x2r) + 4m−1∑
s=0
f(x2s+1) + f(x2m)
).
Fur den Quadraturfehler der zusammengesetzten Cavalieri-Simpson-Formel folgt
E2,m(f) = −b− a5760
H4f (4)(ξ), ξ ∈ (a, b).
4.2 Newton-Cotes-Formeln
Die Newton-Cotes-Formeln sind interpolatorische Quadraturformeln, die auf eineaquidistante Verteilung der Stutzstellen basieren. Man unterscheidet
(i) abgeschlossene Newton-Cotes-Formeln (a, b sind Stutzstellen)
xi = a + ih, i = 0, . . . , n, h =b− an
,
(ii) offene Newton-Cotes-Formeln (a, b sind keine Stutzstellen)
xi = a+ (i+ 1)h, i = 0, . . . , n, h =b− an+ 2
,
In beiden Fallen wird der Integrand f durch das entsprechende Interpolations-polynom vom Grade n ersetzt, also
In(f) =
n∑
i=0
f(xi)
∫ b
a
L(n)i (x) dx.
Fur die Gewichte der abgeschlossenen Formeln erhalt man nach Koordinaten-transformation x = a+ th
ωi =1
b− a
∫ b
a
L(n)i (x) dx =
1
n
∫ n
0
L(n)i (a + th) dt =
1
n
∫ n
0
n∏
j=0j 6=i
(t− ji− j
)dt
Diese Gewichte werden ein fur allemal berechnet und tabelliert. Im Fall der of-fenen Newton-Cotes-Formeln verfahrt man analog. Man beachte, dass bei denabgeschlossenen Newton-Cotes-Formeln ab n = 8 und bei den offenen ab n = 2negative Gewichte auftreten. Die fur das Integral geltende Positivitatseigenschaft
f(x) ≥ 0 fur x ∈ [a, b] ⇒ I(f) ≥ 0
kann im Fall negativer Gewichte fur die Quadraturformel In nicht mehr gesi-chert werden. Zusatzlich erhoht sich die Rundungsfehleranfalligkeit der Formeln(Ausloschungsgefahr). Um die Genauigkeit zu erhohen, ist es daher ratsam zuzusammengesetzten Formeln und nicht zu hoheren Werten von n uberzugehen.
64 KAPITEL 4. NUMERISCHE INTEGRATION
Tabelle 4.1: Gewichte abgeschlossener Newton-Cotes Formeln
n ω0 . . . , ωn Name
1 12
12
Trapezregel
2 16
46
16
Cavalieri-Simpson
3 18
38
38
18
Newtonsche 3/8-Regel
4 790
3290
1290
3290
790
Milne-Regel
5 19288
75288
50288
50288
75288
19288
6 41840
216840
27840
272840
27840
216840
41840
7 75117280
357717280
132317280
298917280
298917280
132317280
357717280
75117280
8 98928350
588828350
− 92828350
1049628350
− 454028350
1049628350
− 92828350
588828350
98928350
Tabelle 4.2: Gewichte offener Newton-Cotes Formeln
n ω0 . . . , ωn Name
0 1 Mittelpunktsregel
1 12
12
2 23−1
323
3 1124
124
124
1124
4 1120−14
202620−14
201120
4.3. GAUSSSCHE QUADRATURFORMELN 65
4.3 Gaußsche Quadraturformeln
Die interpolatorischen Quadraturformeln zu den Stutzstellen x0, x1, . . . , xn sindnach Konstruktion mindestens exakt fur Polynome vom Grade kleiner oder gleichn. Wir haben gesehen, dass die Rechteckregel (n = 0) Polynome vom Grade klei-ner oder gleich 1, die Cavalieri-Simpson-Formel (n = 2) Polynome vom Gradekleiner oder gleich 3 exakt integriert. Es stellt sich die Frage, die Stutzstellenx0, x1, . . . , xn und die Gewichte ω0, ω1, . . . , ωn so zu wahlen, dass Polynome mog-lichst hohen Grades exakt integriert werden.
Eine obere Schranke fur den maximalen Grad der Polynome, die von derQuadraturformel
In(f) = (b− a)n∑
i=0
ωif(xi)
exakt integriert werden, ergibt sich aus der Uberlegung, dass fur das Polynomvom Grade 2n+ 2
p(x) =n∏
i=0
(x− xi)2
In(p) verschwindet, jedoch I(p) nicht.
Wir wollen im folgenden zeigen, dass es tatsachlich interpolatorische Quadra-turformeln zu n+1 Stutzstellen gibt, die Polynome vom Grade kleiner oder gleich2n+ 1 exakt integrieren. Sie heißen Gaußsche Quadraturformeln.
Seien pn ∈ Pn und p2n+1 ∈ P2n+1 die Lagrangeschen Interpolationspolynomeeiner Funktion f ∈ C[a, b] zu den n + 1 bzw. 2n+ 2 Stutzstellen x0, . . . , nn bzw.x0, . . . , xn, xn+1, . . . , x2n+1 ∈ [a, b]. Fur die zugeordneten Quadraturformeln giltdann
I(f)− I2n+1(f) = I(f)−2n+1∑
i=0
f [x0, . . . , xi]
∫ b
a
i−1∏
j=0
(x− xj) dx
= I(f)− In(f)−2n+1∑
i=n+1
f [x0, . . . , xi]
∫ b
a
i−1∏
j=0
(x− xj) dx.
Wir schreiben fur i = n + 1, . . . , 2n+ 1:∫ b
a
i−1∏
j=0
(x− xj) dx =
∫ b
a
n∏
j=0
(x− xj)i−1∏
j=n+1
(x− xj) dx
Die Polynome
1, (x− xn+1), (x− xn+1)(x− xn+2), . . . ,
2n∏
j=n+1
(x− xj)
66 KAPITEL 4. NUMERISCHE INTEGRATION
bilden eine Basis des Pn. Wahlen wir nun die ersten n+ 1 Stutzstellen x0, . . . , xn
aus [a, b] derart, dass
∫ b
a
q(x)n∏
j=0
(x− xj) dx = 0 ∀q ∈ Pn,
so folgtI(f)− In(f) = I(f)− I2n+1(f),
d.h. die interpolatorische Quadraturformel In ist exakt fur Polynome aus P2n+1.Damit ergibt sich die Frage nach der Existenz eines Polynoms n + 1-ten Gradesder Form
p(x) = xn+1 + r(x), r ∈ Pn,
das im L2(a, b) orthogonal zum Pn ist und reelle Nullstellen besitzt, die im Inter-vall [a, b] liegen.
Wir betrachten im folgenden die Aufgabe etwas allgemeiner und legen eingewichtetes Skalarprodukt der Form
(f, g)ω :=
∫ b
a
ω(x)f(x)g(x) dx
mit einer integrablen Gewichtsfunktion ω > 0, x ∈ (a, b) zugrunde. Wie ublichsei dann
‖f‖ω :=
(∫ b
a
ω(x)f 2(x) dx
)1/2
.
Mit Hilfe des Schmidtschen Orthogonalisierungsverfahrens gewinnen wir aus denPolynomen 1, x, x2, . . . Orthogonalpolynome pn, n = 0, 1, . . . :
p0(x) = 1 p0(x) = ‖p0‖−1ω p0(x),
pk(x) = xk −k−1∑
j=0
(xk, pj)ω pj(x), pk(x) = ‖pk‖−1ω pk(x), k = 1, . . . , n+ 1.
Dann ist p0, . . . , pn+1 ein Orthogonalsystem und p0, . . . , pn+1 ein Orthonor-malsystem in Pn+1.
Theorem 4.3.1 Die bezuglich des gewichteten Skalarproduktes (·, ·)ω orthogona-len Polynome pn, n ≥ 1, besitzen reelle, einfache Nullstellen, die alle im Innerndes Intervalls [a, b] liegen.
Beweis. Wir definieren die Menge
Nn := λ ∈ (a, b) : λ Nullstelle ungerader Vielfachheit von pn
4.3. GAUSSSCHE QUADRATURFORMELN 67
und setzen
q(x) := 1 fur Nn = ∅
q(x) :=
m∏
i=1
(x− λi) fur Nn = λ1, . . . , λm.
Dann ist das Polynom q · pn ∈ Pn+m reell und hat in (a, b) keinen Vorzeichen-wechsel. Somit gilt
(pn, q) =
∫ b
a
ω(x)pn(x)q(x) dx 6= 0.
Angenommen die Anzahl m der reellen Nullstellen in (a, b) ist kleiner als n. Dannist q ∈ Pn−1 und (pn, q) = 0 im Widerspruch zur obigen Beziehung.
Die orthogonalen Polynome pn bezuglich des Gewichtes ω(x) = 1 auf [a, b]heißen Legendre-Polynome; ublicherweise betrachtet man diese auf dem Referenz-intervall [−1,+1] und bezeichnet sie mit Ln. Wir konnen nun die Nullstellenλ0, . . . , λn des (n + 1)-ten Legendre-Polynoms Ln+1 als Stutzstellen einer inter-polierenden Quadraturformel auf dem Intervall [−1,+1] verwenden:
∫ +1
−1
f(x) dx ≈ In(f) :=n∑
i=0
αif(λi), αi =
∫ +1
−1
n∏
j=0j 6=i
x− λj
λi − λjdx
Theorem 4.3.2 Es gibt genau eine interpolatorische Quadraturformel zu n + 1paarweise verschiedenen Stutzstellen uber dem Intervall [−1, 1], die Polynomevom Grade kleiner oder gleich 2n + 1 exakt integriert. Ihre Stutzstellen sind dieNullstellen λ0, . . . , λn ∈ (−1, 1) des (n+ 1)-ten Legendre-Polynoms Ln+1 ∈ Pn+1,und ihre Gewichte genugen der Beziehung
αi =
∫ +1
−1
n∏
j=0j 6=i
(x− λj
λi − λj
)2
dx > 0, i = 0, . . . , n.
Fur f ∈ C2n+2[−1, 1] gilt die Restglieddarstellung
Rn(f) =f (2n+2)(ξ)
(2n + 2)!
∫ +1
−1
n∏
j=0
(x− λj)2 dx, ξ ∈ (−1, 1).
Beweis. Das Legendre-Polynom Ln+1 ist orthogonal zu Pn und hat mit seinen(reellen) Nullstellen λ0, . . . , λn die Darstellung
Ln+1(x) =
n∏
i=0
(x− λi).
68 KAPITEL 4. NUMERISCHE INTEGRATION
Nach den obigen Vorbetrachtungen integriert die interpolatorische Quadraturfor-mel dann Polynome vom Grade kleiner gleich 2n+ 1 exakt. Zur Bestimmung derGewichte αi betrachten wir die Polynome
li(x) =n∏
j=0j 6=i
(x− λj
λi − λj
), i = 0, . . . , n.
Da l2i ∈ P2n folgt
0 <
∫ 1
−1
l2i (x) dx =
n∑
j=0
αjli(λj)2 = αi.
Zur Eindeutigkeit der Gaußschen Quadraturformel sei angenommen, es gabe einezweite Formel
I∗n(f) =
n∑
i=0
α∗i f(λ∗i ),
die Polynome vom Grade kleiner oder gleich 2n+1 exakt integriert. Ersetzt manin der Definition von li die Nullstellen λi durch λ∗i erhalt man Polynome l∗i ∈ Pn,mit denen man wie oben α∗
i > 0 zeigen kann. Somit ware
0 =
∫ 1
−1
1
α∗i
l∗i (x)Ln+1(x) dx =n∑
j=0
α∗j
α∗i
l∗i (λj∗)Ln+1(λj∗) = Ln+1(λi∗).
Aus der eindeutigen Bestimmtheit der Nullstellen λi von Ln+1 folgt damit λi = λ∗iund αi = α∗
i .Wir mussen noch die Darstellung fur das Restglied zeigen. Nach den Aussagen
zur Hermite-Interpolation gibt es ein eindeutig bestimmtes Polynom H ∈ P2n+1,das die Hermitesche Interpolationsaufgabe
H(λi) = f(λi), H ′(λi) = f ′(λi), i = 0, . . . , n,
lost und fur f ∈ C2n+2[−1, 1] die Restglieddarstellung
f(x)−H(x) = f [λ0, λ0, . . . , λn, λn, x]
n∏
i=0
(x− λi)2
hat. Die Anwendung der Gaußschen Quadraturformel auf H ergibt dann
I(f)−In(f) = I(f −H)− In(f −H)
=
∫ 1
−1
f [λ0, λ0, . . . , λn, λn, x]n∏
i=0
(x− λi)2 dx−
n∑
i=0
αi (f(λi)− h(λi))
=f (2n+2)(ξ)
(2n+ 2)!
∫ 1
−1
n∏
i=0
(x− λi)2 dx.
4.3. GAUSSSCHE QUADRATURFORMELN 69
Im letzten Schritt haben wir die Intergaldarstellung dividierter Differenzen undden Mittelwertsatz der Integralrechnung benutzt.
Die Legendre-Polynome Ln ∈ Pn lassen sich auf [−1, 1] in der Form (Formelvon Rodriguez)
Ln(x) =1
2nn!
dn
dxn(x2 − 1)n, n = 0, 1, . . . ,
darstellen und genugen der rekursiven Beziehung
L0(x) = 1, L1(x) = x, Ln+1(x) =2n+ 1
n + 1xLn(x)− n
n + 1Ln−1(x), n ≥ 1.
Ihre Nullstellen werden analytisch beziehungsweise (fur n > 3) numerisch be-stimmt und konnen der Tabelle 4.3 entnommen werden. Die Skalierung ergibtsich aus Ln(1) = 1 fur alle n ∈ N. Wir haben insbesondere
L2(x) =3x2 − 1
2λ0 = − 1√
3, λ1 =
1√3,
L3(x) =5x3 − 3x
2λ0 = −
√3
5, λ1 = 0, λ2 =
√3
5.
Zur Bestimmung der Gewichte αi nutzen wir die Eigenschaft, dass
∫ 1
−1
p(x) dx =n∑
j=0
αj p(λj) ∀p ∈ P2n+1.
Setzen wir
p(x) =Ln+1(x)
x− λiLn(x) ∈ P2n,
so folgt wegen p(λj) = 0 fur j 6= i und p(λi) = L′n+1(λi)Ln(λi)
∫ 1
−1
Ln+1(x)
x− λi
Ln(x) dx = αiL′n+1(λi)Ln(λi).
Wir haben mit einem gewissen Polynom Qn−1(x) vom Grade kleiner gleich n− 1
Ln+1(x)
x− λi=
2n+ 1
n+ 1Ln(x) +Qn−1(x),
woraus wegen der Orthogonalitat Ln ⊥ Pn−1 die Beziehung
2n+ 1
n + 1(Ln, Ln) = αiL
′n+1(λi)Ln(λi)
70 KAPITEL 4. NUMERISCHE INTEGRATION
Tabelle 4.3: Stutzstellen und Gewichte der Gaußquadratur auf [−1, 1].
n λ0, . . . , λn−1 α0, . . . , αn−1
1 0 2
2 ±1/√
3 1
3 ±√
3/5 5/90 8/9
4 ±0.8611363116 0.3478548451±0.3399810436 0.6521451549
folgt. Aus der Formel von Rodriguez folgt durch mehrfache partielle Integration
(Ln, Ln) =(−1)n
22n(n!)2
∫ 1
−1
(x2 − 1)n d2n
dx2n(x2 − 1)n dx
=(−1)n(2n)!
22n(n!)2
∫ 1
−1
(x− 1)n(x+ 1)n dx
=(2n)!
22n(n!)2
n!
(n+ 1) · · · (2n)
∫ 1
−1
(x+ 1)2n dx
=2
2n+ 1.
Fur die Gewichte folgt damit
αi =1
L′n+1(λi)Ln(λi)
· 2
n+ 1,
wobei λj , j = 0, 1, . . . , n, die Nullstellen des Legendre-Polynoms Ln+1 sind. Furdas Restglied gilt
Rn(f) =(n+ 1)!4 22n+3
(2n+ 2)!3 (2n+ 3)f (2n+2)(ξ).
Die Gaußschen Quadraturformeln uber einem beliebigen Intervall [a, b] gewinntman durch Anwendung einer Koordinatentransformation, die das Referenzinter-vall [−1, 1] vermoge t 7→ x = (b− a)t/2 + (b+ a)/2 auf [a, b] abbildet.
Theorem 4.3.3 Sei In(f) die n + 1-punktige Gauß-Legendre-Formel zur Inter-gation von f auf [−1, 1]. Dann gilt fur jede Funktion f ∈ C[−1, 1]
∫ 1
−1
f(x) dx− In(f)→ 0 fur n→∞.
4.3. GAUSSSCHE QUADRATURFORMELN 71
Beweis. Betrachten wir die Gauß-Legendre-Formeln fur beliebiges n,
In(f) =
n∑
i=0
α(n)i f(λ
(n)i ), α
(n)i > 0,
n∑
i=0
α(n)i = 2.
Nach dem Weierstraßschen Approximationssatz, gibt es zu jeder vorgegebenenToleranzschranke ε > 0 und jeder auf [−1, 1] stetigen Funktion f ein Polynom pε
mit
‖f − pε‖ = maxx∈[−1,1]
|f(x)− pε(x)| <ε
4.
Wir splitten den Fehler in die drei Anteile
|I(f)− In(f)| ≤ |I(f − pε)|+ |I(pε)− In(pε)|+ |In(pε − f)|.
Der erste wird mittels Intergralabschatzung (Lange des Integrationsintervalls xBetragsmaximum des Integranden) abgeschatzt, also
|I(f − pε)| ≤ 2ε
4.
Der zweite verschwindet fur hinreichend grosse n, da dann das Polynom pε exaktintegriert wird. Fur den dritten gilt
|In(pε − f)| ≤n∑
i=0
α(n)i ‖pε − f‖ ≤
ε
4
n∑
i=0
α(n)i =
ε
2.
Fassen wir die drei Abschatzungen zusammen, erhalten wir die Konvergenz derFolge In(f) gegen I(f).
Die Methode der Konstruktion der Gauß-Legendre-Formeln zur optimalenBerechnung von I(f) kann auf den Fall von Integralen
Iω(f) =
∫ b
a
ω(x)f(x) dx
mit einer integrable Gewichtsfunktion ω(x) > 0 auf (a, b) ubertragen werden.Hierzu verwendet man als Stutzstellen gerade die Nullstellen der bezuglich desgewichteten Skalarproduktes
(f, g)ω =
∫ b
a
ω(x)f(x)g(x) dx
orthogonalen Polynome.
72 KAPITEL 4. NUMERISCHE INTEGRATION
Beispiel: Wir betrachten [a, b] = [−1, 1] und ω(x) = (1− x2)−1/2. Die ortho-gonalen Polynome sind in diesem Fall die Tschebyscheff-Polynome Tn ∈ Pn, diedurch die rekursive Beziehung
T0(x) = 1, T1(x) = x, Tn+1(x) = 2xTn(x)− Tn−1(x), n ≥ 1,
bestimmt sind. Die Stutzstellen und Gewichte der zugehorigen Quadraturformelnsind
λi = cos
(2i+ 1
n+ 1· π
2
), αi =
π
n+ 1, i = 0, . . . , n.
Die Restglieder haben die Form
Rn(f) =2π
(2n+ 2)!
(1
2
)2n+2
f (2n+2)(ξ), ξ ∈ (−1, 1).
Der Fall n = 2 ergibt
∫ 1
−1
ω(x)f(x) dx =π
3
f
(−√
3
2
)+ f(0) + f
(√3
2
)+
π
23040f (4)(ξ).
Kapitel 5
Approximation
Bei der Interpolationsaufgabe wurde eine stetige Funktion durch ein Interpola-tionspolynom angenahert. Die Forderung war dabei, dass Funktion und Interpola-tionspolynom in ausgewahlten Stutzstellen – den Interpolationsknoten – uberein-stimmten. Wir betrachten in diesem Abschnitt die Aufgabe, dass eine stetigeFunktion in gewissem Sinne bestmoglich approximiert wird. Dazu sei im folgen-den die Menge der auf [a, b] stetigen reell- bzw. komplexwertigen Funktionen alsVektorraum uber den Zahlkorper K = R bzw. K = C aufgefasst. Gegeben seieine Funktion f ∈ C[a, b] sowie eine Teilmenge S ⊂ C[a, b], deren Elemente zurApproximation von f dienen sollen.
Die Aufgabe der besten Approximation einer Funktion f ∈ C[a, b] durchElemente aus S lautet:
Finde p ∈ S mit ‖f − p‖ = infq∈S‖f − q‖.
Die Gute der Approximation wird dabei in der Norm ‖ · ‖ gemessen.
5.1 Gauß-Approximation
Bei der Gauß-Approximation bzw. der Approximation im quadratischen Mittelverwendet man die durch das Skalarprodukt
(f, g) =
∫ b
a
f(x)g(x) dx
erzeugte Norm
‖f‖ =
(∫ b
a
|f(x)|2 dx)1/2
.
Versehen mit dem Skalarprodukt (·, ·) wird C[a, b]) zu einem unitaren Raum.
73
74 KAPITEL 5. APPROXIMATION
Theorem 5.1.1 Seien H ein unitarer Raum und S ⊂ H ein endlich dimensio-naler Teilraum. Dann existiert zu jedem f ∈ H eine eindeutig bestimmte besteApproximation p ∈ S.
Beweis. Sei ψ1, . . . , ψn eine Basis von S, also n = dimS. Mit dem Schmidt-schen Orthogonalisierungsverfahren
ϕ1 = ψ1, ϕ1 = ‖ϕ1‖−1ϕ1,
ϕi = ψi −i−1∑
j=1
(ψi, ϕj)ϕj, ϕi = ‖ϕi‖−1ϕi,
erzeugt man ein Orthonormalsystem ϕ1, . . . , ϕn in S, d.h. wir haben
(ϕi, ϕj) = δij , i, j = 1, . . . , n.
In der Tat ist ϕ1 wegen ψ1 6= 0 wohldefiniert. Seien nun ϕ1, . . . , ϕm−1 wohlde-finiert und orthonormal. Ware nun ϕm = 0, so wurde
ψm =
m−1∑
j=1
(ψm, ϕj)ϕj
gelten, d.h. ψm ware linear von ψ1, . . . , ψm−1 linear abhangig im Widerspruch zurAnnahme, dass ψ1, . . . , ψn eine Basis von S sei. Damit ist auch ϕm = ‖ϕm‖−1ϕm
wohldefiniert und fur k = 1, . . . , i− 1 gilt
(ϕm, ϕk) = (ψi, ϕk)−i−1∑
j=1
(ψi, ϕj) (ϕj, ϕk) = (ψi, ϕk)− (ψi, ϕk) = 0.
Als Orthonormalsystem ist ϕ1, . . . , ϕn automatisch eine Basis von S, denn
n∑
j=1
αjϕj = 0 ⇒ αk =n∑
j=1
αj(ϕj , ϕk) = 0.
Jedes p ∈ S kann damit eindeutig in der Form
p =n∑
j=1
αjϕj
dargestellt werden, die Koeffizienten αj sind bestimmt durch
(p, ϕk) =
n∑
j=1
αj(ϕj, ϕk) = αk.
5.1. GAUSS-APPROXIMATION 75
Fur beliebiges
q =
n∑
j=1
βjϕj ∈ S
gilt nun
‖f − q‖2 = (f − q, f − q) =
(f −
n∑
j=1
βjϕj , f −n∑
j=1
βjϕj
)
= (f, f)−n∑
j=1
βj(ϕj, f)−n∑
j=1
βj(f, ϕj) +n∑
i,j=1
βiβj(ϕi, ϕj)
= ‖f‖2 − 2
n∑
j=1
Re[βj(ϕj , f)] +
n∑
j=1
|βj|2.
Wegen
|βj − (f, ϕj)|2 = [βj − (f, ϕj)] [βj − (f, ϕj)] = |βj|2 − 2Re[βj(ϕj , f)] + |(f, ϕj)|2
folgt
‖f − q‖2 = ‖f‖2 −n∑
j=1
|(f, ϕj)|2 +n∑
j=1
|βj − (f, ϕj)|2,
somit ist q ∈ S genau dann beste Approximation von f , wenn fur die Koeffizien-ten βj gilt βj = (f, ϕj), j = 1, . . . , n.
Die beste Approximation p ∈ S von f ∈ C[a, b] im quadratischen Mittel kannalso in der Form
p =n∑
j=1
(f, ϕj)ϕj
dargestellt werden, wobei ϕ1, . . . , ϕn ein Orthonormalsystem von S ist. DerFehler der besten Approximation berechnet sich demnach als
‖f − p‖ =
(‖f‖2 −
n∑
j=1
|(f, ϕj)|2)1/2
.
Theorem 5.1.2 Die beste Approximation p ∈ S von f kann aquivalent durchdie Orthogonalitatseigenschaft f − p ⊥ S, i.e.
(f − p, ϕ) = 0 ∀ϕ ∈ S
charakterisiert werden.
76 KAPITEL 5. APPROXIMATION
Beweis. Sei ϕ1, . . . , ϕn ein Orthonormalsystem in S. Dann ist q genau dannbeste Approximation von f , wenn
(f − q, ϕi) =
(f −
n∑
j=1
βjϕj, ϕi
)= (f, ϕi)− βi = 0, i = 1, . . . , n.
Andererseits ist
(f − q, ϕ) = 0, ∀ϕ ∈ S ⇔ (f − q, ϕi) = 0, i = 1, . . . , n.
Die Orthogonalitatseigenschaft ermoglicht die Bestimmung der besten Appro-ximation in Bezug auf eine beliebige Basis ψ1, . . . , ψn von S. Sei
p =n∑
j=1
γjψj
die eindeutige Darstellung der besten Approximation von f ∈ C[a, b] bezuglichdieser Basis. Dann sind die Koeffizienten γj wegen der Orthogonalitatseigenschaftder besten Approximation Losung des linearen Gleichungssystems
0 =
(f −
n∑
j=1
γjψj , ψi
)= (f, ψi)−
n∑
j=1
γj(ψj , ψi).
Die Koeffizientenmatrix A = (ψj , ψi) ist die Gramsche Matrix der Basis ψ1, . . . , ψn;sie ist hermitisch und positiv definit,
(Aγ, γ) =
n∑
i,j=1
aijγjγi =
n∑
i,j=1
(ψj , ψi)γjγi =
(n∑
j=1
ψj ,
n∑
i=1
ψi
)= (p, p) = ‖p‖2 > 0
fur jedes p 6= 0.
Betrachten wir als Beispiel fur die Menge S die Polynome vom Grade kleinergleich n auf dem Intervall [−1,+1]. Verwenden wir die Standardbasis des Pn, also1, x, x2, . . . , xn, so sind die Eintrage der Gramschen Matrix
aij =
∫ 1
−1
xj xi dx =
0 falls i+ j gerade2
i+ j + 1falls i+ j ungerade.
Die zugehorige Matrix ist zwar regular, aber extrem schlecht konditioniert, furn = 10 berechnet MATLAB cond(A) = 2.2 ·1033. Wesentlich besser geeignet sind
5.1. GAUSS-APPROXIMATION 77
die orthogonalen Legendre-Polynome Ln(x), die wir im Zusammenhang mit derGaußquadratur kennengelernt haben. Wegen
(Ln, Ln) =2
2n + 1, n = 0, 1, . . .
bilden die Polynome
ϕk(x) =
√2k − 1
2Lk−1(x), k = 1, 2, . . .
ein orthognormiertes System. Die Berechnung der besten Approximation p vonf gestaltet sich nun wesentlich einfacher
p(x) =n∑
j=1
∫ 1
−1
f(t)ϕj(t) dtϕj(x).
Die Maximalabweichung
‖f − p‖∞ = max−1≤x≤1
|f(x)− p(x)|
kann jedoch (vor allem zu den Intervallenden hin) sehr groß werden. Um diesenEffekt zu unterdrucken, verwendet man das gewichtete Skalarprodukt
(f, g)ω =
∫ 1
−1
ω(x)f(x)g(x) dx, ω(x) =1√
1− x2,
wodurch eine starkere Bindung an den Intervallenden erreicht wird. Durch Ortho-normierung von 1, x, . . . , xn−1 bezuglich des mit 1/
√1− x2 gewichteten Ska-
larprodukts erhalt man die Polynome
ϕ1 =
√1
π, ϕk =
√2
πTk−1(x), k = 2, . . . , n,
mit den Tschebyscheff-Polynomen Tk.
Theorem 5.1.3 Die Tschebyscheff-Polynome haben fur x ∈ [−1, 1] die Gestalt
Tk(x) = cos (k arccos(x)) , k = 0, 1, 2, . . . ,
und genugen den Beziehungen
T0(x) = 1, T1(x) = x, Tk+1(x) = 2xTk(x)− Tk−1(x)
∫ 1
−1
1√1− x2
Ti(x)Tj(x) dx =
π fur i = j = 0π/2 fur i = j 6= 00 fur i 6= j.
78 KAPITEL 5. APPROXIMATION
Beweis. Wir zeigen zunachst die Rekursionsbeziehung fur gk(x) := cos (k arccos(x)).Wir haben unmittelbar g0(x) = 1 und g1(x) = x. Aus den Additionstheoremenleitet man
cos(α + β) + cos(α− β) = 2 cosα cos β
oder mit x = α + β und y = α− β
cosx+ cos y = 2 cosx+ y
2cos
x− y2
ab. Somit ist
gk+1(x) + gk−1(x) = cos((k + 1) arccos(x)) + cos((k − 1) arccos(x))
= 2 cos(k arccos(x)) cos(arccos(x)) = 2xgk(x).
Damit ist gk genau ein Polynom k-ten Grades und genugt der obigen Rekursi-onsbeziehung. Weiter gilt mit der Substitution x = cos t
∫ 1
−1
1√1− x2
gj(x)gi(x) dx =
∫ π
0
cos jt cos it dt
=1
2
∫ π
0
cos(i+ j)t+ cos(i− j)t dt
woraus obige Orthogonalitatseigenschaft folgt.
5.2 Tschebyscheff-Approximation
Im folgenden betrachten wir nur reell-wertige Funktionen. Im Unterschied zurGauß-Approximation, die die beste Approximation im quadratischen Mittel bzw.im gewichteten quadratischen Mittel sucht, verwendet die Tschebyscheff-Appro-ximation direkt die Maximumnorm
‖f‖∞ = maxx∈[a,b]
|f(x)|.
Diese Norm wird nicht durch ein Skalarprodukt erzeugt, die Existenz einer bestenApproximation kann daher nicht aus Theorem 5.1.1 gefolgert werden. Tatsachlichist im allgemeinen die beste Approximation nicht einmal eindeutig bestimmt, wiedas folgende Beispiel zeigt:
Beispiel: Seien [a, b] = [0, 1] und f(x) = 1 auf [0, 1]. Wir wollen die Funktion f inder Maximumnorm bestmoglich durch ein Element der eindimensionalen MengeS = pα(x) = αx : α ∈ R approximieren. Wir haben ‖f − pα‖∞ ≥ 1 fur allep ∈ S und ‖f − pα‖∞ = 1 fur 0 ≤ α ≤ 2.
5.2. TSCHEBYSCHEFF-APPROXIMATION 79
Theorem 5.2.1 Sei E ein normierter Vektorraum und S ⊂ E ein endlich di-mensionaler Teilraum. Dann gibt es zu jedem f ∈ E eine beste Approximationp ∈ S:
‖f − p‖ = minq∈S‖f − q‖.
Beweis. Ein Element q0 ∈ S mit ‖q0‖ > 2‖f‖ kann keine beste Approximationsein, denn
‖f − q0‖ ≥ ‖q0‖ − ‖f‖ > ‖f‖ = ‖f − 0‖ ≥ infq∈S‖f − q‖.
Somit haben wirinfq∈S‖f − q‖ = inf
q∈S,‖q‖≤2‖f‖‖f − q‖.
Nun ist die Abbildung F : E → R mit F (q) = ‖f − q‖ auf der kompakten MengeA = q ∈ S : ‖q‖ ≤ 2‖f‖ stetig, denn
|F (q1)− F (q2)| =∣∣∣‖f − q1‖ − ‖f − q2‖
∣∣∣ ≤ ‖q1 − q2‖.
Nach dem Satz von Weiserstraß nimmt F auf A das Infimum als Funktionswertan, es gibt also eine beste Approximation p von f .
Im Hinblick auf die Eindeutigkeit fuhren wir folgende Definition ein.
Definition: Wir nennen einen linearen normierten Raum E streng normiert,wenn in der Dreiecksungleichung
‖q1 + q2‖ ≤ ‖q1‖+ ‖q2‖
das Gleichheitszeichen fur q1 6= 0, q2 6= 0, nur im Fall q2 = αq1 mit positivem αgilt. Beachte, dass der R3 versehen mit der euklidischen Norm streng normiertist. Es ist sogar jeder unitare Raum streng normiert. Zum Beweis nehmen wiran, dass
‖q1 + q2‖ = ‖q1‖+ ‖q2‖, q1 6= 0, q2 6= 0
gelte. Durch Quadrieren erhalten wir
‖q1‖2+2(q1, q2)+‖q2‖2 = (q1+q2, q1+q2) = ‖q1+q2‖2 = ‖q1‖2+2‖q1‖ ‖q2‖+‖q2‖2.
Also gibt es zwei Elemente yi = qi/‖qi‖, i = 1, 2, der Lange 1, deren Skalarprodukt(y1, y2) = 1 ergibt. Wegen
‖y1 − y2‖2 = (y1 − y2, y1 − y2) = ‖y1‖2 − 2(y1, y2) + ‖y2‖2 = 1− 2 + 1 = 0
muss aber y1 = y2 bzw. q2 = (‖q2‖/‖q1‖) q1 sein.
80 KAPITEL 5. APPROXIMATION
Theorem 5.2.2 Ist E streng normiert, gibt es hochstens ein Element p ∈ S derbesten Approximation von f ∈ E.
Beweis. Nehmen wir an, dass zwei verschiedene Elemente q1, q2 der besten Ap-proximation von f in S existieren. Dann ist
‖f − q1‖ = ‖f − q2‖ = m > 0,
da andernfalls q1 = q2 = f ware. Nun ist
m ≤ ‖f − q1 + q22‖ = ‖1
2(f − q1) +
1
2(f − q2)‖
≤ 1
2‖(f − q1)‖+
1
2‖(f − q2)‖ = m.
Nun war E streng normiert, also muss f−q1 = α(f−q2) bzw. (1−α)f = q1−αq2gelten. Ware α 6= 1, so ist f eine Linearkombination von Elementen aus S undm = 0 im Widerspruch zur Annahme. Somit ist α = 1 und folglich q1 = q2.
Leider ist der lineare normierte Raum C[a, b] versehen mit der Supremums-norm nicht streng normiert. Hierzu betrachten wir [a, b] = [0, 1] mit q1 = 1 undq2 = x. Es gilt namlich
2 = ‖q1 + q2‖∞ = ‖q1‖∞ + ‖q2‖∞ = 1 + 1.
Die Eindeutigkeit der Tschebyscheff-Approximation wird durch die HaarscheBedingung (H) an den Ansatzraum S ⊂ C[a, b] mit dimS = n garantiert. Wirbetrachten zunachst die Losbarkeit der allgemeinen Interpolationsaufgabe
Finde p ∈ S = span (ϕ1, . . . , ϕn) mit dim S = n, so dass zu paarweiseverschiedenen Knoten xj , j = 1, . . . , n, und Werten yj, j = 1, . . . , n,die Interpolationsbedingung p(xj) = yj, j = 1, . . . , n, gilt.
Im Unterschied zur im Abschnitt 3.1 betrachteten Polynominterpolation hat dieallgemeine Interpolationsaufgabe nicht notwendig eine Losung. So verschwin-det beispielsweise jede Linearkombination der Funktionen ϕ1(x) = sin x undϕ2(x) = sin 2x an der Stelle x = x1 = 0. Somit ist die Interpolationsaufga-be fur alle y1 6= 0 unlosbar. Im allgemeinen Fall ist die Interpolationsaufgabegenau dann fur beliebige Werte yj, j = 1, . . . , n losbar, wenn die Haar-MatrixH = (ϕi(xj)) nicht singular ist. Dies ist, wie wir oben gesehen haben, nichtnotwendig fur jeden Satz paarweise verschiedener Knoten x1, . . . , xn der Fall. Dielineare Unabhangigkeit der Funktionen ϕ1, . . . , ϕn sichert aber, dass es zumindesteinen Satz paarweise verschiedener Knoten x1, . . . , xn gibt, fur den det H 6= 0 gilt.Im oben angegebenen Beispiel S = span (sin x, sin 2x) konnte man beispielsweisex1 = π/4 und x2 = π/2 nehmen.
5.2. TSCHEBYSCHEFF-APPROXIMATION 81
Theorem 5.2.3 (Haar) Zu jeder gegebenen Funktion f ∈ C[a, b] existiert genaudann ein eindeutig bestimmtes p ∈ S,
(H) wenn jedes Element q ∈ S\0 des Unterraumes S ⊂ C[a, b] mit dimS = nmaximal n−1 Nullstellen in [a, b] hat. Das S erzeugende System ϕ1, . . . , ϕn
heißt dann Tschebycheff-System.
Beweis. Siehe z.B. I.S. Beresin, N.P. Shidkow: Numerische Methoden I, Deut-scher Verlag der Wissenschaften, Berlin 1970.
So ist das System
S = span(sin x, sin 2x) ⊂ C[0,π
3
], dim S = 2,
kein Tschebycheff-System, denn q(x) = sin x− sin 2x verschwindet beispielsweisebei x = 0 und x = π/3.
Aquivalent zur Haarschen Bedingung ist die Forderung
(H’) Fur jeden Satz paarweise verschiedener Knotenwerte a ≤ x1 < x2 < · · · <xn ≤ b ist die Interpolationsaufgabe p(xi) = yi, i = 1, . . . , n mit beliebigenWerten y1, . . . , yn ∈ R eindeutig durch ein p ∈ S losbar.
Beachte, dass (H’) gleichbedeutend mit det H 6= 0 fur jeden Satz paarweise ver-schiedener Knotenwerte x1, . . . , xn ist.
Beweis der Aquivalenz. Wir zeigen zunachst (H ′) ⇒ (H): Sei q ∈ S\0und nehmen wir indirekt an, q habe n Nullstellen xi, i = 1, . . . , n. Die eindeutigeLosbarkeit der Interpolationsaufgabe q(xi) = 0 liefert q(x) = 0 fur alle x ∈ [a, b]im Widerspruch zu q ∈ S\0.
Umgekehrt gilt ¬(H ′)⇒ ¬(H): Wir haben also einen Satz paarweise verschie-dener Knotenwerte x1, . . . , xn und Werte y1, . . . , yn, fur den die Interpolations-aufgabe q(xi) = yi, i = 1, . . . , n nicht eindeutig in S losbar ist. Dann muss dieDeterminante der Haar-Matrix verschwinden und die Interpolationsaufgabe mithomogene Werten yi = 0, i = 1, . . . , n ist nicht eindeutig losbar. Neben q(x) = 0gibt es also ein q ∈ S\0 mit q(xi) = 0, also gibt es ein Element q ∈ S\0, dasmindestens n Nullstellen besitzt.
Beispiele: Die Polynomraume S = Pn−1 erfullen die Haarsche Bedingung aufjedem Intervall. Im Fall 0 ∈ [a, b] und S = span (x, x2, . . . , xn) ist die Bedingungjedoch nicht erfullt.
Theorem 5.2.4 (Alternantensatz von Tschebyscheff) Fur den TeilraumS ⊂ C[a, b] mit dimS = n sei die Haarsche Bedingung erfullt. Dann ist die
82 KAPITEL 5. APPROXIMATION
Tschebyscheff–Approximation p einer Funktion f ∈ C[a, b] durch folgende Eigen-schaft charakterisiert. Es gibt mindestens n + 1 Stellen a ≤ x1 < x2 < · · · <xn+1 ≤ b (Alternante genannt), so dass fur die Fehlerfunktion e(x) = f(x)−p(x)gilt:
|e(xi)| = ‖e‖∞, e(xi) = −e(xi+1); i = 1, . . . , n
Insbesondere ist die beste Approximation eindeutig bestimmt.
Beweis. Siehe etwa G. Maeß Vorlesungen uber Numerische Mathematik II,Akademie-Verlag Berlin 1988. Wir wollen hier nur die Eindeutigkeit der bestenApproximation fur den Spezialfall S = Pn−1 zeigen. Seien qi, i = 1, 2, zwei besteApproximationen mit den Fehlerfunktionen ei = f − qi. Fur λ ∈ (0, 1) ist dann‖λe2‖∞ < ‖e1‖∞. Aufgrund der Existenz einer Alternante schneidet der Graphvon e1 den von λe2 mindestens n-mal. Jede der Funktionen ϕλ(x) = e1(x)−λe2(x)hat somit mindestens n (paarweise verschiedene) Nullstellen. Gehen wir zur Gren-ze λ→ 1 uber, so muss ϕ1(x) = e1(x)− e2(x) = q2(x)− q1(x) ∈ Pn−1 mindestensn (ihrer Vielfachheit entsprechend oft gezahlte) Nullstellen haben, dies bedeutetzwangslaufig q2 = q1.
Zur Anwendung der Tschebyscheff–Approximation betrachten wir nun dasProblem der “optimalen” Wahl der Stutzstellen bei der Lagrange-Interpolation.Fur das Lagrange–Interpolationspolynom pn ∈ Pn einer Funktion f ∈ Cn+1[a, b]zu den Stutzstellen a ≤ x0 < x1 < · · · < xn ≤ b gilt die Fehlerdarstellung
f(x)− pn(x) =f (n+1)(ξx)
(n+ 1)!
n∏
j=0
(x− xj).
Die Stutzstellen sollen nun so gewahlt werden, dass das Maximum
maxx∈[a,b]
∣∣∣∣∣
n∏
j=0
(x− xj)
∣∣∣∣∣
minimal wird. Aus der Darstellung
n∏
j=0
(x− xj) = xn+1 − p, p ∈ Pn
sieht man, dass diese Aufgabe aquivalent dazu ist, die Tschebyscheff–Approxi-mation von f(x) = xn+1 bezuglich des Raumes S = Pn zu bestimmen. Nach demAlternantensatz ist der Betrag der Fehlerfunktion e = xn+1 − p an mindestensn+ 2 Stellen im Intervall [a, b] gleich dem Minimum.
Theorem 5.2.5 Auf dem Intervall [a, b] = [−1, 1] ist die Tschebyscheff–Approxi-mation p ∈ Pn zu f(x) = xn+1 gegeben durch
p(x) = xn+1 − 2−nTn+1(x)
5.2. TSCHEBYSCHEFF-APPROXIMATION 83
mit dem (n+ 1)-ten Tschebyscheff–Polynom
Tn+1(x) = cos[(n + 1) arccos(x)].
Die Nullstellen
xk = cos
(π
2· 2k + 1
n+ 1
), k = 0, . . . , n
von Tn+1 sind die “optimalen” Stutzstellen der Lagrange–Interpolation auf [−1, 1].
Bemerkung. Die “optimale” Wahl der Stutzstellen als Nullstellen des Tscheby-scheff–Polynoms Tn+1 bedeutet nicht, dass das zugeordnete Lagrangesche Inter-polationspolynom, ein Polynom bester Approximation in der Maximumnorm ist,denn in der Fehlerdarstellung
f(x)− pn(x) =f (n+1)(ξx)
2n(n+ 1)!Tn+1(x), x ∈ [−1, 1]
ist ξx noch von x abhangig.
Beweis von Theorem 5.2.5. Das Polynom Tn+1 ∈ Pn+1 hat n + 1 Nullstellen
xk = cos
(π
2· 2k + 1
n + 1
), k = 0, . . . , n.
Aus der rekursiven Beziehung
T0(x) = 1, T1(x) = x, Tn+1(x) = 2xTn(x)− Tn−1(x)
folgt, dass Tn+1(x) = 2nxn+1+q(x) mit einem geeigneten q ∈ Pn ist. Damit habenwir die Darstellung
2−nTn+1(x) =n∏
k=0
(x− xk),
aus der insbesondere
max−1≤x≤1
n∏
k=0
|x− xk| = max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− xk)
∣∣∣∣∣ = 2−n max−1≤x≤1
|Tn+1(x)| = 2−n
folgt. Nun nimmt Tn+1(x) = cos[(n + 1) arccos(x)] im Intervall [−1,+1] genaun+ 2 mal einen Extremwert an, abwechselnd ±1. Diese n+ 2 Extremalstellen
x∗k = coskπ
n+ 1, k = 0, . . . , n+ 1,
bilden eine Alternante fur die Approximation p(x) = xn+1 − 2−nTn+1(x) ∈ Pn
zu xn+1. Nach dem Alternantensatz ist somit p die eindeutig bestimmte beste
84 KAPITEL 5. APPROXIMATION
Approximation zu xn+1.
Es bleibt noch die Optimalitatseigenschaft
max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− xk)
∣∣∣∣∣ = min−1≤ζ0<ζ1<···<ζn≤1
max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− ζk)∣∣∣∣∣
zu zeigen. Mit den Nullstellen xk von Tn+1(x) gilt
max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− xk)
∣∣∣∣∣ = 2−n max−1≤x≤1
|Tn+1(x)|
= max−1≤x≤1
∣∣∣xn+1 − [xn+1 − 2−nTn+1(x)]∣∣∣ = min
q∈Pn
max−1≤x≤1
|xn+1 − q(x)|,
also
max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− xk)
∣∣∣∣∣ = minq∈Pn
max−1≤x≤1
|xn+1 − q(x)|
= min−1≤ζ0<ζ1<···<ζn≤1
max−1≤x≤1
∣∣∣∣∣
n∏
k=0
(x− ζk)∣∣∣∣∣ .
und damit die Optimalitatseigenschaft.
Kapitel 6
Nichtlineare Gleichungen
6.1 Nullstellen reeller Funktionen
Sei f : [a, b]→ R eine auf dem Intervall [a, b] stetige Funktion. Das einfachste Ver-fahren zur Bestimmung von Nullstellen von f beruht auf dem Zwischenwertsatzfur stetige Funktionen: Gibt es ein Teilintervall [a0, b0] ⊂ [a, b] mit f(a0)f(b0) < 0,so hat f in (a0, b0) mindestens eine Nullstelle.
Das Intervallschachtelungs-Verfahren erzeugt ausgehend von einen In-tervall [a0, b0] mit f(a0)f(b0) < 0 eine Folge von Intervallen [ai, bi], i = 0, 1, . . .die mindestens eine Nullstelle von f besitzen durch die Vorschrift
xi =1
2(ai + bi), f(xi) = 0 ⇒ STOP
f(ai)f(xi) < 0 ⇒ ai+1 = ai, bi+1 = xi
f(ai)f(xi) > 0 ⇒ ai+1 = xi, bi+1 = bi.
Im Fall, dass der Algorithmus nicht nach endlich vielen Schritten mit einer Null-stelle abbricht, gilt ai ≤ ai+1 ≤ bi+1 ≤ bi und
|bi+1 − ai+1| =1
2|bi − ai| = 2−(i+1)|b0 − a0|.
Die monotonen nach oben bzw. nach unten beschrankten Zahlenfolgen (an)n∈N
und (bn)n∈N konvergieren gegen ein x0 ∈ R, das wegen
f(x0)2 = lim
n→∞f(an)f(bn) ≤ 0
eine Nullstelle von f ist. Hat man also ein Ausgangsintervall [a0, b0] ⊂ [a, b] mitf(a0)f(b0) < 0 gefunden, so liefert das Verfahren fur stetige Funktionen immereine Nullstelle. Nachteil des Verfahrens ist die sehr langsame Konvergenz.
85
86 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
Ist die gegebene Funktion auf [a, b] stetig differenzierbar, so kann diese Zu-satzinformation zur effizienteren Berechnung einer Nullstelle verwendet werden.Das klassische Newton-Verfahren basiert auf Approximation der Funktion inUmgebung einer gesuchten Nullstelle durch das Taylorpolynom ersten Grades(graphisch durch die Tangente). Sei also xn eine Naherung fur die gesuchte Null-stelle x0. Die Gleichung f(x) = 0 wird nun ersetzt durch
f(x) ≈ f(xn) + f ′(xn)(x− xn) = 0
mit der Losung
xn+1 = xn −f(xn)
f ′(xn), n = 1, 2, . . .
Die Iteration ist durchfuhrbar, wenn die Ableitungswerte f ′(xn) nicht zu kleinwerden.
Theorem 6.1.1 Die Funktion f ∈ C2[a, b] habe im Innern von [a, b] eine Null-stelle x0 und es gelte
m := mina≤x≤b
|f ′(x)| > 0, M := maxa≤x≤b
|f ′′(x)| <∞.
Sei der Radius ρ der abgeschlossenen Kugel B(x0, ρ) klein genug gewahlt, so dass
B(x0, ρ) ⊂ [a, b] und q :=M
2mρ < 1.
Dann ist fur jeden Startwert x1 ∈ B(x0, ρ) die Newton-Iteration durchfuhrbarund die Folge (xn)n∈N der Newton-Iterierten konvergiert gegen die Nullstelle x0.Daruber hinaus gelten die a-priori Fehlerabschatzung
|xn+1 − x0| ≤2m
Mq2n
n ∈ N
sowie die a-posteriori Fehlerabschatzung
|xn+1 − x0| ≤1
m|f(xn+1)| ≤
M
m|xn+1 − xn|2.
Beweis. Wir beginnen mit einigen Vorbereitungen. Fur x, y ∈ [a, b] mit x 6= yfolgt aus dem Mittelwertsatz der Differentialrechnung mit einer Zwischenstelleξ ∈ (a, b) ∣∣∣∣
f(x)− f(y)
x− y
∣∣∣∣ = |f ′(ξ)| ≥ m > 0
und somit
|x− y| ≤ 1
m|f(x)− f(y)|. (6.1.1)
6.1. NULLSTELLEN REELLER FUNKTIONEN 87
Hieraus folgt insbesondere, dass die Nullstelle x0 ∈ [a, b] eindeutig bestimmt ist.Wir betrachten nun die Abbildung P : B(x0, ρ)→ R, gegeben durch
P (x) = x− f(x)
f ′(x),
mit deren Hilfe das Newton-Verfahren in der Form
xn+1 = P (xn) n = 1, 2, . . .
beschrieben werden kann. Mit Hilfe der Taylorschen Formel zur Entwicklung vonf(x0) an der Entwicklungsstelle x folgt
0 = f(x0) = f(x) + f ′(x)(x0 − x) +f ′′(ξ)
2(x0 − x)2,
wobei ξ ∈ (a, b) wieder eine geeignete Zwischenstelle bezeichnet. Fur x ∈ B(x0, ρ)folgt
P (x)− x0 = x− f(x)
f ′(x)− x0 = − 1
f ′(x)f(x) + f ′(x)(x0 − x)
=1
f ′(x)
f ′′(ξ)
2(x0 − x)2
|P (x)− x0| ≤M
2m|x− x0|2 ≤
(M
2mρ
)ρ < ρ, (6.1.2)
also liegt P (x) in der Kugel B(x0, ρ), vorausgesetzt dass bereits x in dieser Kugellag, d.h. xn ∈ B(x0, ρ) fur alle n ∈ N. Setzt man ρn = |xn−x0|M/(2m), so erhaltman aus (6.1.2) fur alle n ∈ N
ρn+1 =M
2m|xn+1 − x0| =
M
2m|P (xn+1)− x0| ≤
(M
2m|xn − x0|
)2
= ρ2n.
Somit haben wir
ρn+1 ≤ (ρ1)2n ⇒ |xn+1 − x0| ≤
2m
M(ρ1)
2n
,
woraus wegen ρ1 = |x1 − x0|M/(2m) ≤ ρM/(2m) = q < 1 die Konvergenzvon xn gegen die Nullstelle x0 folgt. Die letzte Abschatzung beinhaltet zugleichdie im Satz behauptete a-priori Fehlerabschatzung. Zum Beweis der a-posterioriFehlerabschatzung nutzt man die Taylorentwicklung von f(xn+1) an der Entwick-lungsstelle xn:
f(xn+1) = f(xn) + f ′(xn)(xn+1 − xn) +f ′′(ξ)
2(xn+1 − xn)2,
88 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
woraus unter Beachtung von f(xn) + f ′(xn)(xn+1 − xn) = 0 und (6.1.1)
|xn+1 − x0| ≤1
m|f(xn+1)− f(x0)| ≤
M
2m|xn+1 − xn|2
folgt.
Fur zweimal stetig differenzierbare Funktionen f gibt es zu jeder einfachenNullstelle x0, d.h. f(x0) = 0 und f ′(x0) 6= 0 eine abgeschlossenen Kugel B(x0, ρ),so dass die Voraussetzungen des Satzes erfullt sind. Das Problem beim Newton-Verfahren ist die Bestimmung eines im Einzugsbereich der Nullstelle liegendenStartwertes x1.
Beispiel: Newton-Verfahren zur WurzelberechnungDie Quadratwurzel einer Zahl a > 0 ist die Nullstelle der Funktion f(x) = x2−a.Die Newton-Iteration ist gegeben durch
xn+1 = xn −x2
n − a2xn
=1
2
(xn +
a
xn
).
Wir untersuchen die Konvergenz der Folge und stellen zunachst fest, dass fureinen beliebigen positiven Startwert x1 > 0 die Positivitat aller Folgengliederxn+1 folgt. Die Ungleichung zum arithmetischen und geometrischen Mittel ergibt
xn+1 =1
2
(xn +
a
xn
)≥ √a
wobei das Gleichheitszeichen nur fur xn =√a gilt. Startet man also mit einem
positiven x1 6=√a, so ist xn >
√a fur alle n ≥ 2. Eine derartige Folge (xn)n≥2
fallt streng monoton, denn
xn+1 − xn =1
2xn(a− x2
n) < 0.
Als streng fallende Folge, die nach unten beschrankt ist, konvergiert (xn)n≥2 gegeng, wobei g der Gleichung
g =1
2
(g +
a
g
)⇒ g =
√a
genugt. Fur hinreichend große n ist xn im Einzugsbereichs der Nullstelle x0 =√a.
Im konkreten Fall haben wir
xn+1 −√a =
1
2
(xn +
a
xn
)−√a =
1
2xn
(x2
n + a− 2√axn
)=
1
2xn
(xn −
√a)2,
|xn+1 −√a| ≤ 1
2√a|xn −
√a|2 n ≥ 2.
6.1. NULLSTELLEN REELLER FUNKTIONEN 89
Quadratische Konvergenz liegt bereits ab dem zweiten Folgeglied vor.
Wir betrachten nun den Fall, dass eine mehrfache Nullstelle mit dem Newton-Verfahren berechnet werden soll. Wir beschranken uns auf den Fall einer zweifa-chen Nullstelle x0, d.h. f(x0) = f ′(x0) = 0 und f ′′(x0) 6= 0. Wir betrachten eineUmgebung von x0 in der |f ′′(x)| ≥ m > 0 gilt. Fur die Newton-Iteration gilt nachdem Mittelwertsatz
xn+1 = xn −f(xn)− f(x0)
f ′(xn)− f ′(x0)= xn −
f ′(ξn)
f ′′(ηn)
fur geeignete Zwischenpunkte ξn, ηn. Der Quotient f(xn)/f ′(xn) bleibt also furxn gegen x0 wohldefiniert. Hinsichtlich der quadratischen Konvergenz ist aber zubeachten, dass wegen (Nachweis durch partielle Integration!)
f(x) = (x− x0)2
∫ 1
0
(1− t)f ′′(x0 + t(x− x0)) dt = (x− x0)2Q(x, x0)
f ′(x) = 2(x− x0)Q(x, x0) + (x− x0)2Q′(x, x0)
der Quotient f(x)/f ′(x) sich verhalt wie
f(x)
f ′(x)=
(x− x0)2Q(x, x0)
2(x− x0)Q(x, x0) + (x− x0)2Q′(x, x0)
=(x− x0)
2− (x− x0)
2
Q′(x, x0)
2Q(x, x0) + (x− x0)Q′(x, x0)
,
wobei
limx→x0
Q(x, x0) =f ′′(x0)
2, lim
x→x0
Q′(x, x0) =f ′′′(x0)
6.
Quadratische Konvergenz kann durch die modifizierte Iteration
xn+1 = xn −2f(xn)
f ′(xn), n ≥ 1,
erreicht werden. In der Tat gilt fur die modifizierte Iteration
xn+1 − x0 = xn − x0 −2f(xn)
f ′(xn)
= (x− x0)2
Q′(x, x0)
2Q(x, x0) + (x− x0)Q′(x, x0)
|xn+1 − x0| ≤ C |xn − x0|2.
Das Newton-Verfahren benotigt in jedem Iterationsschritt die Auswertungder ersten Ableitung f ′(xn), was bei komplizierten (moglicherweise nur implizit
90 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
durch ein Rechenprogramm definierter) Funktionen viel Aufwand erfordert. Inderartigen Fallen geht man zum vereinfachten Newton-Verfahren
xn+1 = xn −f(xn)
f ′(s), n ≥ 1
mit einem fest gewahlten Punkt s uber. Das Verfahren ist ein Spezialfall desallgemeineren Fixpunktverfahrens, σ 6= 0,
xn+1 = xn + σf(xn), n ≥ 1
zur Berechnung einer Nullstelle von f (oder aquivalent eines Fixpunktes der Ab-bildung P (x) = x + σf(x)). Konvergenzaussagen folgen aus dem BanachschenFixpunktsatz.
6.2 Konvergenzverhalten iterativer Verfahren
Das Newton-Verfahren zur Bestimmung einer einfachen Nullstelle besitzt lokaldie Eigenschaft
|xn+1 − x0| ≤ C|xn − x0|2.Man nennt es deshalb quadratisch oder von zweiter Ordnung konvergent. Allge-mein konvergiert ein Iterationsverfahren zur Berechnung einer Große x0 mit derOrdnung p, wenn es eine Konstante C > 0 gibt mit
|xn+1 − x0| ≤ C|xn − x0|p.
Im Falle p > 1 impliziert diee Abschatzung die Konvergenz des Verfahrens furStartwerte x1, die genugend nahe an x0 liegen (analoge Argumentationskette,wie beim Newton-Verfahren!). Im Falle linearer Konvergenz, d.h. p = 1, heißt diebestmogliche Konstante lineare Konvergenzrate und wir haben Konvergenz imFalle C < 1:
|xn+1 − x0| ≤ C|xn − x0| ≤ · · · ≤ cn|x1 − x0| → 0 (n→∞).
Gilt die Abschatzung|xn+1 − x0| ≤ Cn|xn − x0|
mit einer Nullfolge cn → 0 (n → ∞), so nennt man das Verfahren superlinearkonvergent. Bei Fixpunktiterationen xn+1 = P (xn) mit stetig differenzierbarerAbbildung P gilt
∣∣∣∣xn+1 − x0
xn − x0
∣∣∣∣ =
∣∣∣∣P (xn)− P (x0)
xn − x0
∣∣∣∣→ |P ′(x0)| (n→∞),
die lineare Konvergenzrate ist asymptotisch also gerade |P ′(x0)|, fur |P ′(x0)| = 0liegt mindestens superlineare Konvergenz vor.
6.3. INTERPOLATIONSVERFAHREN 91
Theorem 6.2.1 Die Funktion P sei in Umgebung des Fixpunktes x0 p-mal stetigdifferenzierbar mit p ≥ 2. Die Fixpunktiteration xn+1 = P (xn) hat genau danndie Ordnung p, wenn
P ′(x0) = P ′′(x0) = · · · = P (p−1)(x0) = 0 und P (p)(x0) 6= 0.
Beweis. Anwendung der Taylorschen Formel.
In Anwendung des Satzes betrachten wir das Newton-Verfahren zur Bestim-mung einer einfachen Nullstelle. Mit P (x) = x− f(x)/f ′(x) folgt im allgemeinen
P ′(x0) = 1− f ′(x0)2 − f(x0)f
′′(x0)
f ′(x0)2= 0 und P ′′(x0) =
f ′′(x0)
f ′(x0)6= 0,
es ist also von zweiter Ordnung.
6.3 Interpolationsverfahren
Ziel dieser Klasse von Verfahren ist die Nullstellenbestimmung ohne Auswertungvon Ableitungen, die jedoch effizienter als das Intervallschachtelungs-Verfahrenbzw. die einfache sukzessive Approximation sind.
Die Sekantenmethode berechnet ausgehend von einem Wertepaar (xn−1, xn)die neue Iterierte als Nullstelle der Geraden durch die Punkte (xn−1, f(xn−1))und (xn, f(xn)). Dies fuhrt zur Iteration:
xn+1 = xn − f(xn)xn − xn−1
f(xn)− f(xn−1), n ≥ 1.
Bei diesem Verfahren handelt es sich um ein 2-Schrittverfahren, da zur Berech-nung die beiden vorangegangenen Iterierten benotigt (und damit auch gespei-chert) werden. Ahnlich wie fur das Newton-Verfahren haben wir auch eine lokaleKonvergenzaussage fur die Sekantenmethode. Hierzu benotigt man die durch
γ0 = γ1 = 1, γn+1 = γn − γn−1, n ≥ 1
definierten Fibonacci-Zahlen.
Theorem 6.3.1 Die Funktion f ∈ C2[a, b] habe im Innern des Intervalls [a, b]eine Nullstelle x0 und es seien
m := minx∈[a,b]
|f ′(x)| > 0, M := maxx∈[a,b]
|f ′′(x)| <∞.
Sei ferner ρ > 0 klein genug gewahlt, so dass q := ρM/(2m) < 1 gilt und dieabgeschlossene Kugel B(x0, ρ) Teilmenge von [a, b] ist. Dann sind fur jedes Paar
92 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
von Startwerten x1, x2 ∈ B(x0, ρ) mit x1 6= x2 die Iterierten xn ∈ B(x0, ρ) derSekantenmethode wohl definiert und konvergieren gegen die Nullstelle x0. Dabeigelten die a-priori Fehlerabschatzung
|xn+1 − x0| ≤2m
Mqγn , n ≥ 1,
und die a-posteriori Fehlerabschatzung
|xn+1 − x0| ≤1
m|f(xn+1)| ≤
M
2m|xn+1 − xn| |xn+1 − xn−1|, n ≥ 1.
Beweis. Wie im Beweis zum Newton-Verfahren haben wir zunachst als Folgerungaus dem Mittelwertsatz fur zwei Punkte x, y ∈ [a, b] mit x 6= y
|x− y| ≤ 1
m|f(x)− f(y)|,
woraus die Eindeutigkeit der Nullstelle im Intervall [a, b] folgt. Weiter ist
f(x)− f(y)
x− y =
∫ 1
0
f ′(x+ t(y − x)) dt.
Mit einem dritten Punkt z ∈ [a, b], z 6= x, ergibt sich
f(x)− f(y)
x− y − f(x)− f(z)
x− z =
∫ 1
0
(f ′(x+ t(y − x))− f ′(x+ t(z − x))) dt.
Wegen
f ′(x+ t(y − x))− f ′(x+ t(z − x)) = (y − z)∫ t
0
f ′′(x+ t(y − x) + s(z − y)) ds
haben wir fur t ∈ (0, 1)
|f ′(x+ t(y − x))− f ′(x+ t(z − x))| ≤ M |y − z|t,
woraus ∣∣∣∣f(x)− f(y)
x− y − f(x)− f(z)
x− z
∣∣∣∣ ≤M
2|y − z|
folgt. Fur Punkte x, y ∈ B(x0, ρ) mit x 6= y, x 6= x0, y 6= x0 definieren wir
P (x, y) := x− f(x)x− y
f(x)− f(y).
6.3. INTERPOLATIONSVERFAHREN 93
Dann gilt
P (x, y)− x0 = x− x0 − f(x)x− y
f(x)− f(y)
=x− y
f(x)− f(y)
(x− x0)
f(x)− f(y)
x− y − f(x) + f(x0)
|P (x, y)− x0| ≤1
m|x− x0|
∣∣∣∣f(x)− f(y)
x− y − f(x)− f(x0)
x− x0
∣∣∣∣
≤ M
2m|x− x0| |y − x0| ≤
M
2mρ2 < ρ.
Die Iterierten der Sekantenmethode xn bleiben also in der Menge B(x0, ρ) undes gilt die Abschatzung
|xn+1 − x0| ≤M
2m|xn − x0| |xn−1 − x0|.
Setzt man ρn := |xn − x0|M/(2m), so folgt
ρn+1 ≤ ρnρn−1, n ≥ 2,
d.h. mit ρ1 ≤ q, ρ2 ≤ q gilt ρn+1 ≤ qγn. Wegen γn → ∞ (n → ∞) und q < 1konvergiert damit die Folge und wir haben die a-priori Abschatzung
|xn+1 − x0| =2m
Mρn+1 ≤
2m
Mqγn → 0 (n→∞).
Zum Nachweis der a-posteriori Fehlerabschatzung setzen wir oben x = xn, y =xn+1 und z = xn−1 und erhalten
|xn+1 − x0| ≤1
m|f(xn+1)− f(x0)| =
1
m|f(xn+1)|
≤ 1
m
∣∣∣∣f(xn) + (xn+1 − xn)f(xn+1)− f(xn)
xn+1 − xn
∣∣∣∣
≤ 1
m|xn+1 − xn|
∣∣∣∣f(xn+1)− f(xn)
xn+1 − xn− f(xn)− f(xn−1)
xn − xn−1
∣∣∣∣
≤ M
2m|xn+1 − xn| |xn − xn−1|
durch Ausnutzung der Bildungsvorschrift fur xn+1.
Um die Konvergenzgeschwindigkeit der Sekantenmethode zu beurteilen, benotigenwir Informationen uber das Anwachsen der Fibonacci-Zahlen γn fur n→∞. Die-se genugen der homogenen Differenzengleichung
γn+2 − γn+1 − γn = 0, n ≥ 0.
94 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
Zur Losung machen wir den Ansatz γn = λn und erhalten
λn(λ2 − λ− 1) = 0
zur Bestimmung von λ. Mit den Wurzeln λ1,2 = (1 ±√
5)/2 der quadratischenGleichung λ2 − λ− 1 = 0 lautet die allgemeine Losung der Differenzengleichung
γn = c1λn1 + c2λ
n2 .
Die unbekannten Koeffizienten bestimmen wir aus den Forderungen γ0 = γ1 = 1.Die Losung des zugeordneten Gleichungssystems und Einsetzen in die allgemeineDarstellung der Losung ergibt
γn =1√5
(λn+1
1 − λn+12
)=λn+1
1√5
(1−
(λ2
λ1
)n+1)∼ λn+1
1√5
(n→∞).
Die Sekantenmethode konvergiert also mindestens so schnell wie ein 1-Schrittverfahrender Ordnung p = 1.6. Fasst man jedoch zwei Schritte des Sekantenverfahrens zueinem Makro-Schritt zusammen, so erhalt man wegen
|x2n+1 − x0| ≤2m
Mq2γn , γ2n ∼
1√5λ2n
1 =1√5(λ1 + 1)n
ein Verfahren der Ordnung p ≥ 2.6. Da ein Schritt im Newton-Verfahren zweiFunktionsauswertungen erfordert und zwei Schritte des Sekantenverfahrens eben-so, ist das Sekantenverfahren asymptotisch bei gleichem Aufwand schneller. Demtheoretischen Vorteil steht aber die Gefahr von Ausloschungseffekten im Sekan-tenschritt bei monotoner Konvergenz von f(xn)→ 0 gegenuber. Man stabilisiertdie Methode deshalb auch mit der Intervallschachtelungsidee zur regula falsi:
Seien an < bn derart bestimmt, dass f(an)f(bn) < 0, d.h. f hat eineNullstelle x0 ∈ (an, bn). Beim Sekantenschritt
xn = an − f(an)an − bn
f(an)− f(bn)
tritt nun keine Ausloschung im Term f(an) − f(bn) auf, solange dieIntervalllange bn−an ≫ eps. Das neue Intervall [an+1, bn+1] bestimmtsich wie beim Intervallschachtelungsverfahren durch die Vorschrift:
f(xn) = 0 ⇒ STOP
f(an)f(xn) < 0 ⇒ an+1 = an, bn+1 = xn
f(an)f(xn) > 0 ⇒ an+1 = xn, bn+1 = bn.
6.4. NEWTON-VERFAHREN IM RD 95
Die hohere Stabilitat der regula falsi gegenuber dem Sekantenverfahrens wirddurch eine geringere Konvergenzgeschwindigkeit erkauft.
Die Methode der quadratischen Interpolation ist eine Weiterentwicklung derregula falsi. Die neue Iterierte wird als Nullstelle des quadratischen Interpolati-onspolynoms zu den Knoten an, bn und (an + bn)/2 bestimmt. Man kann zeigen,dass es im Fall f(an)f(bn) < 0 genau eine Nullstelle xn innerhalb von (an, bn) gibt,mit der dann wie beim Intervallschachtelungsverfahren weitergearbeitet wird.
6.4 Newton-Verfahren im Rd
Wir betrachten das Newton-Verfahren zur Losung nichtlinearer Gleichungssys-teme der Form f(x) = 0 mit stetig differenzierbarer Abbildung f : Ω ⊂ R
d →Rd. Formal erhalten wir das Newton-Verfahren durch Approximation der Glei-chung f(x) = 0 durch das lineare Taylor-Polynom an einer bereits bekanntenNaherungsstelle xn:
f(x) ≈ f(xn) + f ′(xn)(x− xn) = 0.
Damit lautet die Iterationsvorschrift
xn+1 = xn − (f ′(xn))−1f(xn), n ≥ 1,
mit der Jacobi-Matrix f ′(xn) von f an der Stelle xn. Wir konnen mit dem Ba-nachschen Fixpunktsatz lokale Konvergenz fur Losungen x0 zeigen, fur die dieJacobi-Matrix nicht verschwindet. Bei der Anwendung hat man vor allem mitzwei Problemen zu tun:
(a) hoher Aufwand pro Iterationsschritt(Losung eines linearen Gleichungssystems) und
(b) Wahl geeigneter Startwerte.
Zur Uberwindung dieser Probleme verwendet man gegebenenfalls das vereinfachteNewton-Verfahren
xn+1 = xn − (f ′(s))−1f(xn), n ≥ 1,
mit einem geeigneten s ∈ Rd, zum Beispiel s = x1. In diesem Fall haben die injedem Schritt zu losenden Gleichungssysteme alle die gleiche Koeffizientenmatrixf ′(s) und mit Hilfe einer einmal berechneten LR-Zerlegung kann der Aufwandbetrachtlich reduziert werden. Zur Vergroßerung des Konvergenzbereiches bietetsich das gedampfte Newton-Verfahren an:
xn+1 = xn − ωn (f ′(xn))−1f(xn), n ≥ 1,
bei dem ωn ∈ (0, 1] zu Beginn zunachst klein gesetzt und im Verlaufe der Rech-nung nach einer geeigneten Dampfungsstrategie schrittweise bis auf ωn = 1 erhohtwird.
96 KAPITEL 6. NICHTLINEARE GLEICHUNGEN
Kapitel 7
Lineare Gleichungssysteme II
In der Praxis kommen haufig dunn besetzte Bandmatrizen hoher Dimension n≫106 (z.B. bei der Diskretisierung von Differentialgleichungen) vor. Fur derartgroße Gleichungssysteme ist das Gaußsche Eliminationsverfahren wenig geeignet,da durch das fill in bei der LR-Zerlegung ein sehr hoher Speicherplatzbedarfentsteht, der sogar die Große des Kernspeichers des Rechners ubertreffen kann.Die im folgenden betrachteten iterativen Verfahren benotigen zur naherungsweiseLosung des Gleichungssystems Ax = b nicht viel mehr Speicherplatz, als zurSpeicherung von A erforderlich ist. Außerdem benotigt eine Iteration deutlichweniger als n3/3 arithmetische Operationen.
7.1 Einzelschritt- und Gesamtschrittverfahren
Wir betrachten das Gleichungssystem Ax = b, das ausgeschrieben unter Hervor-hebung des Diagonalelementes
aiixi +n∑
j=1j 6=i
aijxj = bi, i = 1, . . . , n
lautet. Im Falle aii 6= 0 erhalten wir
xi =1
aii
bi −
n∑
j=1j 6=i
aijxj
, i = 1, . . . , n.
Das Gesamtschritt- oder Jacobi-Verfahren erzeugt nun aus einer Iterierten x(k)
eine neue Iterierte x(k+1) durch die Vorschrift
x(k+1)i =
1
aii
bi −
n∑
j=1j 6=i
aijx(k)j
, i = 1, . . . , n.
97
98 KAPITEL 7. LINEARE GLEICHUNGSSYSTEME II
Zum Zeitpunkt der Berechnung von x(k+1)i sind die neuen Komponenten x
(k+1)l
mit l < i berechnet und konnten bereits in die Rechnung eingehen. Diese Ideefuhrt zum Einzelschritt- oder Gauß-Seidel-Verfahren
x(k+1)i =
1
aii
bi −
n∑
j=1j<i
aijx(k+1)j −
n∑
j=1j>i
aijx(k)j
, i = 1, . . . , n.
Beide Iterationen x(k) 7→ x(k+1) konnen als Fixpunktiteration zur Losung des Glei-chungssystems Ax = b aufgefaßt werden. Zur kompakten Schreibweise zerlegenwir A in A = D + L+R mit
D =
a11 0
. . .
0 ann
L =
0 · · · · · · 0
a21. . .
......
. . .. . .
...an1 · · · an,n−1 0
R =
0 a12 · · · a1n...
. . .. . .
......
. . . an−1,n
0 · · · · · · 0
.
Das Jacobi-Verfahren kann damit in der Form
x(k+1) = D−1b− (L+R)x(k)
= −D−1(L+R)x(k) +D−1b
geschrieben werden. Entsprechend gilt fur das Gauß-Seidel-Verfahren
x(k+1) = D−1b− Lx(k+1) − Rx(k)
= −(D + L)−1Rx(k) + (D + L)−1b.
Beide Verfahren lassen sich also in der Form
x(k+1) = Bx(k) + c
mit den Iterationsmatrizen B = −D−1(L+R) und B = −(D+L)−1R schreiben.Konvergiert die Folge, so konvergiert sie gegen einen Fixpunkt der AbbildungP (x) = Bx+ c. Beim Jacobi- und Gauß-Seidel-Verfahren sind per KonstruktionFixpunkte von P Losungen des Gleichungssystems. Zur Konstruktion allgemeineriterativer Verfahren diesen Typs wahlt man eine einfach zu invertierende MatrixC und iteriert ausgehend von
Ax = b ⇔ Cx = Cx− Ax+ b ⇔ x = x+ C−1(b− Ax)
in der Formx(k+1) = (I − C−1A)x(k) + C−1b.
Fur den Fall C = I heißt die Iteration auch Richardson-Verfahren. Nach demBanachschen Fixpunktsatz ist ein hinreichendes Kriterium fur die Konvergenzdieser Iteration, dass
‖B‖ = ‖I − C−1A‖ < 1
7.1. EINZELSCHRITT- UND GESAMTSCHRITTVERFAHREN 99
fur irgendeine Matrixnorm ‖ · ‖ auf Rn×n. Die Gultigkeit dieser Beziehung ist beigegebener Iterationsmatrix B aber von der gewahlten Matrixnorm abhangig. ZurCharakterisierung verwendet man daher besser den Spektralradius der Iterations-matrix B
spr(B) = max|λ| : λ Eigenwert von B.Aus
|λ|‖x‖ = ‖λx‖ = ‖Bx‖ ≤ ‖B‖ ‖x‖, x 6= 0 Eigenvektor zu λ,
folgt, dass fur jede naturliche Matrixnorm ‖ · ‖
spr(B) ≤ ‖B‖,
fur symmetrische Matrizen gilt sogar
spr(B) = ‖B‖2 = supx∈Rn×n
‖Bx‖2‖x‖2
.
Im allgemeinen ist jedoch der Spektralradius keine Norm, wie das Beispiel
A =
[0 10 0
]6= 0 ⇒ spr(A) = 0
zeigt. Es gilt aber das folgende Resultat.
Theorem 7.1.1 Fur jede Matrix B ∈ Rn×n gibt es zu jedem ε > 0 eine naturliche
Matrizennorm, so dass
spr(B) ≤ ‖B‖ε ≤ spr(B) + ε.
Beweis. Siehe z.B. J. Stoer, R. Bulirsch: Numerische Mathematik 2. Springer-Verlag 1990.
Theorem 7.1.2 Die Matrizen A ∈ Rn×n und C ∈ Rn×n seien regular. Die durch
x(k+1) = (I − C−1A)x(k) + C−1b
erzeugten Iterierten x(k) ∈ Rn konvergieren genau dann fur jeden Startwert x(1) ∈Rn gegen die Losung x0 ∈ Rn des Gleichungssystems Ax = b, wenn fur denSpektralradius der Iterationsmatrix B = I − C−1A die Beziehung spr (B) < 1gilt.
Beweis. Wir bezeichnen die Fehlervektoren durch e(k) := x(k) − x0. Da dieLosung x0 ∈ Rn Fixpunkt von P mit P (x) = Bx + c ist, folgt fur die Folge derFehlervektoren
e(k+1) = x(k+1) − x0 = Bx(k) + c− Bx0 − c = B(x(k) − x0) = Be(k).
100 KAPITEL 7. LINEARE GLEICHUNGSSYSTEME II
Damit haben wir e(k+1) = Bke(1). Im Fall spr (B) < 1 gibt es ein ε > 0 mitspr (B) + ε < 1 und eine naturliche Matrizennorm ‖ · ‖ε, so dass
‖e(k+1)‖ε = ‖Bke(1)‖ε ≤ ‖Bk‖ε‖e(1)‖ε ≤ ‖B‖kε‖e(1)‖ε → 0.
Da alle Normen im Rn aquivalent sind, konvergiert die Folge x(k) gegen x0.
Sei umgekehrt die Iteration konvergent fur jeden Startwert. Wir setzen alsStartwert x(1) = w + x0, wobei w ∈ Rn\0 Eigenvektor zum betragsmaßiggroßten Eigenwert λ von B ist. Aus
λkw = Bkw = Bk(x(1) − x0) = e(k+1) → 0 (k →∞)
folgt notwendig |λ| < 1 fur jeden Eigenwert von B, d.h. spr (B) < 1.
Bei der Konstruktion von Iterationsverfahren, d.h. der Wahl einer geeignetenMatrix C, sind zwei Ziele zu berucksichtigen:
(a) spr (I − C−1A) soll moglichst klein sein und
(b) die Gleichungssysteme Cx(k+1) = (C−A)x(k) +b sollen moglichst leicht undmit wenig zusatzlichem Speicherplatzbedarf losbar sein.
Ein Extremfall ist C = A mit spr (I − C−1A) = 0, aber die Losung der Glei-chungssysteme unter (b) entspricht der Losung des Ausgangssystems. Wir wer-den also gewisse Kompromisse eingehen mussen. Beim Jacobi- bzw. Gauß-Seidel-Verfahren lauten die Iterationsmatrizen J = −D−1(L + R) bzw. H = −(D +L)−1R, es sind also nur Dreieckssysteme zu losen. Hinsichtlich der Konvergenzgilt folgendes Resultat.
Theorem 7.1.3 Die Matrix A ∈ Rn×n sei strikt diagonal dominant, d.h.
|aii| >n∑
j=1j 6=i
|aij |.
Dann ist spr J < 1 und spr H < 1, d.h. das Jacobi- und das Gauß-Seidel-Verfahren konvergieren.
Beweis. Sei v Eigenvektor zum Eigenwert λ der Iterationsmatrix J des Jacobi-Verfahrens, d.h.
λv = Jv = −D−1(L+R)v.
In der Maximumnorm ‖ · ‖∞ folgt fur ‖v‖∞ = 1
|λ| = ‖λv‖∞ ≤ ‖D−1(L+R)‖∞‖v‖∞ = ‖D−1(L+R)‖∞
= max1≤i≤n
1
|aii|
n∑
j=1j 6=i
|aij |
< 1.
7.1. EINZELSCHRITT- UND GESAMTSCHRITTVERFAHREN 101
Ahnlich erhalten wir fur den Eigenvektor w zum Eigenwert µ der IterationsmatrixH des Gauß-Seidel-Verfahrens
µw = Hw = −(D + L)−1Rw ⇔ µw = −D−1(µL+R)w,
woraus fur ‖w‖∞ = 1
|µ| = ‖µw‖∞ ≤ ‖D−1(µL+R)‖∞‖w‖∞ = ‖D−1(µL+R)‖∞
= max1≤i≤n
1
|aii|
n∑
j=1j<i
|µ| |aij|+n∑
j=1j>i
|aij |
folgt. Angenommen 1 ≤ |µ|, dann folgt der Widerspruch
|µ| ≤ |µ|‖D−1(L+R)‖∞ < |µ|.
Also kann nur |µ| < 1 gelten.
In der Praxis vorkommende große Matrizen besitzen leider Spektralradiennahe bei 1 und konvergieren daher viel zu langsam. Man versucht daher die Kon-vergenz durch Einfuhrung von Relaxationsparametern zu beschleunigen. BeimSOR-Verfahren (Successive-OverRelaxation method) berechnet man im k-ten Ite-rationsschritt ausgehend vom Gauß-Seidel-Zwischenwert
x(k+1)i =
1
aii
bi −
n∑
j=1j<i
aijx(k+1)j −
n∑
j=1j>i
aijx(k)j
, i = 1, . . . , n
den nachsten Wert x(k+1)i als Linearkombination
x(k+1)i = ωx
(k+1)i + (1− ω)x
(k)i .
mit einem positiven Relaxationsparameter ω. Im Falle ω = 1 erhalten wir geradedas Gauß-Seidel-Verfahren. Ist ω ∈ (0, 1) spricht man von Unterrelaxation, furω > 1 von Uberrelaxation. Die Iterationsmatrix Hω des Relaxationsverfahrenserhalt man aus der Beziehung
x(k+1) = ωD−1b− Lx(k+1) − Rx(k)+ (1− ω)x(k)
in der Darstellung
Hω = (D + ωL)−1 [(1− ω)D − ωR] .
Der folgende Satz zeigt, dass das SOR-Verfahren hochstens fur Werte umω = 1 konvergieren kann.
102 KAPITEL 7. LINEARE GLEICHUNGSSYSTEME II
Theorem 7.1.4 Notwendig fur die Konvergenz des SOR-Verfahrens ist die Be-dingung 0 < ω < 2.
Beweis. Seien λ1, . . . , λn die Eigenwerte der Iterationsmatrix
Hω = (D + ωL)−1[(1− ω)D − ωR].
Die Iterationsmatrix ist das Produkt zweier Dreiecksmatrizen, daher gilt
det Hω = det [(D + ωL)−1] · det [(1− ω)D − ωR]
=
n∏
i=1
1
aii·
n∏
j=1
(1− ω)ajj = (1− ω)n.
Andererseits ist die Determinante einer Matrix das Produkt ihrer Eigenwerte,woraus
|1− ω|n =∣∣(1− ω)n
∣∣ =
∣∣∣∣∣
n∏
i=1
λi
∣∣∣∣∣ =n∏
i=1
|λi| ≤ (spr(Hω))n
folgt. Nun ist |1− ω| < 1 aquivalent zu 0 < ω < 2.
Theorem 7.1.5 Fur eine symmetrische, positiv definite Matrix A ∈ Rn×n gilt
spr (Hω) < 1 fur 0 < ω < 2.
Beweis. Wegen der Symmetrie von A haben wir R = LT , d.h.
A = L+D + LT .
Sei λ Eigenwert von Hω mit zugehorigem Eigenvektor v ∈ Rn. Es gilt
Hωv = λv ⇔[(1− ω)D − ωLT
]v = λ(D + ωL)v
und damit auchω(D + LT )v = (1− λ)Dv − λωLv.
Hieraus berechnen wir
ωAv = ω(D + LT )v + ωLv = (1− λ)Dv + (1− λ)ωLv
und
λωAv = λω(D + LT )v + λωLv
= λω(D + LT )v + (1− λ)Dv − ω(D + LT )v
= (λ− 1)ω(D + LT )v + (1− λ)Dv
= (1− λ)(1− ω)Dv − (1− λ)ωLTv.
7.2. ABSTIEGSVERFAHREN 103
Multiplikation beider Gleichungen mit vT und Addition ergibt
(1 + λ)ωvTAv = (1− λ)(2− ω)vTDv + ω(1− λ)(vTLv − vTLTv)
= (1− λ)(2− ω)vTDv,
woraus wegen der positiven Definitheit von A und D folgt, dass λ 6= ±1 unddamit
µ :=1 + λ
1− λ =2− ωω
vTDv
vTAv> 0.
Nach λ umgestellt, sehen wir |λ| =∣∣(µ− 1)/(µ+ 1)
∣∣ < 1.
Die Frage nach optimalen Relaxationsparametern ist nicht leicht zu beant-worten. Wir betrachten hier als Beispiel nur das Richardson-Verfahren
x(k+1) = x(k) + (b− Ax(k)), x(k+1) = ωx(k+1) + (1− ω)x(k)
mit der Iterationsmatrix B = I − ωA. Wir setzen A symmetrisch und positivdefinit voraus. Dann ist auch die Iterationsmatrix symmetrisch. Ist λ Eigenwertvon A mit zugeordnetem Eigenvektor v ∈ R
n, so ist 1−ωλ Eigenwert von I−ωA,denn
(I − ωA)v = v − ωλv = (1− ωλ)v.
Fur ω 6= 0 gilt auch die Umkehrung. Der Spektralradius der Iterationsmatrixwird damit
spr (I − ωA) = max (|1− ωλmin(A)|, |1− ωλmax(A)|) .
Minimierung von spr (I − ωA) uber ω fuhrt auf ωopt = 2/(λmin(A) + λmin(A)).
7.2 Abstiegsverfahren
Wir betrachten in diesem Abschnitt ausschließlich symmetrische, positiv definiteMatrizen A ∈ Rn×n. Dann sind alle Eigenwerte von A positiv, d.h.
0 < λmin ≤ λi ≤ λmax,
und durch‖x‖A :=
√(Ax, x), x ∈ R
n
wird eine Norm definiert, die wir A−Norm nennen wollen.
Aus der Analysis wissen wir, dass die notwendige Bedingung fur ein lokalesMinimum einer skalaren Funktion mehrerer Veranderlicher Q : Rn → R (odereines Funktionals) das Verschwinden der ersten partiellen Ableitungen ist. Wir
104 KAPITEL 7. LINEARE GLEICHUNGSSYSTEME II
suchen nun ein Funktional Q : Rn → R, so dass die notwendige Bedingung furMinima gerade der Aufgabe Ax = b entspricht. Sei
Q(y) :=1
2(Ay, y)− (b, y).
Theorem 7.2.1 Sei A s.p.d. Dann ist x Losung von Ax = b genau dann, wenn
Q(x) = miny∈Rn
Q(y).
Ferner ist das Minumum eindeutig bestimmt, also Q(x) < Q(y) ∀y 6= x.
Beweis. Notwendig fur ein Minimum von Q in x ∈ Rn ist das Verschwinden derRichtungsableitung in einer beliebigen Richtung r ∈ Rn, ‖r‖ = 1, d.h.
∂Q
∂r= lim
t→0
Q(x+ tr)−Q(x)
t=
1
2limt→0
[(Ax, r) + (Ar, x)− 2(b, r) + t(Ar, r)] = 0.
Mit der Symmetrie der Matrix A folgt hieraus (Ax− b, r) = 0 fur jede Richtung,also ist jede Minimalstelle von Q Losung des Gleichungssystems. Sei nun x dieLosung des Gleichungssystems Ax = b. Dann gilt aufgrund der Symmetrie von A
‖y − x‖2A = (A(y − x), y − x) = (Ay, y)− (Ax, y)− (Ay, x) + (Ax, x)
= (Ay, y)− 2(Ax, y)− (Ax, x) + 2(Ax, x)
= (Ay, y)− 2(b, y)− (Ax, x)− 2(b, x)= 2Q(y)− 2Q(x).
Somit ist Q(y) > Q(x) fur alle y 6= x.
Die Richtungsableitung von Q (siehe Beweis oben) wird maximal fur denGradienten
grad Q(y) = Ay − b.Sei r(k) eine zunachst beliebige Richtung. Wir wollen in Richtung von r(k) fort-schreiten und suchen den bestmoglichen Abstieg, also x(k+1) = x(k) + αkr
(k), sodass (line search)
Q(x(k+1)) = minα∈R
Q(x(k) + αr(k)).
Fur das optimale α ergibt sich
αk = −(Ax(k) − b, r(k))
(Ar(k), r(k)),
das berechenbar ist, solange r(k) 6= 0. Die Idee des Gradienten-Verfahrens bestehtnun darin, als Richtung den steilsten Abstieg im Punkt x(k) zu wahlen, also
r(k) = − grad Q(x(k)) = −(Ax(k) − b),
7.2. ABSTIEGSVERFAHREN 105
und in dieser Richtung bestmoglich abzusteigen. Das Gradienten-Verfahrenlautet demnach:
Start: x(1), g(1) := Ax(1) − bk→ k + 1 : w(k) = Ag(k)
αk = − (g(k), g(k))
(w(k), g(k))
x(k+1) = x(k) − αkg(k)
g(k+1) = g(k) − αkw(k)
Zur Durchfuhrung des Verfahrens benotigt man also die drei Vektoren x, g undw. Man beachte, dass fur g(k) 6= 0 der Nenner in der Berechnung von αk wegen
(w(k), g(k)) = (Ag(k), g(k)) > 0
nicht verschwindet. Hinsichtlich der Konvergenz des Verfahrens gilt
Theorem 7.2.2 Sei A s.p.d. und κ := cond2(A) = λmax/λmin. Dann gilt
‖x(k+1) − x‖A ≤(κ− 1
κ+ 1
)k
‖x(1) − x‖A k ≥ 1,
d.h. das Gradientenverfahren konvergiert fur beliebigen Startwert x(1) gegen dieLosung x des Gleichungssystems Ax = b.
Ist κ nahe bei Eins, die Eigenwerte also nahe beieinander, konvergiert das Gradi-entenverfahren ziemlich schnell. Fur schlecht konditionierte Matrizen konvergiertim allgemeinen das Gradientenverfahren sehr langsam.
Das Verfahren der konjugierten Gradienten (conjugate gradient method) oderkurz CG-Verfahren von Hestenes und Stiefel (1952) zielt auf eine Verbesse-rung des Gradientenverfahrens durch Modifikation der Abstiegsrichtungen. Hier-zu uberlegen wir uns zunachst, wann ein Punkt x optimal bezuglich einer Such-richtung d ist, d.h. wann
Q(x) = minα∈R
Q(x+ αd).
Anwendung der Extremwerttheorie ergibt, dass x optimal bezuglich der Such-richtung d genau dann ist, wenn (g, d) = 0 mit g = Ax − b gilt. Wir sagen, dassx optimal bezuglich des Unterraums Bk := span (d(1), . . . , d(k)) ist, wenn
Q(x) = miny∈x+Bk
Q(y) = mind∈Bk
Q(x+ d).
106 KAPITEL 7. LINEARE GLEICHUNGSSYSTEME II
Obige Argumentation zeigt, dass dies aquivalent zu
(g, d(i)) = 0 i = 1, . . . , k, g = Ax− b,
ist. Sei x optimal bezuglich Bk. Wir suchen eine Richtung d, so dass x′ = x + doptimal bezuglich Bk bleibt. Die Beziehung
0 = (A(x+ d)− b, d(i)) = (Ax− b, d(i)) + (Ad, d(i)) = (Ad, d(i))
zeigt, dass notwendig und hinreichend ist, dass d konjugiert (oder A-orthogonal)zu allen d(i) sein muss, also
(d, d(i))A = (Ad, d(i)) = 0 i = 1, . . . , k.
Die Idee des CG-Verfahrens besteht nun darin, anstelle des steilsten Abstiegs
−g(k+1) = b− Ax(k+1)
eine neue Richtung
d(k+1) ⊥A Bk = span (d(1), . . . , d(k))
zu wahlen. Dazu setzen wir mit unbekannten Koeffizienten β(k+1)j
d(k+1) = −g(k+1) +
k∑
j=1
β(k+1)j d(j)
an und erhalten aus der Forderung der A-Orthogonalitat von d(k+1) zu d(i)
0 = (d(k+1), d(i))A = −(g(k+1), d(i))A + β(k+1)i (d(i), d(i))A
beziehungsweise
β(k+1)i =
(g(k+1), Ad(i))
(d(i), Ad(i))i = 1, . . . , k.
Es stellt sich heraus, dass bei dieser Wahl
Bk = span (d(1), . . . , d(k)) = Kk(d(1);A) := span (d(1), Ad(1), . . . , Ak−1d(1)),
wobei Kk Krylov-Raum genannt wird. Ferner gilt
β(k+1)i = 0 i < k, β
(k+1)k 6= 0 falls g(k+1) 6= 0.
Wie beim Gradienten-Verfahren erhalt man aus der neuen Suchrichtung d(k) dieneue Iterierte x(k+1) sowie g(k+1) = Ax(k+1)−b. Die neue Iterierte ist dann optimalbezuglich d(k+1) und Bk, d.h. bezuglich Bk+1:
Q(x(k+1)) = mind∈Bk+1
Q(x(k+1) + d) = miny∈x(1)+Kk+1
Q(y).
7.2. ABSTIEGSVERFAHREN 107
Das CG-Verfahren lautet damit:
Start: x(1), g(1) = Ax(1) − b, d(1) = −g(1)
k→ k + 1 : w(k) = Ad(k)
γk = (g(k), g(k)) αk = γk/(d(k), w(k))
x(k+1) = x(k) + αkd(k) g(k+1) = g(k) + αkw
(k)
βk = (g(k+1), g(k+1))/γk d(k+1) = −g(k) + βkd(k)
Man braucht also 4 Vektoren x, g, d und w. Da die Richtungen d(i) paarwei-se A-orthogonal sind, sind sie linear unabhangig (solange das Verfahren (nichtabbricht); damit gilt spatestens nach n Schritten
span (d(1), . . . , d(n)) = Rn.
und bei exakter Arithmetik bricht das CG-Verfahren spatestens nach n Schrittenmit der exakten Losung ab. Hinsichtlich der Konvergenz gilt
Theorem 7.2.3 Sei A s.p.d. und κ := cond2(A) = λmax/λmin. Dann gilt
‖x(k+1) − x‖A ≤ 2
(√κ− 1√κ+ 1
)k
‖x(1) − x‖A k ≥ 1,
d.h. das CG-Verfahren konvergiert fur beliebigen Startwert x(1) gegen die Losungx des Gleichungssystems Ax = b.
Beweis. Deuflhard, Hohmann: Numerische Mathematik I, 1993
top related