Numerisches Rechnen - IGPM€¦ · dsdas Linienintegral uber den Weg bezeichnet) = logI(kb ak)...
Transcript of Numerisches Rechnen - IGPM€¦ · dsdas Linienintegral uber den Weg bezeichnet) = logI(kb ak)...
Numerisches Rechnen
Vorlage zur Vorlesung im WS 2009/2010
Prof. Henning Esser
Layout und LATEX-Fassung: Andrea EmondsLATEX-Fassung des Anhangs: Dagmar RostekKorrektur: Dr. Normann Pankratz
Dipl. Math. Christian PleskenUberarbeitung: Dipl.-Math. Dipl.-Phys. Kolja Brix
Institut fur Geometrie und Praktische Mathematikan der Rheinisch-Westfalischen Technischen Hochschule Aachen
Dies ist die Vorlage zur Vorlesung ”Numerisches Rechnen“ von Prof. Esser im SS2007 und kein Skriptum. Diese Vorlage soll den Studenten die Mitschrift bei derVorlesung erleichtern und ist kein Ersatz der Vorlesung. Wir weisen ausdrucklichdarauf hin, dass die vorliegende Fassung noch fehlerbehaftet sein kann. Daher istfur die Klausuren in diesem Semester einzig die Vorlesung als Grundlage anzusehen,Berufungen auf diesen Text sind nicht moglich. Bitte weisen Sie uns auf eventuelleFehler hin ([email protected]), damit diese Vorlage zu einem Skriptausgearbeitet werden kann.
Nachdruck und Vervielfaltigung, auch in Teilen, ist nicht gestattet.Copyright 2007, IGPM RWTH Aachen
Inhaltsverzeichnis
1 Lineare Gleichungssysteme, LR- und Cholesky-Zerlegung,iterative Verfahren 2
1.1 LR- und Cholesky-Zerlegung . . . . . . . . . . . . . . . . . . 4
1.2 Iterative Verfahren . . . . . . . . . . . . . . . . . . . . . . . . 20
2 Approximation im quadratischen Mittel 26
2.1 Normalgleichungen . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2 Householder-Verfahren . . . . . . . . . . . . . . . . . . . . . . 33
3 Nichtlineare Gleichungssysteme 43
4 Interpolation, Integration und Differentiation 50
5 Gewohnliche und (wenige) partielle Differentialgleichungen,numerische Methoden 65
6 Optimierung 99
A Lineare Differentialgleichungen 112
Literatur 123
Index 126
1 Lineare Gleichungssysteme 2
1 Lineare Gleichungssysteme, LR- und Cholesky-Zerlegung, iterative Verfahren
Bei einem der wichtigsten bildgebenden Verfahren der Medizintechnik, derComputertomographie, ist das schnelle Losen großer linearer Gleichungs-systeme von zentraler Bedeutung: Wird z.B. ein Kopf untersucht, so werdenschichtweise Bilder erstellt (”Kopf wird geschichtet“), aus denen dann ein2-d oder auch 3-d Bild der inneren Strukturen erstellt werden kann.
Draufsicht einer Schicht(ohne Nase und Ohren)
ab
Γ
Ω
Auf jeder Schicht werden parallel zur Strecke [a, b] = Γ Rontgenstrahlen(oder ahnliches) durch den Kopf geschickt. Die Intensitat der Strahlung ander Stelle x (von Γ), I(x), genugt fur kleine Strecken ∆s (auf Γ) in ersterNaherung dem Gesetz
0
Γx
x+∆s n
n
‖n‖=1
(1)I(x+ n∆s)− I(x)
∆s= −µ(x)I(x)
mit dem Schwachungskoeffizienten µ(x) bzw. der Absorptionsfunktion. D.h.der durchschnittliche Strahlungsverlust pro Langeneinheit auf der Strecke[x, x+ n∆s] ist proportional zur Intensitat der Strahlung an der Stelle x.Gesucht ist naturlich die Funktion µ(x), die man direkt inneren Struktu-ren (Knochen, Tumor, ...) zuordnen kann. Messen kann man aber nur I(a)(gegeben) und I(b).
Wir wenden (1) auf kleine Teilstrecken von Γ an:
xi = a+ i∆s n, i = 0, . . . , N, xN = a+N ∆s n = b.
1 Lineare Gleichungssysteme 3
0
b
axi
xi+1
n
(1) ergibt
I(a+ (i+ 1) ∆s n)− I(a+ i∆s n) = −µ(a+ i∆s n)I(a+ i∆s n)∆s
oder mit der Abkurzung I(s) = I(a+ sn)
(2) I((i+ 1) ∆s)− I(i∆s) = −µ(a+ i∆s n)I(i∆s), i = 0, 1, . . . , N − 1
Wir nehmen I ∈ C1[0, ‖b− a‖] an. Dann ist
I((i+ 1) ∆s)− I(i∆s) = I ′(i∆s+ Oi ∆s)∆s = I ′(i∆s)∆s+O(∆s)∆s
und somit
I ′(i∆s)∆s = −µ(a+ i∆s n)I(i∆s)∆s+O(∆s)∆s.
Wegen I ≥ C > 0 folgt dann
d
dslog I(s)
∣∣∣∣s=i∆s
∆s = −µ(a+ i∆s n)∆s+O(∆s)∆s,
N−1∑i=0
d
dslog I(s)
∣∣∣∣s=i∆s
∆s = −N−1∑i=0
µ(a+ i∆s n)∆s+O(∆s).
Der Grenzubergang ”∆s→ 0“ (wie in der Analysis) ergibt daher∫ ‖b−a‖0
d
dslog I(s) ds = −
∫ ‖b−a‖0
µ(a+ sn) ds ≡ −∫
Γµ(x) ds
(wobei∫
Γ µds das Linienintegral uber den Weg Γ bezeichnet)
= log I(‖b− a‖)− log(I(0)).
Also haben wir das fundamentale Gesetz
(3) logI(a)I(b)
=∫
Γµ(x) ds.
Uberstreicht man mit Γ das Gebiet ”richtig“, so ist durch (3) eine sogenannteIntegralgleichung fur µ gegeben, die Radon-Transformation von µ.
1 Lineare Gleichungssysteme 4
Zur naherungsweisen Losung von (3) uberdecken wir Ω mit einem aquidi-stanten Gitter bestehend aus den Quadraten Qi, i = 1, 2, . . . , n.
Q1 Q2 Q3
Qn
. . .
Ω
Ist ϕi(x) =
1 x ∈ Qi0 sonst
, so machen wir fur µ den Naherungsansatz
(4) µn(x) =n∑i=1
xiϕi(x)(⇒ µn(x) = xi, x ∈
Qi
)mit den ”Grauwerten“ xi in Pixel i.Ferner mussen wir die Schar Γ diskretisieren, d.h. wir nehmen endlich vie-le Gruppen paralleler Geraden, die gegeneinander verdreht sind (s. Folie):Γj , j = 1, 2, . . . ,m, wobei i.a. die Zahl der Rontgenstrahlen m wesentlichgroßer als n ist, m n.Fuhrt man die Messungen bj = log I(aj)
I(bj)durch, so ergibt sich statt (3) fur
j = 1, . . . ,m:∫Γj
µn(x) ds = bj =n∑i=1
(∫Γj
ϕi(x) ds
)xi =:
n∑i=1
ajixi,
d.h. ein lineares Gleichungssystem
(5) Ax = b, A ∈ Rm×n mit i.a. m n.
1.1 LR- und Cholesky-Zerlegung
Fur ein solches System sind nur fur m = n die Losungsverhaltnisse uber-sichtlich (womit wir uns gleich befassen). Fur m > n wird (5) i.a. nichtlosbar sein. In diesem Fall werden wir ein x bestimmen, das wenigstens dasResiduum ‖Ax− b‖2 minimiert (siehe Kap. 2).
1 Lineare Gleichungssysteme 5
Wir befassen uns jetzt mit dem Fall m = n, d.h.
(1.1)a11x1 + a12x2 + . . . + a1nxn = b1
......
an1x1 + an2x2 + . . . + annxn = bn
Bekanntlich ist (1.1) bzw. Ax = b fur jedes b ∈ Rn (eindeutig) losbar genaudann, wenn detA 6= 0, x = A−1b. Man kann sogar eine explizite Formelangeben:
Ist A = (a1, . . . , an) =
a1...an
mit den Spalten- bzw. Zeilenvektoren ai bzw.
ai, so gilt (Cramersche Regel)
(1.2) xi =det(a1, . . . , ai−1, b, ai+1, . . . , an)
detA, i = 1, . . . , n.
Diese fur theoretische Zwecke sehr nutzliche Formel ist praktisch nicht an-wendbar, da die Zahl der Operationen wie n! wachst (eine untere Schrankefur die Zahl der Operationen, um (1.1) zu losen ist sicher n, da man x jahinschreiben muss). Fur ein kleines System mit n = 50 wurde eine heutigeMaschine mit 1/46 · 10−12 sec/Operation 2 · 1043 Jahre brauchen.Bei dem jetzt behandelten Gauß-Algorithmus (GA) ist die Zahl der Flops∼ n3, womit dieses Problem in 2.7 · 10−9 sec gelost wird.Der GA beruht auf der Tatsache, dass die Losungsverhaltnisse von (1.1)nicht verandert werden, wenn man ein Vielfaches der (z.B.) 1. Gleichungvon der (z.B.) 2. Gleichung subtrahiert:Schreibt man (1.1) in der Form
(1.3) ai · x = bi, i = 1, . . . , n,
(ai · x = (ati, x) mit dem Skalarprodukt ( · , · )) so haben (1.3) und
a1x = b1
a2x− αa1x = b2 − αb1(1.4’)aix = bi, i = 3, . . . , n,
bzw.
(1.4)a1x = b1
aix− αia1x = bi − αib1, i = 2, . . . , n
=: A(2)x = b(2)
dieselben Losungen (Beweis?).
A(2) entsteht aus A(1) = A durch Subtraktion von
0
α2a1...
αna1
=: g(1)A(1);
und es ist b(2) = b− g(1)b.
1 Lineare Gleichungssysteme 6
Offenbar lautet die Matrix g(1)
(1.5) g(1) =
0 0 · · · 0α2 0 · · · 0...
.... . .
...αn 0 · · · 0
(Es ist: (g(1)A(1))t = Atg(1)t = (0, α2a
t1, . . . , αna
t1)
= (at1, . . . , atn)
0 α2 . . . αn0 0 . . . 0...
.... . .
...0 0 . . . 0
⇒ (1.5)).
Es ist daher
(1.6)A(2) = A(1) − g(1)A(1) = (1− g(1))A(1) =: G(1)A(1)
b(2) = G(1)b(1) (b(1) = b)
mit der Gauß-Matrix
(1.7) G(1) = 1− g(1).
Es ist detG(1) = 1.
Der erste Schritt des GA besteht darin, die αi in (1.4) so zu wahlen, dassdie Unbekannte x1 nur in der ersten Gleichung vorkommt:
(1.8) αi =ai1a11
=: li1, i = 2, . . . , n, im Fall a11 6= 0.
Es ist daher
(1.9) g(1) =
0 0 . . . 0l21 0 . . . 0...
.... . .
...ln1 0 . . . 0
, G(1) = 1− g(1) =
1−l21 1
.... . .
−ln1 1
und G(1)A(1) und G(1)b(1) notieren wir in der Form
(1.10)
a11 a12 . . . a1n b1
0 a(2)22 . . . a
(2)2n b
(2)2
......
. . ....
...0 a
(2)n2 . . . a
(2)nn b
(2)n
=
a11 a12 . . . a1n b1
0 b(2)2
... A(2)...
0 b(2)n .
Auf das (n− 1)× (n− 1) Restsystem (mit der Matrix A(2)) wenden wir denersten Schritt an, falls a(2)
22 6= 0, usw. .
1 Lineare Gleichungssysteme 7
Ist der Algorithmus durchfuhrbar, so erhalten wir nach k − 1 SchrittenA(k), b(k)
(1.11)
a(1)11 · · · b
(1)1
0. . . · · ·
... 0 a(k)kk · · · a
(k)kn
......
......
. . ....
0 0 a(k)nk · · · a
(k)nn b
(k)n
=
a(1)11 b
(1)1
. . ....
A(k)
b(k)n
Mit (s. (1.8))
(1.12) lik :=a
(k)ik
a(k)kk
, i = k + 1, . . . , n, falls a(k)kk 6= 0
g(k) =
0 · · · 0 0 0 · · · 0...
... 0...
... lk+1,k...
...0 · · · 0 ln,k 0 · · · 0
← k
, G(k) := 1− g(k)(1.13)
↑k
ist dann
(1.14’) A(k+1) = G(k)A(k), b(k+1) = G(k)b(k), k = 1, 2, . . . , n− 1, d.h.
Sei a(k)kk 6= 0
(1.14)
Fur i = k + 1, . . . , n :Fur j = k + 1, . . . , n :
a(k+1)ij = a
(k)ij −
a(k)ik
a(k)kk
akj
Ende fur
b(k+1)i = b
(k)i −
a(k)ik
a(k)kk
b(k)k
Ende fur
Es ist A(n) = R =
a(1)11 · · · a
(1)1n
. . ....
0 a(n)nn
eine obere Dreiecksmatrix , und das
System
(1.15) A(n)x = b(n)
1 Lineare Gleichungssysteme 8
ist einfach durch Ruckwartseinsetzen zu losen.
(1.16)
Fur i = n, n− 1, . . . , 1 :
xi =1
a(i)ii
b(i)i − n∑j=i+1
a(i)ij xj
Ende fur
Man kann die Großenordnung des Aufwands leicht abschatzen. Die Rech-nung (1.16) ist ∼ n2 und im k−ten Gauß-Schritt sind (n − k)2 Matrixele-mente zu berechnen.
∑n−1k=1(n− k)2 =
∑n−1k=1 k
2 ∼ n3.Also ist der Aufwand des GA ∼ n3 (genauer 1
3n3 + . . .)(n 1).
Ist der GA durchfuhrbar, so ist offenbar G(n−1) · · ·G(1)A = R und diesimpliziert eine wichtige Zerlegung der Matrix A, die LR-Zerlegung A = L ·Rmit einer normierten (lii = 1) linken unteren Dreiecksmatrix L.
Dazu benotigen wir einige Eigenschaften der Matrizen g(k), G(k).Sei B ∈ R(n×n) mit den Zeilen bi, i = 1, . . . , n. Ist
g(k) =
0 · · · 0 0 0 · · · 0...
... 0...
... αk+1,k...
...0 · · · 0 αn,k 0 · · · 0
← k
↑k
so ist offenbar g(k)B =
0...0
k
αk+1bk...
αnbk
und daher (1 + g(k))(1 − g(k))B = B,
d.h. (1 + g(k)) = (1− g(k))−1.
Lemma 1.1 Fur die Matrizen g(k), G(k) in (1.13) gilt
(i) g(k)g(l) = 0, k ≤ l
(ii) (1 + g(k)) = (1− g(k))−1
1 Lineare Gleichungssysteme 9
Beweis:
(i) Durch Inspektion 2
Es ist jetzt
A = (G(n−1) · · ·G(1))−1R = G(1)−1 · · ·G(n−1)−1R
= (1 + g(1))(1 + g(2)) · · · (1 + g(n−1))R,
und
(1 + g(1))(1 + g(2)) = 1 + g(1) + g(2) + g(1)g(2) = 1 + g(1) + g(2).
Also ist (1 + g(1)) · · · (1 + g(n−1)) = 1 +n−1∑i=1
g(i) =: L.
Damit haben wir
Satz 1.1 Der GA sei durchfuhrbar, d.h. a(k)kk 6= 0, k = 1, 2, . . . , n− 1. Dann
erzeugt der GA die LR-Zerlegung der Matrix A, A = L ·R, mit der rechtenoberen Dreiecksmatrix R
R = A(n) =
r11 · · · r1n
. . ....rnn
und der linken unteren Dreiecksmatrix L
(1.17) L = 1 +n−1∑i=1
g(i) =
1l21 1 0... l32
. . ....
.... . . . . .
ln1 ln2 · · · ln,n−1 1
.
Ist detA 6= 0, so ist diese Zerlegung eindeutig, wenn sie existiert. D.h. giltA = L′R′ = LR mit l′ii = lii = 1, so ist L = L′, R = R′. (Ubung)
(1.18) Beispiel:
2 1 1 74 1 0 6−2 2 1 5
1−g(1)=G(1)−→ 2 1 1 70 −1 −2 −80 3 2 12
1−g(2)=G(2)−→ 2 1 1 70 −1 −2 −80 0 −4 −12
g(1) =
0 0 02 0 0−1 0 0
, g(2) =
0 0 00 0 00 −3 0
1 Lineare Gleichungssysteme 10
⇒ x3 = 3⇒ x2 = 2⇒ x1 = 1 ; detA = 8 und
R =
2 1 10 −1 −20 0 −4
, L =
1 0 02 1 0−1 −3 1
.
Die Nullen im obigen Schema konnen wir bis auf lii = 1 durch die Eintragevon L ersetzen:
2 1 14 1 0−2 2 1
−→2 1 1
4/2 −1 −2−2/2 3 2
−→2 1 12 −1 −2−1 3/− 1 −4
Die Losung von Ax = b (fur beliebiges b) mit der LR-Zerlegung erfolgt so:
Ax = b⇔ LRx = b.
Sei Rx = y ⇒ Ly = b ⇒ y = L−1b ≡ b(n) (Vorwartseinsetzen); und manlost Rx = y durch Ruckwartseinsetzen.
Der GA ist jedoch in dieser Form nicht immer durchfuhrbar.
(1.19) Beispiel: Sei
A =
1 3 22 6 92 8 8
, B =
1 3 2 42 6 9 52 6 8 62 6 9 7
.
Es ist detA 6= 0, detB = 0 (s.u.).
1 3 22 6 92 8 8
G(1)
−→1 3 20 0 50 2 4
(⇒ detA = −10)
Abhilfe: Vertausche 2. und 3. Gleichung. D.h. hatte man die 2. und 3. Zeilevon A vertauscht, beschrieben durch PA, so gilt PA = LR! Wir kommendarauf zuruck.
1 3 2 42 6 9 52 6 8 62 6 9 7
G(1)
−→
1 3 2 40 0 5 10 0 4 −20 0 5 −1
A(2)
G(2)=1=
1 3 2 40 0 5 10 0 4 −20 0 5 −1
A(3)
G(3)
−→
1 3 2 40 0 5 10 0 4 −20 0 0 3
2
mit:
g(1) =
0 0 0 02 0 0 02 0 0 02 0 0 0
, g(2) = 0, g(3) =
0 0 0 00 0 0 00 0 0 00 0 5
4 0
,
1 Lineare Gleichungssysteme 11
R =
1 3 2 40 0 5 10 0 4 −20 0 0 3
2
, L =
1 0 0 02 1 0 02 0 1 02 0 5
4 1
.
Und die Losungen von Bx = b sind die Losungen von Rx = y, wobei y dieeindeutige Losung von Ly = b ist.
Es gilt daher: Der GA mit Zeilenvertauschung (ZV) angewendet auf Ax =b liefert nach (n−1) Schritten das aquivalente System A(n)x = b(n) mit eineroberen Dreiecksmatrix A(n) = R.
Bemerkung: Fur den GA mit ZV gilt:detA = 0 ⇔ fur mindestens ein k ∈ 1, . . . , n ist die 1. Spalte der auftre-tenden Restmatrix A(k) ∈ R(n−(k−1))×(n−(k−1)) die Nullspalte.
Es ist naturlich zu vermuten, dass der GA mit ZV eine LR-Zerlegung einerMatrix PA liefert, deren Zeilen eine Permutation der Zeilen von A sind. Wirgeben jetzt die Matrix P an, die ein Vertauschen der i-ten mit der k-ten Zei-le bewirkt.Mit den Einheitsvektoren ei, i = 1, . . . , n, gilt fur jede MatrixA = (a1, . . . , an)
A = 1 ·A =
(e1, a1) · · · (e1, an)... · · ·
...(ei, a1) · · · (ei, an)
... · · ·...
(ek, a1) · · · (ek, an)... · · ·
...(en, a1) · · · (en, an)
← i
← k
(1.20) P =
e1t
...ekt
...eit
...ent
← i
← k
PA bewirkt also das Vertauschen der Zeile i mit der Zeile k. P ist einePermutationsmatrix , fur die offenbar gilt P tP = 1 (orthogonale Matrix,detP = ±1). Speziell gilt noch fur (1.20), da P (PA) = A fur jedes AP 2 = 1, (⇒) P t = P .
1 Lineare Gleichungssysteme 12
Die Operation AP bewirkt ein Vertauschen der Spalte i mit der Spalte k vonA, da AP = (P tAt)t = (PAt)t = (a1, . . . , ak, . . . , ai, . . . , an). Daher liefertder GA mit ZV das Ergebnis
(1.21) (1− g(n−1))Pn−1 · · · (1− g(2))P2(1− g(1))P1A = R,
wobei Pi die i-te Zeile mit der Zeile pi ≥ i vertauscht. Wir notieren unsdaher die Matrix Pi uber die Zahl pi;
(1.22) P := Pn−1 · · ·P1 = (p1, . . . , pn−1).
Aus (1.21) erhalten wir aber die LR-Zerlegung von PA.Z.B. ist (n = 3)
(1− g(2))P2(1− g(1))P1 = (1− g(2))(P2 − P2g(1))P1
= (1− g(2))(P2 − P2g(1) P2P2︸ ︷︷ ︸
=1
)P1 = (1− g(2))(1− P2g(1)P2)P2P1
= (1− g(2))(1− P2g(1))P2P1,
da g(1)P2 = g(1), bzw. allgemein g(i)Pk = g(i) fur k > i gilt.Man zeigt daher durch Induktion die Beziehung
(1.23)(1− g(n−1))Pn−1 · · · (1− g(2))P2(1− g(1))P1 =
(1− g(n−1))(1− Pn−1g(n−2)) · · ·
(1− Pn−1 · · ·P3g(2))(1− Pn−1 · · ·P2g
(1))P1.
Da die Matrix Pn−1 · · ·Plg(l−1) von derselben Struktur wie g(l−1) ist, erhal-ten wir mit Lemma 1.1
Satz 1.2 Sei A ∈ Rn×n. Dann existieren Permutationsmatrizen Pi = pi, diedie Zeile i mit der Zeile pi ≥ i (i = 1, . . . , n − 1) vertauschen , so dass mitP := Pn−1 · · ·P1 gilt
PA = LR
mit der linken unteren (normierten) Dreiecksmatrix
(1.24) L = 1 + Pn−1 · · ·P2g(1) + Pn−1 · · ·P3g
(2) + . . .+ g(n−1)
und der rechten oberen Dreiecksmatrix R, die der GA mit der Zeilenvertau-schungsmatrix P liefert.
Mit δi = sgn (pi − i) ist detA = (−1)Pn−1i=1 δi
∏ni=1 rii, wobei rii die Diago-
nalelemente von R sind.
1 Lineare Gleichungssysteme 13
Wir werden gleich sehen, dass man aus Grunden der ”numerischen Stabi-litat“ i.A. den GA mit einer Zeilenvertauschungsstrategie durchfuhren muss.Es gibt jedoch eine große Klasse von Matrizen fur die der GA ohne ZVdurchfuhrbar und auch ”numerisch stabil “ ist. Dies sind die symmetrischpositiv definiten (spd) Matrizen.
Bekanntlich ist A ∈ Rn×n eine spd Matrix, falls At = A und (x,Ax) > 0fur x 6= 0. Eine symmetrische Matrix ist eine spd Matrix genau dann, wennalle Eigenwerte positiv sind; oder alle Hauptabschnittsdeterminanten positivsind.Ist A spd ⇒ aii > 0, denn mit x = ei ist (ei, Aei) = aii. Daher ist der 1.Schritt des GA (siehe (1.10)) durchfuhrbar, und es ist
G(1)A = (1− g(1))A =
a11 a12 · · · a1n
0... A(2)
0
= A(2)
An dem Formelsatz (1.14) sieht man, dass A(2)t = A(2) gilt. Ferner ist
(1.25) G(1)A(1)G(1)t =
a11 0 · · · 00... A(2)
0
= A(2)G(1)t,
denn A(2)G(1)t subtrahiert die Vielfachen αi aus (1.8) der ersten Spaltevon A(2) von der i-ten Spalte von A(2); nur die ersten Komponenten sindbetroffen, a1i = ai1.Sei x = (0, y)t, y 6= 0, y ∈ Rn−1.Dann ist (x,G(1)A(1)G(1)tx) = (y, A(2)y) = (G(1)tx,A(1)G(1)tx) > 0, daG(1)tx 6= 0.
Also ist A(2) spd, und wegen a(2)22 > 0 ist der GA ohne ZV fortsetzbar.
Offenbar folgt
G(2)A(2)G(1)t =
a
(1)11 0 · · · 00 a
(2)22 · · · a
(2)2n
... 0
...... A(3)
0 0
,
1 Lineare Gleichungssysteme 14
und daher
G(2)G(1)A(1)G(1)tG(2)t =
a
(1)11 0 · · · 00 a
(2)22 · · · 0
0 0...
... A(3)
0 0
,
usw. bis
(1.26) G(n−1) · · ·G(1)AG(1)t · · ·G(n−1)t = D = diag(r11, . . . , rnn)
Daher erhalten wir (vgl. Satz 1.1)
Satz 1.3 (Cholesky-Zerlegung) Mit den Bezeichnungen aus Satz 1.1 lau-tet die LR-Zerlegung einer spd Matrix A
(1.27) A = LDLt
mit D = diag(r11, . . . , rnn).
Beweis: Dies folgt sofort wie bei Satz 1.1 aus (1.26). 2
Das entscheidende ist, dass man die Cholesky-Zerlegung ungefahr mit derHalfte des Aufwands des GA herstellen kann:Es ist lkj = 0 fur j > k. Also folgt aus (1.27)
aik =mini,k∑j=1
lijrjjlkj .
Fur i ≥ k (nur diese brauchen wir) ist daher
aik =k∑j=1
lijrjjlkj ,
und es gilt
(1.28)
rkk = akk −k−1∑j=1
rjjl2kj
lik = (aik −k−1∑j=1
lijrjjlkj)1rkk
, i > k
1 Lineare Gleichungssysteme 15
Kennt man daher die Spalten 1, . . . , k − 1 von L und die Elemente r11, . . . ,rk−1,k−1, so ergibt sich rkk und die k-te Spalte von L aus (1.28)
(1.29) Cholesky-Verfahren
Fur k = 1, 2, . . . , n
rkk = akk −k−1∑j=1
rjjl2kj
Fur i = k + 1, . . . , n
lik = (aik −k−1∑j=1
lijrjjlkj)/rkk
Ende furEnde fur.
Bemerkung: Ist At = A und das Cholesky-Verfahren durchfuhrbar, so giltmit lii = 1 die Darstellung A = LDLt. A ist dann spd genau dann, fallsrkk > 0, k = 1, 2, . . . , n.(Denn es ist A
∣∣ik
= LDLt∣∣ik
fur i ≥ k, und, da beide Matrizen symmetrischsind, A = LDLt)
Bis jetzt haben wir immer in R gerechnet. Wir untersuchen jetzt, wie sichder GA auf dem Rechner verhalt. Denn dieser kann nur mit endlich vielenZahlen, den Maschinenzahlen , rechnen. Diese lauten in Gleitkommadarstel-lung bei einer Maschine mit m stelliger Mantisse (bei Verwendung der Basis10) ± 0.d1 . . . dm · 10e, d1 6= 0, di ∈ 0, . . . , 9 mit einem l-stelligen Exponen-ten e = ±
∑l−1k=0 ek · 10k, ek ∈ 0, . . . , 9.
Fast jeder reellen Zahl x, die betragsmaßig nicht zu klein und nicht zu großist, wird bekanntlich eine Maschinenzahl fl(x) zugeordnet, die von x denkleinsten Abstand hat, i.e. Runden. Es gilt
|fl(x)− x||x|
≤ eps = 5 · 10−m,
wobei eps die Maschinengenauigkeit ist. Es gilt stets
fl(x) = (1 + ε(x))x mit |ε(x)| ≤ eps.
Definition 1.1 Sei x, x ∈ R, x 6= 0,m ∈ N. x und x sind auf (m−1) Stellengenau, wenn |x−x||x| ≤ 10−m gilt.
1 Lineare Gleichungssysteme 16
(Sei x =∑∞
i=1 xi10−i · 10e und x = (1 + 10−m)x das schlechteste x, das dieUngleichung erfullt, so enthalt x eine Korrektur in der m-ten Stelle.)
Wir untersuchen daher zunachst das Verhalten der Losung von Ax = b(detA 6= 0 vorausgesetzt) bei kleinen Storungen der Daten (z.B. das Einlesenderselben in den Rechner). Es ist so, dass die Losung dramatisch verandertwerden kann!
(1.30) Beispiel:
A =(
1 11 1.0001
),detA 6= 0.
Die Losung von Ax = b mit b = (2, 2.0001)t lautet x = (1, 1)t. Nimmt manb = (2, 2.0002)t, so lautet die Losung x = (0, 2)t! Eine Anderung der Datenin der 5. Stelle ergibt ein Ergebnis, das mit x in keiner Stelle ubereinstimmt.
A ist schlecht konditioniert!
Woran kann man dies messen?Setzt man b = b + δb, so ist mit x = x + δx, A(x + δx) = b + δb bzw.Aδx = δb. Wegen Definition 1.1 interessiert der relative Fehler ‖δx‖
‖x‖ . Es giltδx = A−1δb, ‖δx‖ ≤ ‖A−1‖‖δb‖. Also ist
‖δx‖‖x‖
≤ ‖A−1‖‖δb‖‖x‖
= ‖A−1‖ ‖b‖‖x‖· ‖δb‖‖b‖
, ‖b‖ = ‖Ax‖ ≤ ‖A‖‖x‖
Damit erhalten wir
(1.31)‖δx‖‖x‖
≤ ‖A‖‖A−1‖‖δb‖‖b‖
.
Fur beliebige Matrizen A,B ist ‖AB‖ ≤ ‖A‖‖B‖ (Beweis?). Daher gilt‖I‖ = 1 = ‖AA−1‖ ≤ ‖A‖‖A−1‖.
Definition 1.2 Fur eine nicht singulare Matrix A ∈ Rn×n ist die Kondi-tion (bzgl. ‖ · ‖) von A, cond‖·‖(A) definiert durch
cond‖·‖(A) := ‖A‖‖A−1‖.
Die Kondition von A ist daher (die kleinste) obere Schranke fur den Ver-starkungsfaktor des relativen Eingabefehlers von b.
Fur Beispiel (1.30) ist cond∞(A) = 4 · 104, ‖δb‖‖b‖ ∼ 10−5, so dass sich mit
(1.31) ‖x−x‖‖x‖ ∼ 10−1 ergibt; keine Stelle ist genau.
1 Lineare Gleichungssysteme 17
Definition 1.3 Die Matrix A ∈ Rn×n ist schlecht konditioniert , wenn
cond‖·‖(A) 1.
Die Kondition cond∞(A) lasst sich verbessern durch Multiplikation der i-tenZeile mit di = (
∑k |aik|)−1, i.e. Zeilenaquilibrierung (Ubung).
Es ist zu erwarten, dass bei einer schlecht konditionierten Matrix jeder Algo-rithmus auf dem Rechner (d.h. endliche Stellenzahl) zum Losen von Ax = bein schlechtes Ergebnis liefert.Aber ein Grundproblem des Numerischen Rechnens (auf dem Rechner) ist,dass selbst bei gut konditionierten Problemen ”vernunftige“ (in R) Algo-rithmen vollig falsche Ergebnisse liefern konnen, i.e. numerische Instabilitat.Der GA ohne ZV ist ein Beispiel dafur.
(1.32) Beispiel:
A =(
10−5 11 1
), b =
(12
), cond∞A = 4.00004 . . .
A ist gut konditioniert. Bei exakter Rechnung (in R) erhalt man
x =(
1.000010 . . .0.9999899 . . .
).
Nimmt man 10−5 als Pivotelement , so folgt im ersten Gauß-Schritt dasSchema
10−5 1 1
105 1− 105 2− 105
Benutzt man jetzt zum Ausrechnen eine Maschine mit 4-stelliger Mantisse,so erhalt man aber
10−5 1 1
105 −105 −105,
woraus durch Ruckwartseinsetzen x =(
01
)folgt, was nichts mit dem
ursprunglichen Ergebnis zu tun hat.Der GA (ohne ZV) ist instabil. Durch den großen Multiplikator l21 lost mandas System
10−5 1 11 0 0
,
das eine zu große Storung des ursprunglichen ist.
1 Lineare Gleichungssysteme 18
Die folgende Strategie vermeidet große Multiplikatoren und hat eine stabi-lisierende Wirkung:
Spaltenpivotsuche: Im k-ten Schritt des GA nehme man ein betrags-großtes Element der ersten Spalte von A(k) als Pivotelement. D.h. bestimmep mit |a(k)
pk | = maxj=k,...,n |a(k)jk | und vertausche in A(k) bzw. A(k) die Zeile k
mit der Zeile p. (Es gilt dann |lik| < 1).
Auf unser Beispiel angewendet (k = 1) ergibt dies
1 1 2
10−5 1 1G(1)
−→1 1 2
10−5 1− 10−5 1− 2 · 10−5=
1 1 2
10−5 1 1
bei 4-stelliger Rechnung. Dies ergibt(
11
), was innerhalb der Stellenzahl
genau ist!
Definition 1.4 (Numerische Stabilitat eines Algorithmus)Ein Algorithmus ist stabil (auf dem Rechner), wenn der gesamte Rundungs-fehler durch die Kondition des mathematischen Problems bestimmt wird.
Es gelten die wichtigen Aussagen:
Satz 1.4
(i) Der GA mit Spaltenpivotsuche ist ein stabiler Algorithmus (mit patho-logischen Ausnahmen).
(ii) Die Spaltenpivotsuche ist nicht notig und nicht sinnvoll bei spd Ma-trizen und den diagonaldominanten Matrizen A, d.h.
∑nk=1,k 6=i |aik| ≤
|aii|, i = 1, . . . , n mit detA 6= 0.
Ist daher cond(A) 1, so fuhrt der GA mit Spaltenpivotsuche (wenn ange-bracht) i.a. zu großer Fehlerverstarkung, die im Fall cond(A) ∼ 1 moderatbleibt.Die Zeilenvertauschungen beschreiben wir mit den Matrizen Pi aus Satz 1.2.D.h. der GA mit Spaltenpivotsuche liefert die Zerlegung Pn−1 · · ·P1A = LRmit L = 1+
∑n−2k=1 Pn−1 · · ·Pk+1g
(k) +g(n−1)(n > 2); d.h. die k-te Spalte vonL wird bestimmt uber die k-te Spalte von g(k), k = 1, . . . , n − 1. Hat mang(k) berechnet, k = 1 . . . , n−2, d.h. den k-ten Eliminationsschritt vollzogen,so mussen alle spateren Vertauschungen Pk+1, . . . , Pn−1 auch in gk gebildetwerden, um die k-te Spalte von L zu erzeugen:
1 Lineare Gleichungssysteme 19
(1.33) Beispiel: Fortsetzung von Beispiel (1.18) (das Pivotelement ist je-weils fett gedruckt)
2 1 14 1 0−2 2 1
P1 b=p1=2−→4 1 02 1 1−2 2 1
G(1)
−→4 1 0
1/2 1/2 1−1/2 5/2 1
P2 b=p2=3−→4 1 0−1/2 5/2 11/2 1/2 1
G(2)
−→4 1 0−1/2 5/2 11/2 1/5 4/5
g(1) =
0 0 01/2 0 0−1/2 0 0
, g(2) =
0 0 00 0 00 1/5 0
damit gilt: P2P1A = LR mit
L = 1 + P2g(1) + g(2) =
1 0 0−1/2 1 01/2 1/5 1
, R =
4 1 00 5/2 10 0 4/5
.
Ein Gleichungssystem Ax = b lost man dann so:
(1) Sind die |aik| sehr unterschiedlich in der Großenordnung, so fuhrt maneine Zeilenaquilibrierung durch. Dieses System sei wieder mit Ax = bbezeichnet.
(2) Bestimme mit dem GA mit Spaltenpivotsuche (wenn notig) die Zerle-gung
(1.34) PA = LR, P = Pn−1 · · ·P1 = (p1, p2, . . . , pn−1).
Wegen PAx = LRx = Pb lose
(3) Ly = Pb (y ist eindeutig bestimmt) durch Vorwartseinsetzen (Pb istaus den pi einfach zu bestimmen), und lose Rx = y, wenn es geht,durch Ruckwartseinsetzen.
Ein Pseudoprogramm fur detA 6= 0 sieht etwa so aus:
1 Lineare Gleichungssysteme 20
Berechnung von L,R in A:
(1.35)
Fur k = 1, 2, . . . , n− 1 :bestimme p mit |apk| = maxi=k,...,n |aik|pk = pvertausche akj mit apj , j = 1, . . . , nwj = akj , j = k + 1, . . . , nFur i = k + 1, . . . , n :q = aik
akklik = qFur j = k + 1, . . . , n :aij = aij − qwj
Berechnung von Pb:Fur k = 1, . . . , n :
falls pk 6= kvertausche bk mit bpk
Vorwartseinsetzen:Fur i = 1, 2, . . . , n :yi = biFur j = 1, 2, . . . , i− 1 :yi = yi − aijyj
Ruckwartseinsetzen:
Fur i = n, n− 1, . . . , 1 :xi = yiFur j = i+ 1, . . . , n :xi = xi − aijxj
xi = xiaii
.
1.2 Iterative Verfahren
Fur große n z.B. n = 107 benotigt der GA ∼ 84 Tage, um ein vollbesetztesSystem zu losen. Selbst bei den dunnbesetzten Matrizen (engl. ”sparse ma-trix“), bei denen pro Zeile eine kleine Anzahl von Elementen nicht Null ist(A · x ∼ O(n)), sind direkte Verfahren wegen des ”fill in“, d.h. L und R ha-ben wesentlich mehr Nichtnullelemente, nicht geeignet. Die Sparse-Matrizentauchen bei der Diskretisierung von PDEs auf (siehe auch die Matrix A inder Einleitung).Man verwendet dann iterative Verfahren, die die Losung approximativ be-rechnen, d.h. bis auf eine gewisse Genauigkeit. Eine großere Genauigkeit als
1 Lineare Gleichungssysteme 21
die Modellgenauigkeit von Ax = b zu erzielen, macht naturlich keinen Sinn(Ubung).
Lineare Iterationsverfahren lassen sich so begrunden:Sei detA 6= 0 und x die Losung von Ax = b. Mit einer Ausgangsnaherungx0 erhalt man das Residuum r0 = b − Ax0. Ist r0 6= 0, so gilt mit x =x0 + δ, Aδ = b − Ax0, womit nichts gewonnen ist. Sei daher B eine ”gute“Approximation von A−1, wobei By genauso aufwendig wie Ay ist (A sparse⇒ B sparse). Meistens kennt man A und B nicht explizit sondern nur dieOperationen Ay,By.Man setzt daher δ = B(b−Ax0) und
(1.36)x1 = x0 +B(b−Ax0) bzw.
xj+1 = xj +B(b−Axj), j = 0, 1, . . ..
Es ist
(1.37) xj+1 = Txj +Bb
mit der Iterationsmatrix
(1.38) T = 1−BA (= 0, falls B = A−1).
Ferner gilt
(1.39) xj − x = T j(x0 − x), j = 0, 1, . . . .
Das Iterationsverfahren (1.36) konvergiert daher fur jedes x0 genau dann,falls ‖T jd‖ → 0, j →∞ fur jedes d ∈ Rn. Dies ist aber gleichbedeutend mit(Ubung)
(1.40) ‖T j‖ → 0, j →∞.
Bezeichnen λi(T ), i = 1, 2, . . . , n die Eigenwerte von T , so ist der Spek-tralradius ρ(T ) = maxi=1,...,n |λi(T )|. Wir normieren Cn durch ‖z‖Cn =(‖x‖2 + ‖y‖2)
12 , z = x+ iy. Dann gilt fur jede reelle Matrix A ∈ Rn×n
(1.41) ‖A‖Cn = ‖A‖
(Beweis: supz 6=0‖Az‖Cn‖z‖Cn
≥ supx 6=0‖Ax‖‖x‖ = ‖A‖, und ‖Az‖2Cn = ‖Ax‖2 +
‖Ay‖2 ≤ ‖A‖2(‖x‖2 + ‖y‖2) = ‖A‖2‖z‖2Cn)
Fur Matrizen T ∈ Cn×n gilt die wichtige Formel (Algebra?)
(1.42) ρ(T ) = limj→∞
‖T j‖1/jCn = infj‖T j‖1/jCn .
((1.42) ist sogar unabhangig von ‖ · ‖Cn)
1 Lineare Gleichungssysteme 22
Satz 1.5
(i) Das Iterationsverfahren (1.36) konvergiert (fur jedes x0 ∈ Rn) genaudann, wenn
(1.43) ρ(T ) < 1
gilt.
(ii) Wegen ρ(T ) ≤ ‖T‖ ist ‖T‖ < 1 hinreichend fur die Konvergenz.
Beweis:
(i) ⇒: Sei λ ∈ C ein EW von T zum EV y ∈ Cn mit ‖y‖Cn = 1. Es istwegen (1.41)
‖T jy‖Cn = |λ|j‖y‖Cn = |λ|j ≤ ‖T j‖ · 1→ 0
nach Voraussetzung. Also gilt (1.43).⇐: Mit (1.42) ist z.B. wegen (1.41) ‖T j‖ ≤ (1
2(ρ(T ) + 1))j fur j ≥ N0.
(ii) ist klar. 2
Ein Konvergenzkriterium (z.B. ‖T‖ < 1) kann man auch zum Nachweis vondetA 6= 0 benutzen:
Bemerkung: Ist ρ(T ) < 1 (wenn z.B. ein Konvergenzkriterium erfullt ist)und detB 6= 0, so gilt detA 6= 0.
Beweis: BA = 1− T . 0 ist kein EW von 1− T , denn sonst existiert x 6= 0mit Tx = x und somit ρ(T ) ≥ 1. Also existiert (BA)−1 2
Um die Bedingung ‖T‖ < 1 anzuwenden, benotigen wir
Lemma 1.2 Sei A ∈ Rn×n
(i) Fur ‖x‖1 =n∑i=1|xi| ist ‖A‖1 = max
1≤k≤n
n∑i=1|aik|
(ii) Fur ‖x‖2 =(
n∑i=1|xi|2
) 12
ist ‖A‖2 = λmax(AtA)12
(iii) Fur ‖x‖∞ = maxi=1,...,n
|xi| ist ‖A‖∞ = max1≤i≤n
n∑k=1
|aik|
(‖A‖ = supx 6=0
‖Ax‖‖x‖ )
1 Lineare Gleichungssysteme 23
Beweis: (i) und (ii) Ubung(iii) |
∑k
aikxk| ≤∑k
|aik||xk| ≤∑k
|aik|‖x‖∞ ≤ maxi
∑|aik|‖x‖∞ Also ist
‖Ax‖∞ ≤ maxi
∑k
|aik|‖x‖∞, und somit
(1) ‖A‖∞ ≤ maxi
∑k
|aik|.
Fur jedes x mit ‖x‖∞ = 1 ist ‖Ax‖∞ ≤ ‖A‖∞ · 1.Sei i0 ein Index mit max
i
∑k
|aik| =∑k
|ai0k|, und x definiert durch xk =
sgn ai0k, k = 1, 2, . . . , n. Dann ist ‖x‖∞ = 1 im Fall A 6= 0, und
‖Ax‖∞ ≥
∣∣∣∣∣∑k
ai0ksgn ai0k
∣∣∣∣∣ =∑k
|ai0k| = maxi
∑k
|aik|,
woraus mit (1) die Behauptung folgt. 2
Es folgen zwei grundlegende Iterationsverfahren:Das Gesamtschrittverfahren (GS)Ist aii 6= 0, so erhalt man ausgehend von der aquivalenten ”Fixpunktglei-chung“
(1.44) xi = − 1aii
n∑k=1k 6=i
aikxk − bi
, i = 1, 2, . . . , n
das Verfahren in der Form
(1.45)xj+1i = − 1
aii
n∑k=1k 6=i
aikxjk − bi
,
i = 1, 2, . . . , n; j = 0, 1, . . . ;x0 gegeben.
Der Aufwand pro Schritt ist O(n2) bei einer vollbesetzten Matrix und O(n)bei einer dunnbesetzten.Um das Verfahren in der Form (1.36) zu schreiben, setzten wir D = diag(A),
L =
0 · · · · · · 0
a21. . .
......
. . . . . ....
an1 · · · an,n−1 0
, R =
0 a12 · · · a1n...
. . . . . ....
.... . . an−1,n
0 · · · · · · 0
und es ist daher fur j = 0, 1, . . .
(1.46) xj+1 = −D−1(L+R)xj+D−1b = TGSxj+D−1b = xj+D−1(b−Axj),
1 Lineare Gleichungssysteme 24
wobei L + R = A − D benutzt wurde. Als Approximation von A−1 wirdD−1 genommen, was Sinn macht, wenn die Diagonalelemente uberwiegen,‖TGS‖∞ < 1 oder ‖TGS‖1 < 1. Dann ist das Verfahren konvergent .
Definition 1.5 A ∈ Rn×n mit aii 6= 0 erfullt das Zeilensummenkriteri-um , wenn
(1.47) ‖TGS‖∞ = max1≤i≤n
∑k=1k 6=i
∣∣∣∣aikaii∣∣∣∣ < 1
bzw. das Spaltensummenkriterium , wenn
(1.48) ‖TGS‖1 = max1≤k≤n
∑i=1i 6=k
∣∣∣∣ aikakk∣∣∣∣ < 1.
Bemerkung: Die Kriterien hangen zusammen in dem Sinn, dass (stets)‖TGS‖∞ = ‖T tGS‖1 gilt. T tGS ist aber nicht die Iterationsmatrix fur At. Manzeige (Ubung): Ist detB 6= 0, so gilt ρ(1−BtAt) = ρ(1−BA).
Das GS ist leicht parallel durchfuhrbar (”additives“ Verfahren) wohingegendas folgende ”bessere“ Einzelschrittverfahren (ES) nicht diese Eigenschaftbesitzt (”multiplikatives“ Verfahren):In (1.45) kann man auf die Idee kommen, die bereits ”verbesserten“ Kom-ponenten xj+1
l , l ≤ i− 1, zur Berechnung von xj+1i zu benutzen, d.h.
xj+11 = − 1
a11
(n∑k=2
a1kxjk − b1
)
xj+12 = − 1
a22
(a21x
j+11 +
n∑k=3
a2kxjk − b2
)usw.
(1.49)xj+1i = − 1
aii
(i−1∑k=1
aikxj+1k +
n∑k=i+1
aikxjk − bi
),
i = 1, . . . , n; j = 0, 1, . . . ;x0 gegeben
oder (D + L)xj+1 = −Rxj + b bzw.
(1.50)xj+1 = −(D + L)−1Rxj + (D + L)−1b
= TESxj + (D + L)−1b
= xj + (D + L)−1(b−Axj).
1 Lineare Gleichungssysteme 25
((D + L)−1 ∼ A−1)
Zur Approximation von A−1 werden mehr Informationen von A benutzt.Dies ergibt jedoch nicht notwendig eine bessere Konvergenz.
Wegen (1.39) ‖xj − x‖ ≤ ‖T j‖‖x0 − x‖ und (1.42) gilt ungefahr ‖xi − x‖ ∼ρ(T )j‖x0 − x‖, so dass die Fehlerreduktion von ρ(T ) bestimmt wird. DerAnfangsfehler ‖x0−x‖ wird auf ε‖x0−x‖ reduziert nach ∼ | log ε|
| log ρ| Iterationen.
Bemerkung: Das heute oft benutzte und gefeierte Multigridverfahren mitseiner uberlegenen Konvergenzeigenschaft ist – richtig aufgeschrieben – einES-Verfahren.Die parallelen Unterraumkorrektur-Verfahren sind GS-Verfahren.
Es gilt der nicht ganz einfach zu zeigende Konvergenzsatz
Satz 1.6 Sei A ∈ Rn×n
(i) A erfulle das Zeilensummen- oder Spaltensummenkriterium.Dann konvergiert das GS- (klar) und das ES-Verfahren.
(ii) Ist A eine spd Matrix, so konvergiert das ES-Verfahren.
2 Approximation im quadratischen Mittel 26
2 Approximation im quadratischen Mittel
Das Gleichungssystem Ax = b aus dem Beispiel der Computertomographiewar hoffnungslos uberbestimmt, d.h. A ∈ Rm×n mit m n. Ein solchesSystem ist selten losbar. Man sucht daher ein x, das wenigstens das Residu-um ‖Ax− b‖ minimal macht. Dasselbe Problem ergibt sich, wenn man z. B.ein physikalisches Gesetz, das die Form Φ(t) =
∑ni=1 ϕi(t)xi besitzt, durch
Messung fur ein bestimmtes Objekt bestimmen will; d.h. die Parameter xifestlegen mochte.Fur die Punkte t1 < t2 < . . . < tm macht man in tj αj Messungen f lj , l =1, 2, . . . , αj (die statistisch verteilt sind) und hatte gern
(2.1) Φ(tj) =n∑i=1
ϕi(tj)xi = f lj , l = 1, 2, . . . , αj , j = 1, 2, . . . ,m.
Dies ist ein uberbestimmtes System fur die xi, das i.a. nicht losbar ist. Hiermacht es Sinn, die Fehlerquadratsumme
m∑j=1
αj∑l=1
∣∣∣∣∣f lj −n∑i=1
ϕi(tj)xi
∣∣∣∣∣2
zu minimieren.In Matrixschreibweise bedeutet dies
(2.2) ‖Ax− b‖22 → min
mit offensichtlichen A und b.
y
tt1 t2 tj tm
←−Datenwolke
n∑i=1
ϕi(t)x∗if lj
Das Problem (2.2) werden wir jetzt genauer untersuchen.Mit A = (a1, . . . , an) und R(A) =< a1, . . . , an > ist g = Ax =
∑ni=1 a
ixi ∈< a1, . . . , an >. Damit ist (2.2) gleichbedeutend mit dem Approximations-problem
(2.3)Gesucht ist die beste Approximationg∗ = Ax∗ von b aus< a1, . . . , an >; ‖Ax∗ − b‖2 ≤ ‖Ax− b‖2, x ∈ Rn.
2 Approximation im quadratischen Mittel 27
Wir zeigen, dass (2.3) stets losbar ist (fur jede Norm).Sei RgA = l (≤ n) und o.E.< a1, . . . , an >=< a1, . . . , al >, A = (a1, . . . , al).Dann existiert
infg∈<a1,...,an>
‖b− g‖ = infx∈Rn
‖b−Ax‖ = infg∈<a1,...,al>
‖b− g‖ = infx∈Rl‖b− Ax‖.
Einen Minimalpunkt x∗ brauchen wir nur zu suchen in der Menge G =x | ‖Ax − b‖ ≤ ‖A0 − b‖ = ‖b‖, denn fur x /∈ G ist ‖Ax − b‖ > ‖b‖ =‖A0− b‖ ≥ inf x ‖Ax− b‖.Die Menge G ist abgeschlossen. G ist auch beschrankt, denn fur x ∈ G istwegen ‖Ax‖−‖b‖ ≤ ‖Ax−b‖ ≤ ‖b‖, ‖Ax‖ ≤ 2‖b‖. Und wegen ‖Ax‖ ≥ c‖x‖(‖A · ‖ ist eine Norm) ist G beschrankt. Also ist G kompakt.Da F (x) = ‖Ax− b‖2 stetig ist, existiert minx F (x) = F (x∗). Also ist (2.3)bzw. das Problem (2.2) stets losbar!
2.1 Normalgleichungen
Wie kann man eine Losung erhalten? Hier ist die Wahl der Norm entschei-dend. Fur ‖·‖ = ‖·‖2 ist F (x) := ‖Ax−b‖22 = (x,AtAx)−2(Atb, x)+‖b‖2 (furspatere Einsichten: F ist konvex und stetig differenzierbar). Notwendig fureinen lokalen Minimalpunkt (f.s.E.: der dann global ist) ist gradF (x∗) = 0(i.e. kritischer Punkt von F ), d.h.
(2.4) AtAx∗ = Atb (Normalgleichung (NG))
Jede Losung von (2.2) erfullt daher (2.4). (2.4) ist daher stets losbar (Ax = bunlosbar ⇒ AtAx = Atb losbar!).Die Normalgleichungen (2.4) sind mit g∗ = Ax∗ aquivalent zu (ai, Ax∗) =(ai, b), i = 1, . . . , n, ⇔ g∗ − b ⊥ < a1, . . . , an >, woher der Name kommt;
0
b
·Ax∗
< a1, . . . , an >= R(A)
(2.5) g∗ − b ⊥ < a1, . . . , an > (⇔ (2.4))
Ist detAtA 6= 0 (AtA ist eine quadratische Matrix), so ist die Losung von(2.2) die Losung von (2.4). Es gilt (?) detAtA 6= 0⇔ RgA = n ((AtA)ik =
2 Approximation im quadratischen Mittel 28
(ai, ak)).Ist detAtA = 0, so ist aber jede Losung von (2.4) eine Losung von (2.2) (dennf.s.E.: F konvex, gradF (x∗) = 0⇒ x∗ globaler Minimalpunkt): Wegen (2.5)ist g∗ − b⊥ g fur g ∈< a1, . . . , an >, und daher gilt (nach dem Pythagoras)‖(g∗ − b) + g‖22 = ‖g∗ − b‖22 + ‖g‖22 ≥ ‖g∗ − b‖22 oder ‖b − (g∗ + g)‖22 ≥‖b− g∗‖22, g ∈< a1, . . . , an >, woraus naturlich ‖b− g′‖22 ≥ ‖b− g∗‖22, g′ ∈<a1, . . . , an > folgt.
Damit haben wir das folgende Ergebnis bewiesen:
Satz 2.1 Sei A ∈ Rm×n,m ≥ n, b ∈ Rm.
(i) Das Problem (2.2), d.h. ‖Ax − b‖2 → min ist losbar. x∗ ist Losunggenau dann, wenn x∗ eine Losung der Normalgleichung (2.4)
AtAx = Atb
ist. Die Normalgleichungen sind eindeutig losbar genau dann, wenndetAtA 6= 0⇔ RgA = n gilt.
(ii) Das Element bester Approximation g∗ aus (2.3) ist eindeutig bestimmt(es gilt also g∗ = Ax∗1 = Ax∗2, x
∗1, x∗2 Losungen von (2.4)). g∗ heißt die
orthogonale Projektion von b auf < a1, . . . , an >. Ist RgA = n, sogilt die Darstellung
(2.6) g∗ = Pb, P = A(AtA)−1At,
und es ist
(2.7) ‖P‖2 = 1.
Beweis: Wir brauchen nur (ii) zu zeigen.Sei g1 = Ax∗1, g2 = Ax∗2 mit AtAx∗1 = Atb, AtAx∗2 = Atb,⇒ AtA(x∗1 − x∗2) =0⇒ (x∗1−x∗2, AtA(x∗1−x∗2)) = 0 = (A(x∗1−x∗2), A(x∗1−x∗2)) = ‖Ax∗1−Ax∗2‖2 =‖g∗1 − g∗2‖ ⇒ g∗1 = g∗2.Ist RgA = n, so lautet die Losung der Normalgleichungen x∗ =(AtA)−1Atb, und somit g∗ = Pb.
Es ist wegen (2.5) b−Pb⊥(−Pb). Also gilt (Pythagoras) ‖b−Pb+(−Pb)‖22 =‖b‖22 = ‖b − Pb‖22 + ‖Pb‖22 ≥ ‖Pb‖22, woraus ‖P‖2 ≤ 1 folgt. Anderer-seits ist P (Pb) = Pb fur jedes b, d.h. P 2 = P (Projektionsmatrix), woraus‖P 2‖2 = ‖P‖2 ≤ ‖P‖22 bzw. ‖P‖2 ≥ 1 folgt. 2
Bemerkung: Wichtig war die Wahl der Norm ‖ · ‖p mit p = 2 (Skalarpro-dukt-Struktur). Fur 1 < p < ∞, p 6= 2 erhalt man keine linearen notwendi-gen Bedingungen (2.4). Und die nicht differenzierbaren Falle p = 1,∞, die
2 Approximation im quadratischen Mittel 29
viele wichtige Anwendungen besitzen, fuhren jeweils auf ein lineares Opti-mierungsproblem mit Ungleichungsnebenbedingungen.
Wir untersuchen jetzt das Verhalten der Losung von AtAx = Atb bei kleinenStorungen der Daten b, A. Damit cond(AtA) <∞ ist, nehmen wir RgA = n;und die Kondition von AtA wird eine Rolle spielen.
Zunachst ist ‖Ax‖22 = (Ax,Ax) = (x,AtAx), so dass gilt
(2.8) λmin(AtA)1/2‖x‖2 ≤ ‖Ax‖2 ≤ λmax(AtA)1/2‖x‖2.
Ferner ist bekanntlich
(2.9) cond2(AtA) = λmax(AtA) / λmin(AtA).
Wir untersuchen nun Storungen von b (Storungen von A ergeben ein analo-ges Resultat, s. Literatur ).Sei AtAx = Atb, AtAx = Atb; also Ax = Pb,Ax = P b. Daher gilt ‖Ax −Ax‖ = ‖Pb− P b‖ und wegen (2.8)
‖x− x‖ ≤ λ−1/2min ‖Pb− P b‖,
‖x− x‖‖x‖
≤ λ−1/2min
‖Pb− P b‖‖Pb‖
‖Pb‖‖x‖
.
Wegen ‖Pb‖ = ‖Ax‖ ≤ λ1/2max‖x‖ erhalten wir
(2.10)‖x− x‖‖x‖
≤ cond(AtA)1/2 ‖Pb− P b‖‖Pb‖
.
Ist daher Pb = b, P b = b (⇔ b, b ∈ R(A)), so wird die Storungsempfind-lichkeit des mathematischen Problems bei Storungen von b ∈ R(A) (inR(A)) durch cond2(AtA)1/2 bestimmt (dies gilt auch bei festem b ∈ R(A)bei Storungen von A).
In (2.10) taucht etwas Neues auf: ‖Pb−P b‖‖Pb‖ kann aber im Verhaltnis zu ‖b−b‖‖b‖sehr groß werden (d.h. das math. Problem der Bildung Pb kann schlechtkonditioniert sein).
‖Pb− P b‖‖Pb‖
(2.7)
≤ ‖b− b‖‖b‖
‖b‖‖Pb‖
; (hier kann Gleichheit gelten),
und der Faktor ‖b‖‖Pb‖ kann sehr groß werden. Sei θ = ∠Pb, b(∈ [0, π/2]).
cos θ = (Pb,b)‖Pb‖‖b‖ = (Pb,Pb)
‖Pb‖‖b‖ = ‖Pb‖‖b‖ (da Pb− b⊥Pb).
Damit haben wir
Lemma 2.1 Unter den obigen Voraussetzungen (RgA = n, . . .) gilt
(2.11)‖x− x‖2‖x‖2
≤ 1cos θ
cond2(AtA)1/2 ‖b− b‖2‖b‖2
(Eine ahnliche Abschatzung gilt bei Storungen von A)
2 Approximation im quadratischen Mittel 30
Die Normalgleichungen sind daher schlecht konditioniert (cond2(AtA) ≥ 1),wenn θ ∼ π
2 gilt.
Beispiel:
A =
1 00 10 0
⇒ cond2(AtA) = 1,
b =
εε1
, b =
ε+ δε+ δ
1
, x =(εε
), x =
(ε+ δε+ δ
).
‖x− x‖2‖x‖2
=δ
ε,‖b− b‖2‖b‖2
=√
2δ(1 + 2ε2)1/2
,
cos θ =‖Pb‖2‖b‖2
=‖Ax‖2‖b‖2
=√
2ε(1 + 2ε2)1/2
,
‖Pb− P b‖2‖Pb‖2
=δ
ε=‖b− b‖2‖b‖2
1cos θ
.
Ist cos θ ∼ 1, dann richtet sich die Kondition nach (cond2(AtA))1/2, und wirsuchen jetzt ein numerisches Verfahren (numerisch stabil), dessen gesamterRundungsfehler durch (cond2(AtA))1/2 bestimmt wird.
Oft sind die Matrizen AtA, die das Leben schreibt, schlecht konditioniert:
Beispiel: (Weitere in der Ubung)Ein Fahrzeug wird aus der Geschwindigkeit v zum Stehen gebracht. Fur denBremsweg s nimmt man das folgende Gesetz an:
s = x1v2 + x2v + x3 = s(x1, x2, x3; v)
In einer Versuchsreihe wurden Werte (vi, si), i = 1, 2, . . . ,m = 5 ermittelt
vi (km/h) 9 17 17 25 35si 3 9 5 14 23
Das uberbestimmte System lautet
s(x1, x2, x3; vi) = si, i = 1, 2, . . . ,m
oder 81 9 1289 17 1289 17 1625 25 11225 35 1
x =
3951423
.
2 Approximation im quadratischen Mittel 31
AtA =
20648 69055 250969055 2509 1032509 103 5
, cond2(AtA) = 1.67 . . . · 107
x∗ =
0.0119 . . .0.2590 . . .−0.5247 . . .
, cos θ = 0.993 . . . .
Wendet man jetzt (also im Fall cos θ ∼ 1) zum numerischen Losen vonAtAx = Atb die Cholesky-Zerlegung an (?), so wird die Rundungsfehler-Fortpflanzung beginnend mit den Fehlern in AtA,Atb durch cond(AtA) be-stimmt (und nicht durch cond(AtA)1/2). Diese Vorgehensweise gilt daher (s.Def. 1.4) als numerisch instabil.(∗) Zudem ist die Skalarproduktbildung inAtA oft numerisch instabil: Sei ε = 10−3 und
A =
1 1 1ε 0 00 ε 00 0 ε
.
Es ist RgA = 3 und
AtA =
1 + ε2 1 11 1 + ε2 11 1 1 + ε2
=
1 1 11 1 11 1 1
(4-stellige Rechnung),
die Matrix AtA ist nach 4-stellige Rechnung nun singular. Um ein nume-risch stabiles Verfahren zu finden, gehen wir zu dem ursprunglichen Problemzuruck, namlich eine beste Approximation von b aus< a1, . . . , an > zu finden:
ming∈<a1,...,an>
‖b− g‖2.
Idee: Wahle eine ”bessere“ Basis fur < a1, . . . , an >.Sei q1, . . . , q
n, Q = (q1, . . . , qn) eine andere Basis fur < a1, . . . , an >. Danngilt bekanntlich
(2.12)ai =
n∑l=1
qlrli, i = 1, 2, . . . , n, mit det(rik) 6= 0, oder
A = QR, mit detR 6= 0 (R = (QtQ)−1QtA).
Es ist ming∈<a1,...,an> ‖b − g‖ = ming′∈<q1,...,qn> ‖b − g′‖, und das Problemwird gelost durch g∗ = Qy∗(= Ax∗) mit
(2.13) QtQy∗ = Qtb.
(∗)Man verwendet trotzdem oft die Cholesky-Zerlegung wegen des geringeren Aufwandsfur m n : mn2 + 1
3n3
2 Approximation im quadratischen Mittel 32
Wahlt man daher eine orthonormierte Basis fur R(A) =< a1, . . . , an >, d.h.(qi, qk) = δik oder Q ist eine orthogonale Matrix, QtQ = 1, so gilt y∗ = Qtbund Qy∗ = g∗ = Ax∗ = QRx∗, d.h.
(2.14) Rx∗ = Qtb.
Diese orthonormierte Basis existiert nach dem Gram-Schmidtschen Ortho-gonalisierungsverfahren (Algebra) mit einer rechten oberen DreiecksmatrixR, d.h. es gilt sogar
(2.15) ai =i∑l=1
qlrli mit rii > 0, i = 1, 2, . . . , n.
Und die Darstellung A = QR mit einer orthogonalen m × n Matrix Q undeiner rechten oberen Dreiecksmatrix R heißt QR-Zerlegung von A.
Entscheidend ist nun, dass cond2(R) mit der Kondition des mathematischenProblems (AtAx∗ = Atb) zusammenfallt:
Satz 2.2 Sei A ∈ Rm×n,m ≥ n mit RgA = n und A = QR eine QR-Zerlegung von A. Dann gilt
(2.16) AtAx∗ = Atb⇔ Rx∗ = Qtb
und
(2.17) cond2(R) = cond2(AtA)1/2.
D.h. kann man die QR-Zerlegung numerisch stabil herstellen (was der Fallist), dann ist das Ruckwartseinsetzen in (2.14) ein stabiler Algorithmus.
Beweis: Die Aquivalenz (2.16) zeigt man auch direkt durch Anwendungder QR-Zerlegung auf die Normalgleichungen, da AtA = RtQtQR = RtR,also AtAx∗ = Atb⇔ RtRx∗ = RtQtb⇔ Rx∗ = Qtb, da detR 6= 0.
‖R‖2 = λmax(RtR)1/2 = λmax(AtA)1/2 (RtR = AtA).‖R−1‖2 = λmax(R−tR−1)1/2 = λmax((RRt)−1)1/2.
Fur Eigenwerte gilt λi((RRt)−1) = 1λi(RRt)
, und wegen Rt(RRt)R−t = RtR
besitzen RRt und RtR dieselben Eigenwerte. Daher gilt λmax((RRt)−1)1/2 =1
λmin(RtR)1/2woraus (2.17) folgt. 2
Leider ist das Gram-Schmidtsche Orthogonalisierungsverfahren oft nume-risch instabil (die Auswertung von orthogonalen Projektionen kann, wie wirgesehen haben, instabil sein). Das Verfahren von Householder , mit dem wiruns jetzt befassen, ist numerisch stabil. Zunachst wollen wir die Bedingung
2 Approximation im quadratischen Mittel 33
RgA = n fallen lassen. Man sieht durch Betrachten des Gram-Schmidt Ver-fahrens, dass ein ONS q1, . . . , qn existiert mit < a1, . . . , ak >⊂< q1, . . . , qk >und ai =
∑il=1 rliq
l, i = 1, 2, . . . , k, rii ≥ 0, k = 1, 2, . . . , n. Also existiert all-gemein ein QR-Zerlegung von A:
(2.18) A = Q1R1, Q1 = (q1, . . . , qn), Qt1Q1 = 1, R1 =
r11 ∗. . .
0 rnn
,
die man auch zu einer QR-Zerlegung mit einer m×m Matrix Q = (Q1, Q2)erweitern kann
(2.19)A = (q1, . . . , qn, qn+1, . . . , qm) ·
r11 · · · r1n
. . ....rnn
0 · · · 0...
...0 · · · 0
= (Q1, Q2) ·R = QR
mit einer beliebigen orthogonalen Erganzung Q2.
Es gilt jetzt auch QQt = 1 und Qt ist auch eine orthogonale Matrix. Da-mit gilt, und das ist der Ausgangspunkt des Householder- (und Givens-)Verfahrens
(2.20) QtA = R.
Also gibt es eine orthogonale Matrix (Qt), die A auf obere Dreiecksformtransformiert. Sind B1, B2 ∈ Rn×n orthogonal so ist auch B1 ·B2 orthogonal.Deswegen versuchen wir jetzt, die Form (2.20) schrittweise zu erreichen.
2.2 Householder-Verfahren
Orthogonale Matrizen beschreiben typischerweise Drehungen und Spiege-lungen, denn (Qx,Qy) = (x,QtQy) = (x, y), so dass Winkel und Langenerhalten bleiben. Wir beschreiben das Householder-Verfahren durch Bilder:Sei A ∈ R5×4. Sei Q1 = H1 ∈ R5×5 eine Matrix, die die erste Spalte vonA1 = A in die Richtung von e1 ∈ R5 (stabil) dreht. Dann gilt
Q1A1 =
• · · ·0 · · ·0 · · ·0 · · ·0 · · ·
=
• · · ·00 A2
00
.
2 Approximation im quadratischen Mittel 34
Sei H2 ∈ R4×4 eine Matrix, die die erste Spalte von A2 in die Richtung vone1 ∈ R4 dreht. Dann gilt mit der orthogonalen Matrix Q2
Q2 =
1 0 0 0 000 H2
00
, Q2
(x1
y2
)=(
x1
H2y2
)
Q2Q1A =
• · · ·0 • · ·0 0 · ·0 0 · ·0 0 · ·
=
• · · ·0 • · ·0 00 0 A3
0 0
.Sei H3 ∈ R3×3 eine Matrix, die die erste Spalte von A3 in die Richtung vone1 ∈ R3 dreht. Dann gilt mit der orthogonalen Matrix Q3
Q3 =
1 0 0 0 00 1 0 0 00 00 0 H3
0 0
, Q3
x1
x2
y3
=
x1
x2
H3y3
Q3Q2Q1A =
• · · ·0 • · ·0 0 • ·0 0 0 ·0 0 0 ·
=
• · · ·0 • · ·0 0 • ·0 0 00 0 0 A4
.Sei H4 ∈ R2×2 eine Matrix, die die (erste) Spalte von A4 in die Richtungvon e1 ∈ R2 dreht. Dann gilt mit der orthogonalen Matrix Q4
Q4 =
1 0 0 0 00 1 0 0 00 0 1 0 00 0 00 0 0 H4
, Q4
x1
x2
x3
y4
=
x1
x2
x3
H4y4
Q4Q3Q2Q1A = R =
• · · ·0 • · ·0 0 • ·0 0 0 •0 0 0 0
.
2 Approximation im quadratischen Mittel 35
Allgemein lauten die orthogonalen Matrizen Qk ∈ Rm×m
(2.21) Qk =
1
. . . 01
0 Hk
mit einer (orthogonalen) Matrix
(2.22) Hk ∈ R(m−k)×(m−k),
die die erste Spalte von Ak in die Richtung von e1 ∈ Rm−k dreht.
Eine solche Matrix konstruieren wir jetzt (stabil).Das Householder-Verfahren benutzt Spiegelungen (Das Givens-Verfahrenbenutzt Drehungen). Sei E eine Hyperebene durch den Nullpunkt mit demNormalenvektor v, ‖v‖ = 1, d.h.
(2.23) E = x ∈ Rn | (x, v) = 0
Wir wollen a ∈ Rn(a 6= 0) an dieser Hyperebene spiegeln.
0Ha
← E
a
v
(a, v)v
·
·
Die orthogonale Projektion von a auf < v > ist (a, v)v. Daher ist Ha =a− 2(a, v)v der gespiegelte Vektor.
Ha = a− 2v(a, v) = a− 2vvta = (1− 2vvt)a.
Definition 2.1 Fur v ∈ Rn, ‖v‖2 = 1 ist die Householder Matrix H ∈Rn×n definiert durch
(2.24) H(v) = (1− 2vvt).
Man sieht sofort, Ht = H und H2 = 1, so dass H eine symmetrische,orthogonale Matrix ist.
2 Approximation im quadratischen Mittel 36
In unserem Fall mochten wir Ha = ±‖a‖2e1 haben.
−‖a‖e1 ‖a‖e1
a
v− v+
E+ = x | (x, v+) = 0
E− = x | (x, v−) = 0
0
Dies fuhrt auf die gesuchte Definition von v = w‖w‖2 mit
(2.25) w± = a± ‖a‖2e1.
Nun ist w± = (a1±‖a‖2, a2, . . . , an)t. Um hier ”Ausloschung“ (wir erlauterndies gleich) zu vermeiden – und dies macht das Verfahren numerisch stabil– wahlen wir
(2.26)w = w(a) = a+ sgn (a1)‖a‖2e1 (∗)
v = v(a) = w(a)‖w(a)‖2
H = H(v(a)).
Mit n = sgn (a1)e1 ist H = a − 2(a,w)‖w‖22
w und (a,w) = ‖a‖22 + ‖a‖2(a, n),
‖w‖22 = ‖a‖22 +2‖a‖2(a, n)+‖a‖22 = 2(a,w). Also gilt H(v(a))a = a−w(a) =−‖a‖2n, oder
(2.27) H(v(a))a =−‖a‖2e1, a1 > 0, Spiegelung an E+
‖a‖2e1, a1 ≤ 0, Spiegelung an E−.(∗∗)
(2.28) H(v(a))x = x− (x,w)‖a‖2(‖a‖2 + |a1|)
w
Damit ist klar, wie mit den Matrizen (2.21) und H aus (2.26) das Househol-derverfahren durchfuhrbar ist. Fur eine effektive Durchfuhrung s. Literatur,bzw. [8].
(∗)sgn (a1) =
1, a1 > 0−1, a1 ≤ 0
(∗∗)Bei Gram-Schmidt ist stets rii > 0
2 Approximation im quadratischen Mittel 37
Wir fassen zusammen:
Satz 2.3 (i) Sei A ∈ Rm×n(m ≥ n). Dann existiert eine QR-Zerlegungvon A, die man sich numerisch stabil mit dem Householder-Verfahren(oder dem Givens-Verfahren) erzeugen kann.
(2.29)
A = QR, Q ∈ Rm×m, QtQ = 1, R ∈ Rm×n, R =(R1
0
),
R1 ∈ Rn×n, R1 =
r11 · · · r1n
. . ....
0 rnn
.
Mit Q = (q1, . . . , qn, qn+1, . . . , qm) = (Q1, Q2) gilt
(2.30) A = Q1R1 (kleine QR-Zerlegung)
Der Aufwand fur das Householder-Verfahren ist
m n : 2n2mm ' n : 21
3n3 (2 mal LR-Aufwand, m = n).
(ii) Es ist RgA = RgR = RgR1.Ist RgA = n, so gilt
(2.31) < a1, . . . , ak >=< q1, . . . , qk >, k = 1, 2, . . . , n
(2.32) R(A)⊥ = x | (x, y) = 0, y ∈ R(A) =< qn+1, . . . , qm > .
Q1 liefert daher eine Orthonormalbasis fur R(A) und Q2 eine furR(A)⊥ (Rm = R(A)⊕R(A)⊥).
Beweis: Nur Teil (ii). Wegen A = Q1R1 und RgA = n ist detR1 =∏ni=1 rii 6= 0, und daher < q1, . . . , qk >⊂< a1, . . . , ak >, k = 1, 2, . . . , n.
Daher gilt (2.31).Bekanntlich ist Rm = R(A) ⊕ R(A)⊥, und daher dimR(A)⊥ = m − n.qn+1, . . . , qm ∈ R(A)⊥ und sie sind linear unabhangig, also gilt (2.32). 2
Zur Ausloschung : Dies ist eine haufige Ursache fur die numerische Instabi-litat eines Algorithmus. Die Subtraktion fast gleich großer Maschinenzahlenauf dem Rechner kann zu erheblichem Genauigkeitsverlust fuhren. DieserFehler vererbt sich naturlich auf die weitere (noch so genaue) Rechnung.
Beispiel: Fur f ∈ C3(R) ist
f(t+ h)− f(t− h)2h
= f ′(t) +h2
6f (3)(η), η ∈ (t− h, t+ h), h > 0
2 Approximation im quadratischen Mittel 38
Sei f(x) = ex, t = 1, h = 10−5. Dann ist∣∣∣∣f(1 + h)− f(1− h)2h
− f ′(1)∣∣∣∣ 1|f ′(1)|
< 0.166 . . . · 10−10 ∼ 1.6 · 10−11,
und daher ist f(1+h)−f(1−h)2h bei Rechnung in R auf ∼ 10 Stellen genau.
rechnet man mit 10-stelliger Mantisse, so folgt
0.2718309011 · 10−1 − 0.2718254646 · 10−1
2 · 10−5= 0.27183 · 10−1.
Wegen e = 2.718281828 . . . sind aber nur ∼ 5 Stellen richtig!Hier ist x = fl(x) = (1 + εx)x, |εx| ≤ 5 · 10−10, x ⊕ y = (1 + ε+)(x + y),mit der Maschinenoperation ⊕. Fur δ definiert durch x⊕ y = (1 + δ)(x+ y)findet man daher δ =
(εxx+yx+ εy
x+yy)
(1 + ε+) + ε+.
Wegen x+ y ∼ 10−5 ist xx+y ,
yx+y ∼ 10−5 und also δ ∼ 10−5.
Zur Prazisierung von Satz (2.2) (Losen der Normalgleichung) haben wir nun
Satz 2.4 Sei A ∈ Rm×n(m ≥ n) mit RgA = n, und A = QR die QR-Zerlegung von A nach Householder (oder Givens) mit Q = (Q1, Q2), Q1 ∈
Rm×n, Q2 ∈ Rm×(m−n) und R =(R1
0
)mit R1 ∈ Rn×n. Dann lost man die
Normalgleichung AtAx∗ = Atb mit der QR-Zerlegung A = Q1R1, d.h.
(2.33) R1x∗ = Qt1b.
Ferner gilt
(2.34) ‖Ax∗ − b‖2 = ‖Pb− b‖2 = ‖Q2b‖2
und mit θ = ∠Pb, b, sin2 θ = ‖Q2b‖22‖b‖22
(s. Lemma 2.1).
Beweis: Nur von (2.34).
‖Ax− b‖22 = ‖QRx−QQtb‖22 = ‖Rx−Qtb‖22. Rx =(R1x
0
), Qt =
(Qt1bQt2b
).
Also ist ‖Ax − b‖22 = ‖R1x − Qt1b‖22 + ‖Qt2b‖22, woraus ((2.33)) (2.34) folgt.‖b‖22 = ‖Pb− b− Pb‖22 = ‖Pb− b‖22 + ‖Pb‖22; cos θ = ‖Pb‖2
‖b‖2 2
Bis jetzt haben wir bei der Losung der Normalgleichungen RgA = n ange-nommen. Zum Abschluss des Kapitels machen wir noch einige Ausfuhrungen
2 Approximation im quadratischen Mittel 39
zum Fall RgA = n− l, l ≥ 1, d.h. detAtA = 0.
Die Normalgleichungen waren ja stets losbar (was man auch so sehen kann:b = b1 + b2, b1 ∈ R(A), b2 ∈ R(A)⊥ = N(At) (?) ⇒ Atb = Atb1, b1 = Ax⇒x ist Losung der Normalgleichung). Sie sind aber nicht eindeutig losbar.Es gilt
(2.35) N(AtA) = N(A).
Denn: x ∈ N(A), d.h. Ax = 0,⇒ x ∈ N(AtA). Ist x ∈ N(AtA) ⇒ AtAx =0 ⇒ (x,AtAx) = (Ax,Ax) = ‖Ax‖22 = 0 ⇒ x ∈ N(A). Sei x0 eine festeLosung der NG, AtAx0 = Atb, und x∗ eine beliebige. Dann gilt AtA(x0 −x∗) = 0, d.h. (s. (2.35)) x∗ − x0 ∈ N(A) =< u1, . . . , ul > . Also gilt dieDarstellung
(2.36) x∗ = x0 −l∑
i=1
αiui.
Wir wollen b die Losung mit der kleinsten ‖ · ‖2-Norm zuordnen, d.h.
‖x0 −l∑
i=1
αiui‖2 → min .
Dieses Problem ist aber nach Satz 2.1 eindeutig losbar : x∗min = x0−∑l
i=1 α∗i ui,
wobei x∗min charakterisiert ist (Normalgleichungen) durch (2.5),
(2.37) x∗min⊥N(A).
Wir haben daher die Aussage:Die Losung x∗min der NG mit der kleinsten ‖ · ‖2-Norm ist die Losung von
(2.38)AtAx = Atb
x⊥N(A), d.h. x ∈ N(A)⊥.
Wir konnen daher jedem b ∈ Rm die Losung x von (2.38) zuordnen(∗):
(2.39) x∗min = A+b.
Man zeigt leicht (Ubung), dass die Abbildung A+ linear ist, also ist A+ eineMatrix.
(∗)Die Eindeutigkeit der Losung sieht man auch so: Sind x1, x2 Losungen von (2.38)⇒ x1 − x2 ∈ N(A) und x1 − x2 ∈ N(A)⊥ ⇒ x1 − x2 = 0.
2 Approximation im quadratischen Mittel 40
Definition 2.2 Die Matrix A+ ∈ Rn×m in (2.39) heißt die Pseudoinversevon A.
Man sieht sofort die Eigenschaften
(2.40) N(A) = 0 (⇔ RgA = n)⇒ A+ = (AtA)−1At,
und ist zusatzlich A ∈ Rn×n, so gilt A+ = A−1.
Die numerische Berechnung der Pseudoinversen erfolgt naturlich nicht uberdas Losen von (2.38), sondern uber eine weitere wichtige Zerlegung der Ma-trix A, die Singularwertzerlegung SVD (Singular Value Decomposition). Wirgeben hier nur die wichtigsten Tatsachen an und verweisen fur die numeri-sche Realisierung wie immer auf [8].
Definition 2.3 Sei A ∈ Rm×n. Die Singularwerte von A, σi(A), sind defi-niert durch σi(A) = λi(AtA)1/2, i = 1, 2, . . . , n.
Setzt man mit den Eigenvektoren vi (zu den Eigenwerten λi(AtA)) V =(v1, . . . , vn), so ist (bekanntlich) V tAtAV = diag(σ2
1, . . . , σ2n), und daher
R(AtAV ) = R(AtA) =< σ21v
1, . . . , σ2nv
n >, woraus
(2.41) RgA(= Rg (AtA)) = #i |σi > 0
folgt.
Sei RgA = r und o.E. σ1 ≥ σ2 ≥ . . . ≥ σr > σr+1 = 0 = . . . = σn. Wir
bestimmen zunachst die Pseudoinverse der Diagonalmatrix D =(D1
0
)∈
Rm×n(m ≥ n) mit D1 = diag(σ1, . . . , σn).
(2.42) D =
m
σ1
. . . 0σr
0
0 . . .0
0
n
n
ergibt
2 Approximation im quadratischen Mittel 41
D+ = (D+1 , 0) =
n m
n
1σ1
. . .1σr
0 0. . .
0
Denn ‖Dx − b‖22 =
∑ri=1(σixi − bi)2 +
∑mi=r+1 b
2i ist minimal fur jedes x∗
der Form x∗ = ( 1σ1b1, . . . ,
1σrbr, xr+1, . . . , xn)t. Die Losung mit der kleinsten
‖ · ‖2-Norm ist also
(2.43) x∗min = (1σ1b1, . . . ,
1σrbr, 0, . . . , 0)t ≡ D+b,
d.h. es gilt (2.42).
Satz 2.5 (Singularwertzerlegung, SVD) Sei A ∈ Rm×n(m ≥ n) mitRgA = r. Dann existieren orthogonale Matrizen U = (u1, . . . , um) ∈ Rm×m,
V = (v1, . . . , vn) ∈ Rn×n und eine Diagonalmatrix D =(D1
0
)∈ Rm×n mit
D1 = diag(σ1, . . . , σn), wobei σ1 ≥ σ2 ≥ . . . ≥ σr > σr+1 = 0 = . . . = σn, sodass gilt
(2.44) A = UDV t (⇒ AAtU = UDDt).
Mit U1 = (u1, . . . , un) gilt die ”kleine SVD“
(2.45) A = U1D1Vt (AtAV = V D2
1).
Aus Satz (2.5) ergibt sich sofort (nur f.s.E.)
Lemma 2.2 Mit den Bezeichnungen von Satz (2.5) gilt
(i) Ax =r∑i=1
σi(x, vi)ui, x ∈ Rn, Aty =r∑i=1
σi(y, ui)vi, y ∈ Rm
(ii) R(A) =< u1, . . . , ur >, R(At) =< v1, . . . , vr >,N(A) =< vr+1, . . . , vn >, N(At) =< ur+1, . . . , un > .
Mit Hilfe von (2.44) (⇔ (i)) erhalten wir jetzt eine Darstellung der Pseudoin-versen A+: (‖·‖ = ‖·‖2) ‖Ax−b‖2 = ‖UDV tx−b‖2 = ‖U(DV tx−U tb)‖2 =
2 Approximation im quadratischen Mittel 42
‖DV tx − U tb‖2 = ‖Dy − U tb‖2 mit y = V tx. Dy =(D1y
0
), also ist
‖Ax − b‖2 =∑r
i=1(σiyi − (ui, b))2 +∑m
i=r+1(ui, b)2. Dieser Ausdruck wirdminimal fur jedes y∗ der Form
y∗ =(
1σ1
(u1, b), . . . ,1σr
(ur, b), yr+1, . . . , yn
)t,
bzw. x∗ = V y∗.
Wegen ‖x∗‖ = ‖y∗‖ ist die Losung mit der kleinsten Norm
y∗min =(
1σ1
(u1, b), . . . ,1σr
(ur, b), 0, . . . , 0)t
= D+U tb = D+1 U
t1b,
x∗min = V D+U tb = V D+1 U
t1b.
Satz 2.6 Sei A ∈ Rm×n(m ≥ n) mit RgA = r und der SVD A = UDV t =U1D1V
t. Dann gilt
(2.46) A+ = V D+U t = V D+1 U
t1
mit der Pseudoinversen D+, D+1 nach (2.42).
Bei Rangabfall (was nicht einfach festzustellen ist) verwendet man daherzur Losung der Normalgleichung die SVD Zerlegung A = U1D1V
t. DerAlgorithmus ist fur m n so teuer wie die QR-Zerlegung (s. [8]).
(2.46) erhalt man auch direkt durch Anwendung der SVD auf die Normal-gleichungen: A = U1D1V
t, c := U t1b. AtAx∗ = Atb ⇔ V D2
1Vtx∗ = V D1c ⇔
D21V
tx∗ = D1c. Setzte y∗ = V tx∗ und lose das Diagonalsystem
(2.47) D21y∗ = D1c⇔ y∗i =
1σici, i = 1, 2, . . . , r.
y∗min = ( 1σ1c1, . . . ,
1σrcr, 0, . . . , 0), x∗min = V y∗min =
∑ri=1
1σi
(ui, b)vi.
3 Nichtlineare Gleichungssysteme 43
3 Nichtlineare Gleichungssysteme
Wir beginnen mit einem Beispiel. Gesucht sind Losungen des Systems
(3.1)x1 +
140x2
1 +150x2
2 −12
=: f1(x1, x2) = 0
x2 +140x2
1 +120x2
2 −12
=: f2(x1, x2) = 0.
Zunachst ist uberhaupt nicht klar, ob Losungen existieren; und es gibt imGegensatz zu linearen Gleichungen kein allgemeines Konzept, dies zu ent-scheiden. Hier sieht man, dass man zwei Ellipsen zum Schnitt bringen muss:Denn die erste Gleichung in (3.1) lautet
(x1 + 20)2
420+
x22
525= 1, und die zweite
x21
220+
(x2 + 10)2
110= 1.
Man sieht, dass das Nullstellenproblem(3.1) genau zwei Losungen besitzt.(Durch Verschieben der Ellipsen kann man sich leicht ein Problem konstru-ieren, dass eindeutig losbar, bzw. nicht losbar ist.)
Wie kann man (3.1) numerisch losen? Es bietet sich ein iteratives Vorge-hen an. Lost man (s. Herleitung des GS bei linearen Gleichungen) die ersteGleichung nach x1 und die zweite nach x2, so erhalt man eine aquivalenteFixpunktgleichung (davon gibt es beliebig viele)
(3.2)x1 = − 1
40x2
1 −150x2
2 +12
=: ϕ1(x1, x2)
x2 = − 140x2
1 −120x2
2 +12
=: ϕ2(x1, x2),
oder
(3.3) x = ϕ(x), ϕ = (ϕ1, ϕ2)t.
Jede Losung von (3.3) heißt Fixpunkt .
Wie bei dem Gesamtschrittverfahren kann man jetzt die Fixpunktiteration
(3.4) xk+1 = ϕ(xk), k = 0, 1, . . .
versuchen. Der Ausgangsnaherung x0 kommt hier aber eine besondere Be-deutung zu; auch bei nur einer Losung. Es ist daher keine ”globale Konver-genz“ zu erwarten.
3 Nichtlineare Gleichungssysteme 44
Fur (3.1) sieht man, dass eine Nullstelle von f = (f1, f2)t bzw. ein Fixpunktin der abgeschlossenen (und beschrankten) Menge
(3.5) G = [0, 0.5]× [0, 1]
liegt. Wegen 0 < ϕ1(x1, x2) < 12 , 0 < ϕ2(x1, x2) < 1
2 , (x1, x2) ∈ G gilt
(3.6) ϕ(G) ⊂ G.
(F.s.E.: Da G kompakt und konvex, ϕ stetig, folgt aus dem Fixpunktsatzvon F.E. Browder, dass wegen (3.6) ein Fixpunkt existiert.) Also liegt dieFolge (3.4) in G, falls x0 ∈ G.
Mit x0 = (0.4, 0.4)t ”steht“ die Iteration (bei 10 stelliger Mantisse) beix9 = 0.4893593998, 0.48237877228)t. x9 ist aber eine Approximation desFixpunktes in G. Denn aus limk→∞ x
k = x folgt x ∈ G (da G abge-schlossen), und wegen (3.4) und der Stetigkeit von ϕ limk→∞ x
k+1 = x =limk→∞ ϕ(xk) = ϕ(x).
Wie kann man nun garantieren (durch eine Bedingung an ϕ), dass die Folgeaus (3.4) konvergiert? Dies wurde dann auch die Existenz eines Fixpunktesliefern.
Nach Stefan Banach lautet eine Bedingung: (Bild fur n=1?)ϕ ist kontrahierend auf G, d.h. es ex. 0 < κ < 1 mit
(3.7) ‖ϕ(x)− ϕ(y)‖ ≤ κ‖x− y‖, x, y ∈ G.
(Beachte: κ = κ(‖ · ‖))Ist diese Bedingung erfullt, so kann ϕ nur einen Fixpunkt in G besitzen.Sind namlich x, y zwei Fixpunkte mit ‖x− y‖ > 0, so folgt wegen (3.7) derWiderspruch ‖x− y‖ < ‖x− y‖.
Bemerkung: Fur die Iterationsverfahren (1.36) ist die Bedingung (3.7)gleichbedeutend mit ‖T‖ < 1.
Fur das Beispiel (3.2) ist die Kontraktionsbedingung fur z.B. ‖ · ‖ = ‖ · ‖∞erfullt:
Sei x, y ∈ G
|ϕ1(x)− ϕ1(y)| =∣∣∣∣(− 1
40x2
1 −150x2
2 +12
)−(− 1
40y2
1 −150y2
2 +12
)∣∣∣∣≤ 1
40|x2
1 − y21|+
150|x2
2 − y22|
=140|x1 − y1||x1 + y1|+
150|x2 − y2||x2 + y2|
3 Nichtlineare Gleichungssysteme 45
≤ ‖x− y‖∞(
140
+150
)=
9200‖x− y‖∞.
|ϕ2(x)− ϕ2(y)| ≤ 340‖x− y‖∞;
also gilt
‖ϕ(x)− ϕ(y)‖∞ ≤ κ‖x− y‖∞ , κ =340.
Die Bedingung (3.7) garantiert nun, dass die Folge (3.4) (unter der Voraus-setzung (3.6)) das Cauchy-Kriterium erfullt:Man benutzt den Trick uber die Teleskopreihe. Sei m > n. Dann gilt
‖xm − xn‖ = ‖m−1∑k=n
(xk+1 − xk)‖ ≤m−1∑k=n
‖xk+1 − xk‖
und wegen (3.7) geht ‖xk+1 − xk‖ hinreichend schnell gegen Null;
(3.8) ‖xk+1−xk‖ = ‖ϕ(xk)−ϕ(xk−1)‖ ≤ κ‖xk−xk−1‖ ≤ . . . ≤ κk‖x1−x0‖.
Damit haben wir aber
(3.9) ‖xm − xn‖ ≤ ‖x1 − x0‖m−1∑k=n
κk,
und wegen 0 < κ < 1 ist∑∞
k=0 κk = 1
1−κ , so dass das Cauchy-Kriteriumerfullt ist. Also existiert limk→∞ x
k = x ∈ G und ist dort der einzige Fix-punkt.
Damit haben wir den beruhmten Fixpunktsatz von S. Banach bewiesen, derauch von großer theoretische Bedeutung ist.Wir fassen zusammen:
Satz 3.1 (Banachscher Fixpunktsatz) Sei G ⊂ Rn eine abgeschlosseneMenge, und die Funktion ϕ : G→ Rn erfulle die folgenden Bedingungen:
(i) ϕ(G) ⊂ G
(ii) ‖ϕ(x)− ϕ(y)‖ ≤ κ‖x− y‖, x, y ∈ G mit 0 < κ < 1
Dann gilt: Die Folge xk∞k=0, xk+1 = ϕ(xk), k = 1, 2, . . . , x0 ∈ G konver-
giert gegen den einzigen Fixpunkt x von ϕ in G.
3 Nichtlineare Gleichungssysteme 46
Ferner gelten die a posteriori bzw. a priori Fehlerabschatzungen
(3.10) ‖xk − x‖ ≤ κ
1− κ‖xk − xk−1‖ ≤ κk
1− κ‖x1 − x0‖, k = 0, 1, . . . .
Insbesondere ist daher ‖xk−xk−1‖ ein Maß fur die Genauigkeit von xk, unddie Folge xk konvergiert (mindestens) linear gegen x, da
(3.11) ‖xk+1 − x‖ ≤ κ‖xk − x‖.
Beweis: Wir brauchen nur (3.10) und (3.11) zu zeigen.
‖xk − x‖ = ‖ϕ(xk−1)− ϕ(x)‖ ≤ κ‖xk−1 − x‖ ≤ κ‖xk − xk−1‖+ κ‖xk − x‖,
woraus (1− k)‖xk− x‖ ≤ κ‖xk−xk−1‖ folgt. Setzt man hier noch (3.8) ein,so ergibt sich die a priori Abschatzung. 2
Die lineare Konvergenz in (3.11) kann sehr langsam sein;z.B. x = cosx, xk+1 = cosxk, k = 0, 1, . . . mit x0 = 0.6 konvergiert sehrlangsam gegen den Fixpunkt. Ist namlich κ = 10−α, α > 0 und ‖xk−x‖
‖x‖ ∼10−(m+1), d.h. xk ist auf m Stellen genau, so gewinnt man fur α ≥ 1, α =[α]+β (0 ≤ β < 1) pro Schritt [α] Stellen. Ist α < 1 ( 1
α > 1), so gewinnt man
jeweils nach jα =
1/α, 1/α ∈ N[1/α] + 1, sonst
(wegen jα ·α ≥ 1) Iterationen ∼ 1
Stelle hinzu. ([x] = großte ganze Zahl ≤ x).
Ein ”besser“ gewahltes ϕ als in (3.2) kann fur schnellere Konvergenz sorgen.
Zunachst
Definition 3.1 Eine iterativ definierte Folge xk∞k=0 mit xk+1 = ϕ(xk),k = 0, 1, . . . , limk→∞ x
k = x = ϕ(x), xk 6= x, k ∈ N0 konvergiert mindestensmit der Ordnung p ≥ 1, wenn
limn→∞
‖xk+1 − x‖‖xk − x‖p
= c ≥ 0
und im Fall p = 1, 0 ≤ c < 1 gilt.Im Fall c > 0 (dann ist p eindeutig) heißt p die Konvergenzordnung.
Ein nichtstationares Iterationsverfahren ist daher konvergent mit der Ord-nung p ≥ 1, falls fur k ≥ k0
(3.12) ‖xk+1 − x‖ ≤ c‖xk − x‖p
mit c < 1 im Fall p = 1 gilt.
3 Nichtlineare Gleichungssysteme 47
Ein Verfahren, das oft quadratisch (p = 2) konvergiert, ist das Newton-Verfahren zur Losung von f(x) = 0, d.h.
(3.13)f1(x1, . . . , xn) = 0
...fn(x1, . . . , xn) = 0.
Wir setzen also voraus, dass f eine Nullstelle x besitzt. Sei x0 eine guteAusgangsnaherung von x. Dann kann man f(x) in einer Umgebung von x0
(in der sich x befindet) gut approximieren (Modell) durch ihre Linearisierungin x0:
(3.14) f(x) = f(x0) + f ′(x0)(x− x0) +O(‖x− x0‖), f ∈ C1.
Als eine ”bessere“ Approximation von x nehmen wir die Nullstelle der Li-nearisierung, x1. Also
(3.15) 0 = f(x0) + f ′(x0)(x1 − x0).
usw.
Das Verfahren lautet daher
(3.16)
Fur k = 0, 1, . . .
f ′(xk)∆xk = −f(xk)
xk+1 = xk + ∆xk.
Man hat daher pro Schritt ein lineares Gleichungssystem zu losen. Haltman f ′ konstant fur einige Iterationen, so spricht man vom vereinfachtenNewton-Verfahren. Fur das Beispiel (3.1) erhalt man mit x0 = (0.4, 0.4)t
bei 10 stelliger Mantisse
x1 =
0.489......
0.482...
, x2 =
0.489359......
0.482378...
, x3 =
0.4893593999...
0.4823787230
.
Man sieht die schnelle quadratische Konvergenz in der Verdoppelung dergultigen Ziffern pro Schritt. Denn ist p = 2 in (3.12) und ‖x
k−x‖x‖ ∼ 10−(m+1),
so folgt ‖xk+1−x‖‖x‖ ∼ 10−1 · 10−(2m+1) ∼ 10−(2m+2).
Wir haben den folgenden lokalen (weil x0 hinreichend genau sein muss)Konvergenzsatz.
Satz 3.2 Es sei D ⊂ Rn offen, f : D → Rn sei stetig differenzierbar (d.h.∂fi∂xk
(x) sei stetig in D, d.h. f ∈ C1(D)). f besitze in D eine einfache
3 Nichtlineare Gleichungssysteme 48
Nullstelle x, d.h. f ′(x)−1 existiert. f ′ erfulle auf einer Kugel BR(x) ⊂ D dieLipschitzbedingung
(3.17) ‖f ′−1(z)(f ′(z)− f ′(x))‖ ≤ l‖z − x‖, z, x ∈ BR(x).(∗)
Fur die kleinere Umgebung Br(x) ⊂ BR(x) mit r < R und
(3.18) γ := 2rl < 1
hat man die Aussagen:Fur x0 ∈ Br(x) liegt die Folge xk∞k=0 definiert durch
(3.19)f ′(zk)∆xk = −f(xk), zk ∈ Br(x)
xk+1 = xk + ∆xk, k = 0, 1, . . .
in Br(x) und konvergiert mindestens linear gegen x,
(3.20) ‖xk+1 − x‖ ≤ γ‖xk − x‖, k = 0, 1, . . . .
Fur zk = xk (Newton-Verfahren) konvergiert die Folge quadratisch
(3.21) ‖xk+1 − x‖ ≤ l‖xk − x‖2, k = 0, 1, . . . .
Bemerkung:
(i) (3.18) bedeutet, dass r hinreichend klein sein muss; und x0 ∈ Br(x)bedeutet dann, dass die Ausgangsnaherung hinreichend genau ist.
(ii) Fur (3.17) ist die Lipschitzbedingung
‖f ′(z)− f ′(z)‖ ≤ l‖z − x‖, x, z ∈ BR(x)
hinreichend. Und diese ist erfullt, wenn z.B. f ∈ C2(D) ist.
(iii) Fur die quadratische Konvergenz des Newton-Verfahrens ist die Ein-fachheit der Nullstelle wesentlich.Denn z.B. f(x) = (x− x)2 (f(x) = 0, f ′(x) = 0, f ′′(x) 6= 0) ergibt dielineare Konvergenz xk+1 − x = 1
2(xk − x).
Beweis von Satz 3.2:
Fur festes z ∈ Br(x) sei
(1) gz(x) = x− f ′−1(z)f(x), x ∈ Br1(x), r1 =12
(r +R).
(∗)Ohne Einschrankung existiert f ′−1(z) auf BR(x) (Stetigkeit).
3 Nichtlineare Gleichungssysteme 49
Dann ist g′z(x) = 1 − f ′−1(z)f ′(x) = f ′−1(z)(f ′(z) − f ′(x)) fur x ∈ Br1(x).Daher gilt nach dem Mittelwertsatz fur x, y ∈ Br(x) (Br ⊂ Br1 ⊂ BR) miteinem 0 < θ < 1
‖gz(x)− gz(y)‖ ≤ ‖g′z(θx+ (1− θ)y)(x− y)‖≤ ‖g′z(...)‖‖x− y‖
(3.17)
≤ l‖z − (θx+ (1− θ)y)‖‖x− y‖= l‖θ(z − x) + (1− θ)(z − y)‖‖x− y‖≤ (lθ‖z − x‖+ l(1− θ)‖z − y‖)‖x− y‖≤ (lθ · 2r + l(1− θ)2r)‖x− y‖= 2lr‖x− y‖, (‖z − x‖ ≤ ‖z − x‖+ ‖x− x‖ ≤ 2r).
Daher gilt
(2) ‖gz(x)− gz(y)‖ ≤ γ‖x− y‖, x, y, z ∈ Br(x).
Wegen x = gz(x) ist fur x ∈ Br(x) und z ∈ Br(x)
‖gz(x)− x‖ = ‖gz(x)− gz(x)‖ ≤ γ‖x− x‖ ≤ r;
also gz(Br) ⊂ Br, z ∈ Br. Daher ist die Folge in (3.19) (xk+1 = gzk(xk)!)wohldefiniert und liegt in Br(x).
‖xk+1 − x‖ = ‖gzk(xk)− gzk(x)‖(2)
≤ γ‖xk − x‖ ≤ γk+1‖x0 − x‖,
woraus die lineare Konvergenz gegen x folgt.
Die quadratische Konvergenz (xk+1 = gxk(xk)) ergibt sich so:
‖xk+1 − x‖ = ‖gxk(xk)− gxk(x)‖= ‖g′
xk(θxk + (1− θ)x)(xk − x)‖
≤ ‖g′xk
(θxk + (1− θ)x)‖‖(xk − x)‖
mit einem 0 < θ(xk, x) < 1.Mit (3.17) folgt daher
‖xk+1 − x‖ ≤ l‖xk − (θxk + (1− θ)x)‖‖xk − x‖= l(1− θ)‖xk − x‖2
≤ l‖xk − x‖2
2
4 Interpolation, Integration und Differentiation 50
4 Interpolation, Integration und Differentiation
Im R2 besteht ein typisches Interpolationsproblem darin, durch n verschie-dene Punkte (xi, fi), i = 1, 2, . . . , n, eine geeignete Funktion zu legen, eineInterpolante Ln−1(x), mit der man die diskreten Informationen ”interpolie-ren“ bzw. ”extrapolieren“ kann.
x1 xn
••
••
Ln−1
x
y
(x1,f1)
(xn,fn)
Bei der Polynominterpolation verwendet man ein Polynom niedrigsten Gra-des, das diese Bedingungen erfullt. Es sei jedoch bemerkt, dass Polynomehohen Grades, d.h., wenn n >> 1, nicht geeignet sind zur Reproduktionvon Funktionen, die die Spuren (xi, fi) hinterlassen haben. Zur globalenApproximation durch Interpolation verwendet man stuckweise polynomiale(niederen Grades) Funktionen, Splines.
Die Approximationsgute durch Polynome auf einem beschrankten Intervallrichtet sich nach der globalen Glattheit der zu approximierenden Funkti-on. Lokale Storungen dieser Glattheit konnen mit Splines wesentlich besserberucksichtigt werden.
Da die Polynominterpolation die Grundlage fur viele numerische Verfahrenist (Integration, Differentiation, Differenzenverfahren,...), befassen wir unsdamit.
Da ein Polynom pn−1(x) =∑n−1
k=0 akxk n Koeffizienten besitzt (dimPn−1 =
dim < 1, x, . . . , xn−1 >= n) lautet das polynomiale Interpolationsproblem
Fur n verschiedene Stutzstellen x1, . . . , xn und Daten f1, . . . , fn istein pn−1 ∈ Pn−1 gesucht mit
(4.1) pn−1(xi) =n−1∑k=0
akxki = fi, i = 1, 2, . . . , n.
Dies ist offenbar ein lineares Gleichungssystem fur die ak, k = 0, . . . , n− 1.
4 Interpolation, Integration und Differentiation 51
Bekanntlich ist Ax = b (A ∈ Rn×n) eindeutig losbar, wenn A surjektiv ist,d.h. zu jedem b ∈ Rn ein x ∈ Rn existiert mit Ax = b. Dazu losen wir furfestes j ∈ 1, . . . , n das spezielle Problem
(4.2)n−1∑k=0
akxki = δij =
1, i = j0, i 6= j
.
Dieses Polynom ∈ Pn−1 soll also die n− 1 Nullstellen xi, i 6= j besitzen und
in xj den Wert 1 annehmen. Also lautet dieses Polynom lj(x) = cn∏i=1i 6=j
(x−xi),
bzw. wegen lj(xj) = 1
(4.3) lj(x) =
n∏i=1i6=j
(x− xi)
n∏i=1i6=j
(xj − xi).
Es ist also lj(xi) = δij , und daher lost das Lagrange-Polynom
(4.4) Ln−1(x) =n∑j=1
lj(x)fj
das Interpolationsproblem (4.1). Damit ist die auftretende Matrix in (4.1),die Vandermonde-Matrix A(x1, . . . , xn), nichtsingular (und die Koeffizien-ten von lj(x) bilden die j-te Spalte von A−1, s. (4.2)). Die lj(x) sind dieLagrange-Grundpolynome. Es ist naturlich (Eindeutigkeit)
(4.5) Pn−1 =< l1, . . . , ln > .
Mit dem Stutzstellenpolynom
(4.6) ω(x) =n∏i=1
(x− xi)
giltn∏i=1i 6=j
(xj − xi) = limx→xj
ω(x)(x− xj)
=ω′(xj)
1(L’Hospital),
und somit etwas handlicher
(4.7) lj(x) =ω(x)
(x− xj)ω′(xj)(x 6= xj).
4 Interpolation, Integration und Differentiation 52
Wir fassen zusammen und geben noch eine Fehlerabschatzung in
Satz 4.1 Es seien xi ∈ [a, b], i = 1, 2, . . . , n, n verschiedene Stutzstellen.
(i) Es gibt genau ein Polynom Ln−1 ∈ Pn−1, das die Interpolationsbeding-ungen (4.1) erfullt.
(ii) Die Lagrange-Darstellung lautet Ln−1(x) =∑n
k=1 lk(x)fk.
(iii) Fur f ∈ Cn[a, b], fi = f(xi), i = 1, 2, . . . , n, gilt fur x ∈ [a, b] dieFehlerdarstellung
(4.8) f(x)− Ln−1(f ;x) =f (n)(ξ)n!
ω(x)
mit einem ξ(x) mit minx, x1, . . . , xn < ξ < maxx, x1, . . . , xn.
Beweis: Es ist nur (4.8) zu zeigen fur x 6= xi, i = 1, . . . , n.Sei zusatzlich x fest. Fur t ∈ [a, b] sei
F (t) := f(t)−n∑k=1
lk(t)f(xk)− c ω(t),
wobei die Konstante c gleich festgelegt wird. Es ist, unabhangig von c,F (xi) = 0, i = 1, 2, . . . , n; und wir wahlen c so, dass F (x) = 0. D.h.
(1) c =f(x)− Ln−1(f ;x)
ω(x).
Daher besitzt F (n + 1) verschiedene Nullstellen x, x1, . . . , xn. Nach demSatz von Rolle besitzt F ′ n verschiedene Nullstellen in (minx, x1, . . . , xn,maxx, x1, . . . , xn), bzw. F (j), j = 1, . . . , n, n − j + 1 verschiedene Null-stellen dort. Also hat F (n) eine Nullstelle ξ in diesem Intervall.
F (n)(t) = f (n)(t)− 0− c ω(n)(t) = f (n)(t)− cn!.
Also ist c = f (n)(ξ)n! , woraus mit (1) die Behauptung folgt. 2
Fur die praktische Berechnung des Interpolationspolynoms ist die Lagrange-Darstellung zu aufwendig (und auch nicht sonderlich stabil). Z.B. muss beider Hinzunahme eines weiteren Punktes alles neu berechnet werden. DieDarstellung (4.4) hat viele theoretische Vorteile.
Die Berechnung uber (4.1) ist nicht zu empfehlen, da die Vandermonde-Matrix schon fur kleine n oft sehr schlecht konditioniert ist. Und hinzu
4 Interpolation, Integration und Differentiation 53
kommt, dass die Potenz-Darstellung, also die Verwendung der Basis 1, x, . . . , xn−1
auch oft sehr schlecht konditioniert ist: Sei a = (a0, . . . , an−1)t und
m = min‖a‖∞=1
∥∥∥∥∥n−1∑k=1
akxk
∥∥∥∥∥C[a,b]
, M = max‖a‖∞=1
∥∥∥∥∥n−1∑k=1
akxk
∥∥∥∥∥C[a,b]
(‖f(x)‖C[a,b] = maxx∈[a,b] |f(x)|), so ist die Kondition
cond(1, x, . . . , xn−1; [a, b]) :=M
m.
Denn es gilt ‖P
∆akxk‖
‖Pakxk‖
≤ cond‖∆a‖‖a‖ (Ubung).
Fur die Lagrange-Basis gilt cond = maxx∈[a,b]
∑n−1k=1 |lk(x)|, die – wenn man
die freie Wahl hat – man durch Wahl der xi gunstig beeinflussen kann(Ubung).
Mochte man Ln−1 an wenigen Stellen x auswerten, so ist es nicht notig,das Polynom auszurechnen. Grundlage dafur ist das folgende Lemma, wo-bei Lk(x1, . . . , xk+1;x) das Interpolationspolynom bezeichnet, das in xi dieWerte fi, i = 1, . . . , k + 1 interpoliert.
Lemma 4.1 Es gilt fur k ≥ 1
(4.9)Lk(x1, . . . , xk+1;x) = x−x1
xk+1−x1Lk−1(x2, . . . , xk+1;x)
+ xk+1−xxk+1−x1
Lk−1(x1, . . . , xk;x).
Beweis: Sei qk ∈ Pk das Polynom auf der rechten Seite von (4.9). Offenbarist qk(xi) = fi, i = 1, 2, . . . , k + 1, woraus wegen der Eindeutigkeit (4.9)folgt. 2
Setzt man fur jedes k = 1, 2, . . . , n
(4.10) Lik = Lk−1(xi−k+1, . . . , xi;x), i = k, . . . , n
(Li1 = L0(xi) = fi, i = 1, 2, . . . , n), so gilt mit Lemma 4.1 der Aitken-Neville-Algorithmus
Lik =x− xi−k+1
xi − xi−k+1Li,k−1 +
xi − xxi − xi−k+1
Li−1,k−1,(4.11)
k = 1, . . . , n, i = k, . . . , n
mit Lnn = Ln−1(x1, . . . , xn;x).
L11... L22...
. . .Ln1 · · · · · · Lnn
Dieses Schema bautman sich spaltenweisevon unten nach obenauf und uberschreibtdie erste Spalte.
4 Interpolation, Integration und Differentiation 54
Der Aufwand ist ∼ n2.Die fur die Numerik wichtigste Darstellung des Interpolationspolynoms istdie Newton-Darstellung : Sei
(4.12) ω0(x) = 1, ωk(x) =k∏i=1
(x− xi), k = 1, 2, . . . .
Was muss man zu Lk−1(x1, . . . , xk;x) addieren, um das PolynomLk(x1, . . . , xk, xk+1;x) zu erhalten? (Man kann dies unmittelbar aus (4.9)ablesen.) Es ist fur k ≥ 1
Lk(x1, . . . , xk, xk+1;xi)− Lk−1(x1, . . . , xk, ;xi) = 0, i = 1, 2, . . . , k;
oder
(4.13) Lk(x1, . . . , xk+1;x) = Lk−1(x1, . . . , xk;x) + ckωk(x), k ≥ 1
mit gewissen eindeutigen ck = ck(x1, . . . , xk+1, f1, . . . , fk+1).Also gilt (Teleskopreihe)
n−1∑k=1
(Lk(x1, . . . , xk+1;x)− Lk−1(x1, . . . , xk;x))
= Ln−1(x1, . . . , xn;x)− L0(x1;x)
=n−1∑k=1
ckωk(x);
oder mit
(4.14) c0 = L0(x1;x) = f1
erhalten wir die Newton-Darstellung
(4.15) Ln−1(x) = Ln−1(x1, . . . , xn;x) =n−1∑k=0
ckωk(x)
Definition 4.1 Fur die ck in (4.15) schreiben wir
[x1, . . . , xk+1]f = ck(x1, . . . , xk+1, f1, . . . , fk+1), k ≥ 1
[x1]f = f1.
Bei Kenntnis der ck lasst sich (4.15) schnell (∼ n) auswerten (Horner-Schema): n = 4
L3(x) = c0 + (x− x1)(c1 + (x− x2)(c2 + c3(x− x3)))
4 Interpolation, Integration und Differentiation 55
bzw. allgemein
(4.16)d = cn−1
j = n− 2, . . . , 0d = cj + d(x− xj+1).
Wie die Lk sind naturlich auch die ck rekursiv berechenbar. Aus (4.13)folgt, dass Lk(x1, . . . , xk+1;x) = ckx
k + . . . gilt; es ist sogar (s. (4.4)) ck =∑k+1i=1
fiω′k+1(xi)
. Vergleicht man daher die Koeffizienten von xk in (4.9), sofolgt die Rekursionsformel
(4.17)[x1, . . . , xk+1]f = [x2,...,xk+1]f−[x1,...,xk]f
xk+1−x1, k ≥ 1
[xi]f = fi.
Deswegen heißt [x1, . . . , xk+1]f die k-te dividierte Differenz der Daten fibzgl. xi. Und sie hat auch mit der k-ten Ableitung zu tun, falls fi = f(xi)und f ∈ Ck. Denn setzt man in (4.13) x = xk+1, so folgt mit der Darstellung(4.8)
f(xk+1)− Lk−1(x1, . . . , xk;xk+1) = [x1, . . . , xk+1]fωk(xk+1)
=f (k)(ξ)k!
ωk(xk+1);
oder
(4.18) k![x1, . . . , xk+1]f = f (k)(ξ).
Fur n = 4 sieht das Steigungsschema nach (4.17) so aus:
x1 f1 = c0
x2 f2 [x1, x2]f = c1
x3 f3 [x2, x3]f [x1, x2, x3]f = c2
x4 f4 [x3, x4]f [x2, x3, x4]f [x1, x2, x3, x4]f = c3
Dies baut man sich spaltenweise von unten nach oben auf.
Wie eingangs erwahnt, ist die Approximation von Funktionen durch Poly-nome hohen Grades nur in Spezialfallen sinnvoll. Selbst die C∞-Funktionf(x) = 1
1+x2 , x ∈ [−5, 5] (C. Runge 1901) kann nicht durch Interpolation an
den Knoten x(n)i = −5 + 10(i−1)
n−1 , i = 1, 2, . . . , n;n = 2, 3, . . . approximiertwerden (s. Bild).
Wir besprechen noch kurz die Idee der Spline-Funktionen, bei denen manstuckweise mit Polynomen eines festen Grades interpoliert. Wir behandelnnur den einfachsten Fall, namlich stuckweise linear und stetig (Polygon-Streckenzug). Sei
Th = xi |xi+1 = xi + hi, i = 0, 1, . . . , N − 1, x0 = a, xN = b
4 Interpolation, Integration und Differentiation 56
eine Zerlegung des Intervalls [a, b] und eine Funktion f : [a, b]→ R gegeben.
Es sei die stuckweise lineare Interpolation Lh(f, t) definiert durch
(4.19) Lh(f ;x) = L1(ti, ti+1;x), xi ≤ x ≤ xi+1, i = 0, 1, . . . , N − 1
a=x0 x1 x2 x3 x4 b=x5 x
•
•
•
•
••
y
f0
f1
f2
f3
f4f5
fi=f(xi)
Ist z.B. f ∈ C2[a, b], so folgt sofort aus der Fehlerdarstellung (4.8) fur x ∈[xi, xi+1] die Abschatzung
|f(x)− Lh(f ;x)| ≤ 12
maxx∈[xi,xi+1]
|f (2)(x)| |(x− xi)(x− xi+1)|
≤ 18h2i maxx∈[xi,xi+1]
|f (2)(x)|,
bzw. mit h = maxi hi
(4.20) ‖f − Lh(f)‖∞ ≤18h2‖f (2)‖∞.
Lh(f ;x) ist die Interpolante von f aus dem Splineraum
(4.21) S2,Th = s ∈ C[a, b] | s|(xi,xi+1) ∈ P1, i = 0, 1, . . . , N − 1.
Offenbar ist dimS2,Th = 2 ·Anzahl der Intervalle−Anzahl derStetigkeitsbedingungen in den inneren Knoten = 2N − (N − 1) = N + 1.Nicht jede Basis von S2,Th ist – wie bei Polynomen – gut konditioniert. Undes ist das besondere von Splineraumen, z.B.
Sk,Th = s ∈ Ck−2[a, b] | s|(xi,xi+1) ∈ Pk−1, i = 0, . . . , N − 1 (k ≥ 2),
dass es gut konditionierte Basen gibt, die auch noch lokalisiert sind, i.e.B-Splines. Das letztere ist von Bedeutung fur die ”richtige“ Wiedergabelokaler Eigenschaften, die man von dieser stuckweisen Konstruktion erwar-tet. Wir konstruieren jetzt die B-Spline Basis fur den einfachsten Fall undverweisen sonst auf [8]; hier auch auf [3].
Jedes s ∈ S2,Th ist festgelegt durch die Werte yi = s(xi), i = 0, 1, . . . , N .
4 Interpolation, Integration und Differentiation 57
x0=a xi xi+1 xN=b
•
•
••
yi
yi+1 s(x)
Sei 1 ≤ i ≤ N − 2 und x ∈ [xi, xi+1] (a < xi, xi+1 < b)
s(x) = yixi+1 − xxi+1 − xi
+ yi+1x− xixi+1 − xi
1 • • 1
xi−1 xi xi+1 xi+2
y y
Die beiden Geradenstucke und erganzen wir jeweils zu einerFunktion Hi bzw. Hi+1 ∈ S2,Th mit kleinstem Trager TrHi, wobei TrHi =x |Hi(x) 6= 0.Man erhalt die Hutfunktionen.
Hi(x) • • Hi+1(x)
xi−1 xi xi+1 xi+2
(4.22) Hi(x) =
x− xi−1xi − xi−1
, xi−1 ≤ x ≤ xi,xi+1 − xxi+1 − xi , xi ≤ x ≤ xi+1, i = 1, 2, . . . , N − 1
0, sonst
Fur die Randintervalle [x0, x1], [xN−1, xN ] machen wir eine ahnliche Kon-struktion
H0(x) (/∈C(R))
x0 x1
•HN (x)
xN−1 xN
•
4 Interpolation, Integration und Differentiation 58
Es gilt dann fur x ∈ [xi, xi+1], i = 0, 1, . . . , N − 1
(4.23) s(x) = yiHi(x) + yi+1Hi+1(x).
Daraus erhalten wir aber eine Darstellung, die nicht vom Intervall abhangt:Sei namlich x ∈ [a, b], dann existiert ein i ∈ 0, . . . , N−1 mit x ∈ [xi, xi+1].Wegen Hj(x) = 0, x /∈ TrHj = [xj−1, xj+1] ist daher fur x ∈ [xi, xi+1]Hj(x) = 0 fur j ≤ i− 1 und j ≥ i+ 2. Also gilt mit (4.23) die Darstellung
s(x) =i−1∑j=0
yjHj(x) + yiHi(x) + yi+1Hi+1(x) +N∑
j=i+2
yjHj(x),
(4.24) s(x) =N∑j=0
yjHj(x), x ∈ [a, b]
Die Hutfunktionen Hj , j = 0, . . . , N bilden die sogenannte B-Spline Basisfur S2,Th . Offenbar gilt Hj(xi) = δji und Lh(f ;x) =
∑Ni=0 f(xi)Hi(x).
Ferner besitzt diese Basis eine optimale Kondition, denn es gilt mit s(x) =∑Ni=0 aiHi, s(x) =
∑Ni=0(ai + ∆ai)Hi
(4.25)‖s− s‖∞‖s‖∞
=‖∆a‖∞‖a‖∞
(cond(H0, . . . ,HN , [a, b]) = 1).
Da die Basisfunktionen kompakten Trager besitzen, hat man die Moglich-keit, lokale Eigenschaften von Funktionen gut wiederzugeben.Ist f z.B. ein Signal, von dem nur die Werte in den Knoten abgetastet wer-den konnen, so ubertragt die Polynominterpolation LN (f ;x) =
∑Ni=0 fili(x)
=∑R−1
i=L+1 fili(x) ein Signal, das nicht lokalisiert ist; TrLN = [a, b].Aber Lh(f ;x) =
∑R−1i=L+1 fiHi(x) ist lokalisiert und ubertragt das Signal viel
besser; TrLh = [L,R].
a=x0 L R xN=b
f
Als Anwendung der Polynominterpolation befassen wir uns noch mit nume-rischer Integration und Differentiation.Oft kann man ein Riemann-Integral
∫ ba f(x) dx nicht elementar ausrechnen
(z.B.∫ 100−10 e
−x2dx). Man verwendet dann ein Quadraturverfahren der Form
(es gibt auch andere)
(4.26) Qn(f) =n∑k=1
Akf(xk)
4 Interpolation, Integration und Differentiation 59
mit den verschiedenen Stutzstellen xi und Gewichten Ai, i = 1, 2, . . . Beigegebenen xk entstanden die ersten Quadraturverfahren aus der Forderung
(4.27) Qn(f) =∫ b
af dx, f = 1, x, . . . , xn−1,
d.h. Qn(f) =∫ ba f dx, f ∈ Pn−1. (4.27) ist ein lineares Gleichungssystem fur
die Ak, was man sofort losen kann:(4.27) ⇔ Qn(li) =
∫ ba li(x) dx, k = 1, 2, . . . , n, Qn(li) = Ai. Daher ist (4.27)
gleichbedeutend mit der Wahl Ak =∫ ba lk(x) dx, d.h.
(4.28) Qn(f) =∫ b
aLn−1(f ;x) dx.
Wegen dieser Beziehung heißen die Formeln, die durch die Polynomgenau-igkeit (4.27) (eindeutig) definiert sind interpolatorische Quadraturformeln.
Wie bei der Polynominterpolation ist es nur in Ausnahmefallen sinnvoll, mitgroßen n zu arbeiten. Stattdessen bilden wir die Zerlegung a = t0 < t1 <. . . < tN = b, hi = ti+1 − ti, i = 0, 1, . . . , N − 1, die dem Verlauf von fangepasst ist. Es ist
(4.29)∫ b
af dx =
N−1∑i=0
∫ ti+1
ti
f dx.
Jedes∫ ti+1
tif dx wird mit einem Qn (n fest) ausgewertet, das jedoch von
i abhangt. Deswegen ist es hier zweckmaßig, sich auf ein Referenzintervall[−1, 1] zuruckzuziehen, um dann auf [ti, ti+1] zu transformieren mit
(4.30) x =12
(ti + ti+1) + thi2
= tmi + thi2, t ∈ [−1, 1]
Sei daher Qn eine interpolatorische Quadraturformel auf [−1, 1].
(4.31)∫ 1
−1g dt = Qn(g) +R0(g) =
n∑k=1
Akg(xk) +R0(g)
mit dem Restglied R0(g); R0(g) = 0, g ∈ Pn−1;R0(g) =? fur g ∈ Cn[−1, 1].Es ist
(4.32)
∫ ti+1
ti
f(x) dx =hi2
∫ 1
−1f(tmi + t
hi2
) dt
(4.31)=
hi2
n∑k=1
Akf(tmi + xkhi2
) +hi2R0(f(tmi + t
hi2
)).
Damit haben wir die interpolatorischen Quadraturformel fur [ti, ti+1] mitden Knoten tmi + xk
hi2 , k = 1, 2, . . . , n gefunden. Denn ist f ∈ Pn−1 ⇒
f(tmi + thi2 ) ∈ Pn−1 ⇒ R0(f(. . .)) = 0.
4 Interpolation, Integration und Differentiation 60
Summiert man jetzt gemaß (4.29) auf, so erhalt man die zusammengesetzteQuadraturformel Qh(f), die stuckweise Polynome vom Grad n − 1 exaktintegriert:
(4.33) Qh(f) =12
N−1∑i=0
hi
n∑k=1
Akf(tki + xkhi2
).
Hier gibt es kein Konvergenzproblem fur ”h→ 0“, denn
Qh(f) =12
n∑k=1
Ak
N−1∑i=0
f(tmi + xkhi2
)hi →12
n∑k=1
Ak
∫ b
af(x) dx =
∫ b
af dx
Als Beispiel behandeln wir die Trapezregel : n = 2
g
−1 0 1 t
L1(g;t)
x1 = −1, x2 = 1
QT (g) =1∫−1
L1(g; t) dt
= g(−1) + g(1)
Da der Quadraturfehler der integrierte Interpolationsfehler ist,
RT (g) =∫ 1
−1g dt−QT (g) =
∫ 1
−1(g(t)− L1(g; t)) dt,
erhalten wir mit der Fehlerdarstellung (4.8) fur g ∈ C2[−1, 1]
RT (g) =12
∫ 1
−1g(2)(ξ(t))(t2 − 1) dt =
12g(2)(η)
∫ 1
−1(t2 − 1) dt;
oder
(4.34)∫ 1
−1g(t) dt = QT (g)− 2
3g(2)(η), η ∈ (−1, 1), g ∈ C2[−1, 1].
Ist nun f ∈ C2[ti, ti+1], i = 0, 1, . . . , N − 1, so folgt (s. (4.32))
R0(f(tmi + thi2
)) = −23
(hi2
)2
f (2)(tmi + ηhi2
) = − 212h2i f
(2)(ηi)
mit einem ηi ∈ (ti, ti+1).
Fur das summierte Restglied Rh(f) =∑N−1
i=0hi2 R0(f(tmi + thi2 )),
∫ ba f dx =
Qh(f) +Rh(f), erhalten wir daher
(4.35’) RTh (f) = − 112
N−1∑i=0
hih2i f
(2)(ηi), f ∈ C2[ti, ti+1], i = 0, 1, . . . , N − 1
4 Interpolation, Integration und Differentiation 61
(4.35) |RTh (f)| ≤ 112
(b− a) maxi=0,1,...,N−1
(h2i ‖f (2)‖∞,[ti,ti+1]),
und es bietet sich eine Strategie zur Wahl der Schrittweite hi an.
Die summierte Formel lautet
(4.36) QTh (f) =12
(h0f(a) + hN−1f(b)) +12
N−1∑i=1
(hi + hi−1)f(ti),
Im aquidistanten Fall h = hi = b−aN ergibt sich die Euler-McLaurin-Summen-
formel , f ∈ C2k[a, b]
(4.37)∫ b
af(x) dx = QTh (f)+
k−1∑j=1
cj(h
2)2j(f (2j−1)(b)−f (2j−1)(a))+O(h2k),
die viele Anwendungen hat (Ubung).
Zum Abschluss befassen wir uns noch kurz mit numerischer Differentiation,wobei wir uns auf aquidistante Gitter beschranken. Analog zur Quadraturentstehen die gangigen Formeln durch Differentiation von Interpolationspo-lynomen. Es sei t fest, r, s ≥ 0 und (h > 0)
(4.38) t+ kh, k = −r, . . . , st−rh t t+sh
h︷︸︸︷.......
r + s+ 1 aquidistante Punkte.
Man mochte y(n)(t+ θh), t+ θh ∈ [t− rh, t+ sh] approximieren durch einenfiniten Ausdruck der Form
(4.39) D(n)h y(t) =
s∑k=−r
Ck(t, h)y(t+ kh).
Man nennt w = r+ s die Weite der Formel. Da es w+ 1 Koeffizienten gibt,bestimmen wir Ck aus der Forderung
(4.40) D(n)h y(t) = y(n)(t+ θh), y ∈ Pw.
Sind lk(x) die Lagrange-Grundpolynome bzgl. der Knoten−r, . . . , 0, 1, . . . , s,d.h.
(4.41) lk(x) =ω(x)
(x− k)ω′(k), k = −r, . . . , s, ω(x) =
s∏k=−r
(x− k),
4 Interpolation, Integration und Differentiation 62
so lauten die Lagrange Grundpolynome bzgl. der Knoten t+kh, k =−r, . . . , soffenbar
(4.42) lk(z) = lk
(z − th
), k = −r, . . . , s.
Satz 4.2 (a) (i) (analog zur Quadratur) Es gilt (4.40) genau dann, wenn
(4.43) Ck = h−nAk, Ak = l(n)k (θ), k = −r, . . . , s.
(ii) Es gilt (4.40) auch fur y ∈ Pw+1(†) genau dann, wenn (naturlich
(4.43) und)
(4.44) ω(n)(θ) = 0.
(b) Es sei w ≥ n (s. Bemerkung), und (4.40) gelte fur y ∈ Pq mit q ∈w,w + 1. Dann ist fur y ∈ Cq+1[t− rh, t+ sh]
(4.45) y(n)(t+ θh) = D(n)h y(t) +O(hq+1−n‖y(q+1)‖∞;[t−rh,t+sh]).
Beweis:
(a) (i) Pw =< l−r(z), . . . , ls(z) >. Daher gilt (4.40) genau dann, wennCk = l
(n)k (t+ θh), k = −r, . . . , s.
(ii) Fur p ∈ Pw+1 gilt die Zerlegung p(x) = cω(x) + q(x) mit q ∈Pw, c 6= 0 o.E.Der Rest ist Ubung.
(b) (fur spater) Nach Teil (a) gilt (h = 1, t = 0)
(1)s∑
k=−rAkg(k) = g(n)(θ), g ∈ Pq.
Wir suchen eine Darstellung fur R(g) = g(n)(θ) −∑s
k=−r Akg(k)
(= ( ddx)n(g(x)−
∑sk=−r lk(x)g(k))|x=θ
(3.8)= ( d
dx)n(g(w)(ξ(x))w! ω(x)|x=θ)
Sei g ∈ Cq+1. Nach der Taylorformel mit Integralrestglied gilt
(2) g(x) =q∑j=0
g(j)(θ)j!
(x−θ)j+∫ x
θ
(x− t)q
q!g(q+1)(t) dt = pq(x)+r(x)
(†)Fur y ∈ Pw+2 kann (4.40) nicht gelten
4 Interpolation, Integration und Differentiation 63
mit dem Restglied r(x), fur das deswegen gilt
(3) r(j)(θ) = 0, j = 0, . . . , q.
Daher ist
R(g) = R(Pq) +R(r)
= P (n)q (θ)−
s∑k=−r
AkPq(k) + r(n)(θ)−s∑
k=−rAkr(k)
(1),(3)= −
s∑k=−r
Akr(k).
Also gilt wegen r(x) =∫ sθ
(x−t)q+q! g(q+1)(t) dt
R(g) =∫ s
θ
(− 1q!
s∑k=−r
Ak(k − t)q+
)g(q+1)(t) dt(4)
=:∫ s
θK(t)g(q+1)(t) dt;
und somit
(5) |R(g)| ≤∫ s
θ|K(t)| dt max
t∈[−r,s]|g(q+1)(t)|.
Sei jetzt y ∈ Cq+1[t − rh, t + sh] und g(x) := y(t + xh), x ∈ [−r, s].Dann gilt
|R(g)| = |hny(n)(t+ θh)−s∑
k=−rAky(t+ kh)|
(5)
≤ Chq+1 maxx∈[t−rh,t+sh]
|y(q+1)(x)|,
woraus durch Division durch hn die Behauptung (4.45) folgt. 2
Bemerkung: Sinnvolle Formeln erfullen stets w ≥ n;z.B. y′(t) = y(t+h)−y(t)
h +O(h).Ist w < n in (4.43), so gilt Ck = 0, k = −r, . . . , s.
Wir geben noch einige Formeln an.Mit den Differenzen Dly(t) = D(Dl−1y(t)), Dy(t) = 1
h(y(t+h)−y(t)) erhaltman die folgende Approximation 2. Ordnung mit w = n
(4.46)
D2ky(t− kh) = y(2k)(t) +O(h2‖y(2k+2)‖),w = 2k, q = w + 1
D2k+1y(t− kh) = y(2k+1)(t+ h2 ) +O(h2‖y(2k+3)‖),
w = 2k + 1, q = w + 1
4 Interpolation, Integration und Differentiation 64
Formeln mit w > n sind z.B.
(4.47)
12h(y(t+ h)− y(t− h))
= y′(t) +O(h2‖y(3)‖), w = 2, q = w
112h(y(t− 2h)− 8y(t− h) + 8y(t+ h)− y(t+ 2h))
= y′(t) +O(h4‖y(5)‖), w = 4, q = w
112h(−y(t− 2h) + 16y(t− h)− 30y(t) + 16y(t+ h)− y(t+ 2h))
= y′′(t) +O(h4‖y(6)‖), w = 4, q = w + 1
1h2 (2y(t)− 5y(t+ h) + 4y(t+ 2h)− y(t+ 3h))
= y′′(t) +O(h2‖y(4)‖), w = 3, q = w.
5 Gewohnliche und partielle DGL, numerische Methoden 65
5 Gewohnliche und (wenige) partielle Differential-gleichungen, numerische Methoden
Viele Prozesse der Informationsverarbeitung werden durch Differentialglei-chungen modelliert, z.B. die Antwort auf: Wie schnell wachst Fußpilz?Ist n(t), t ≥ t0, die Zahl der Lebewesen zur Zeit t, so nehmen wir n(t) ∈ C1
an. Dann ist n′(t) die Wachstumsgeschwindigkeit zur Zeit t. Wir machenjetzt die Modellannahme (siehe (1) Kap. 1)
(5.1) n′(t) = αn(t), t ≥ t0, α > 0;
und man hat die Idee, dass die Dynamik des Wachstums ausgehend vonn(t0) durch (5.1) bestimmt ist.(5.1) ist das erste Beispiel einer Differentialgleichung 1. Ordnung . Gibt mann(t0) vor, so spricht man von einem Anfangswertproblem 1. Ordnung . Es istklar, dass durch (5.1) n(t) nicht allein bestimmbar ist, denn mit n(t) wareauch 7 · n(t) eine Losung.
Die Gleichung (5.1) kann man hier noch nach n(t) ”auflosen“, was i.a. nichtexplizit moglich ist.Wir nehmen an, (5.1) ist losbar mit einem n(t) > 0, t ≥ t0. Dann istddt log n(t) = α, t ≥ t0, oder log n(t) = αt+ c, oder
(5.2) n(t) = ceαt, t ≥ t0.
Man sieht hier leicht, dass das Anfangswertproblem (n(t0) = n0) durch
(5.3) n(t) = n(t0)eα(t−t0), t ≥ t0,
eindeutig losbar wird.
Dieses Modell hat die Eigenschaft limt→∞ n(t) = ∞, was nicht realistischist, denn die Wachstumsrate n′(t)
n(t) ist konstant (100α% pro sek). Man kannein besseres Modell formulieren (Ubung), das wieder auf ein komplizierteresProblem der Form
(5.4)y′(t) = f(t, y(t)), t ≥ t0y(t0) = y0 geg.
(AWP)
fuhrt.
Man kann sich auch vorstellen, dass mehrere Sorten von Lebewesen siedelnund sich bekriegen. Man erhalt z.B. ein Rauber (x(t)) – Beute (y(t)) Modell,
5 Gewohnliche und partielle DGL, numerische Methoden 66
wo die Wachstumsraten jeweils durch die andere Sorte beeinflusst wird(∗):
(5.5)
x′(t) = (α− β y(t)) x(t),y′(t) = (−γ + δ x(t)) y(t)
t ≥ t0
mit α, β, γ, δ ≥ 0 und x(t0), y(t0) gegeben.
Dies ist wieder von der Form (5.4), wobei f und y vektorwertig sind. Manspricht von einem Differentialgleichungssystem 1. Ordnung. Es ist keines-wegs klar, dass das AWP (5.5) uberhaupt losbar ist, und die Losung furt ≥ t0 existiert! Deswegen mussen wir eine Existenztheorie entwickeln undnumerische Verfahren, um diese dann wenigstens naherungsweise zu ermit-teln.
Es konnen auf naturliche Weise auch große Systeme (5.4), d.h. f ∈ Rn,n 1, entstehen, z.B. bei der Beschreibung von Vielteilchensystemen oderbei der folgenden (teilweisen) Diskretisierung einer partiellen Differential-gleichung .Dies ist eine Differentialgleichung fur eine Funktion von mehreren Verander-lichen; z.B. u(x, t), x ∈ R, t ∈ R. Gegeben sei ein (eindimensional beschreib-barer) Stab der Lange L, der zur Zeit t = 0 die Temperaturverteilungϕ0(x), 0 ≤ x ≤ L besitzt mit z.B. ϕ0(0) = 0, ϕ0(L) = 0.
0 L
t
x
←− ϕ0(x)
←− ϕ2(t)ϕ1(t) −→
yu(x, t)
Die Enden des Stabes stecken in einem Warmebad mit den Temperaturenϕ1(t) bzw. ϕ2(t) mit naturlich ϕ1(0) = ϕ2(0) = 0. Gesucht ist die Tempera-tur im Stab an der Stelle x, 0 < x < L, zur Zeit t, 0 < t ≤ T ; u(x, t).
Dies fuhrt auf die Warmeleitungsgleichung bzw. auf das partielle Anfangs-randwertproblem fur u(x, t)
(5.6)
∂u
∂t(x, t) =
∂2u(x, t)∂x2
, 0 < x < L, 0 < t ≤ T
mit u(x, 0) = ϕ0(x), 0 ≤ x ≤ L Anfangsbedingung
u(0, t) = ϕ1(t), 0 ≤ t ≤ T Randbedingungenu(L, t) = ϕ2(t), 0 ≤ t ≤ T
(∗) Ist z.B. y(t0) = 0⇒ y(t) ≡ 0⇒ x(t)→∞. x(t0) = 0 ergibt y(t)→ 0.
5 Gewohnliche und partielle DGL, numerische Methoden 67
Wir erhalten eine ”Linienmethode“, durch ”Diskretisierung“, der Ortsvaria-blen x:
Sei h = LN , xi = ih, i = 0, . . . , N . Es ist nach der Differentiationsformel
(4.47), k = 1
(5.7)∂2u
∂x2(ih, t) =
u((i− 1)h, t)− 2u(ih, t) + u((i+ 1)h, t)h2
+O(h2)
falls u hinreichend glatt ist.Daher ist mit (5.6)
(5.8)∂u
∂t(ih, t) =
u((i− 1)h, t)− 2u(ih, t) + u((i+ 1)h, t)h2
+ Fehlerterm
Das numerische Verfahren entsteht durch Streichen des Fehlerterms:
(5.9) u′i(t) =ui−1(t)− 2ui(t) + ui+1(t)
h2, i = 1, . . . , N − 1.
Und wir hoffen, dass ui(t) eine Approximation fur u(ih, t) ist. Wir setzendaher
(5.10)u0(t) = ϕ1(t),uN (t) = ϕ2(t),ui(0) = ϕ0(ih),
0 ≤ t ≤ T ;
i = 1, . . . , N − 1
Mit (5.9) erhalten wir daher das große System erster Ordnung fury(t) = (u1(t), . . . , uN−1(t))
(5.11) y′(t) =1h2
−2 11 −2 1
. . .1 −2
y(t)+1h2
ϕ1(t)
0...0
ϕ2(t)
, 0 ≤ t ≤ T
mit y(0) = (ϕ0(h), . . . , ϕ0((N − 1)h))t.
Probleme der Form (5.6) tauchen auch in der Bildverarbeitung auf.
Man muss sich auch mit Differentialgleichungen hoherer Ordnung befassen:
(5.12) y(m)(x) = f(x, y(x), . . . , y(m−1)(x)),
wobei sogar Systeme auftreten, d.h. f ∈ Rn, y ∈ Rn, so dass man dann einSystem von n Differentialgleichungen m-ter Ordnung hat.
5 Gewohnliche und partielle DGL, numerische Methoden 68
Es ist nun ganz wichtig fur Theorie und Praxis, dass man stets o.E. m = 1annehmen kann.
Beispiel: Sei I ein Intervall, y ∈ C2(I) sei Losung der Differentialgleichung2. Ordnung auf I
(1) y′′(x) = f(x, y(x), y′(x)), x ∈ I.
Setze z(x) =(y(x)y′(x)
)=(z1(x)z2(x)
). Dann gilt
z′(x) =(y′(x)y′′(x)
)=
(z′1(x)z′2(x)
)=
(y′(x)
f(x, y(x), y′(x))
)
=(
z2(x)f(x, z1(x), z2(x))
)=: F (x, z(x))
Also ist z(x) ∈ C1(I) Losung des Systems 1. Ordnung
(2) z′(x) = F (x, z(x)), x ∈ I.
Ist umgekehrt z ∈ C1(I) Losung von (2), so ist z2(x) = z′1(x). Setzt many(x) := z1(x), so ist z2(x) = y′(x) ∈ C1(I), d.h. y ∈ C2(I), und wegenz2(x) = f(x, z1(x), z2(x)) ist y′′(x) = f(x, y(x), y′(x)), x ∈ I.
Allgemein gilt
Lemma 5.1 Sei I ein Intervall.
(a) Sei y ∈ Cm(I) Losung der Differentialgleichung m-ter Ordnung
(5.13) y(m)(x) = f(x, y(x), . . . , y(m−1)(x)), x ∈ I.
Dann ist z(x) = (y(x), . . . , y(m−1)(x))t ∈ C1(I) Losung des Systems1. Ordnung
(5.14) z′(x) = F (x, z(x)) =
z2(x)
...zm(x)
f(x, z1(x), . . . , zm(x))
, x ∈ I.
(b) Ist z(x) ∈ C1(I) Losung des Systems (5.14), so ist y(x) := z1(x)Losung von (5.13), und es gilt zi(x) = y(i−1)(x), i = 1, . . . ,m, x ∈ I.
Beweis:
(a) ist klar.
5 Gewohnliche und partielle DGL, numerische Methoden 69
(b) Sei z(x) ∈ C1(I) Losung von (5.14), so ist mit y(x) = z1(x) wegenz′i(x) = zi+1(x), i = 1, . . . ,m− 1, zi(x) = y(i−1)(x), i = 2, . . . ,m. Unddie letzte Gleichung ergibt z′m(x) = y(m)(x) = f(x, z1(x), . . . , zm(x)) =f(x, y(x), . . . , y(m−1)(x)), x ∈ I 2
Die Funktion f in Lemma 5.1 kann naturlich auch vektorwertig sein.Wir beschranken uns daher auf Systeme 1. Ordnung.
Definition 5.1 Sei n ∈ N und D ⊂ R × Rm offen und sei f : D → Rm,
(x0, y0) ∈ D. I sei ein Intervall mitI 6= ∅ und x0 ∈ I. Dann heißt y ∈ C1(I)
Losung des Anfangswertproblems (AWP)
(5.15)y′(x) = f(x, y(x))y(x0) = y0
auf I, wenn gilt
1. y(x0) = y0
2. (x, y(x)) ∈ D,x ∈ I
3. y′(x) = f(x, y(x)), x ∈ I.x0
y0 •
D
Rm
xI
Wir verlangen also, dass y′(x) stetig auf I ist und die Differentialgleichungpunktweise fur x ∈ I erfullt ist. Dies ist unser Losungsbegriff (es gibt auchandere). Ist y eine Losung des AWP’s auf einem Intervall I mit x0 ∈
I , so
heißt y lokale Losung des AWP’s.
Eine lokale Losung y des AWP’s, y ∈ C1(Imax), heißt maximale Losung,wenn sie jede andere lokale Losung z, z ∈ C1(Iz) des AWP’s fortsetzt; d.h.Iz ⊂ Imax und z(x) = y(x), x ∈ Iz.
Eine maximale Losung ist naturlich eindeutig (falls sie existiert), und wirsind auf der Suche nach Bedingungen fur die Existenz einer maximalenLosung.
Fur das weitere brauchen wir noch eine nutzliche Ungleichung.Ist g : [a, b]→ Rm, g(x) = (g1(x), . . . , gm(x))t, so heißt g Riemann-integrier-bar , wenn gi ∈ R[a, b], i = 1, . . . ,m; und
∫ ba g(x) dx := (
∫ ba g1(x) dx, . . . ,∫ b
a gm(x) dx)t.Ist g ∈ R[a, b], so auch ‖g(x)‖∞ = maxi=1,...,m |gi(x)| ∈ R[a, b] (†).
(†)z.B. ist maxi=1,2 |gi(x)| = 12(g1(x) + g2(x)) + 1
2|g1(x)− g2(x)|
5 Gewohnliche und partielle DGL, numerische Methoden 70
Daher folgt∣∣∣∣∫ b
agi(x) dx
∣∣∣∣ ≤ ∫ b
a|gi(x)| dx ≤
∫ b
a‖gi(x)‖∞ dx oder∥∥∥∥∫ b
ag(x) dx
∥∥∥∥∞≤∫ b
a‖g(x)‖∞ dx.
Ein Ausgangspunkt fur die Existenz lokaler Losungen ist das folgende
Lemma 5.2 Gegeben sei das AWP (5.15) mit einer in D stetigen Funktionf . Dann ist eine auf dem Intervall I (
I 6= ∅) stetige Funktion y(x) Losung
des AWP’s auf I genau dann, wenn
(5.16) y(x) = y0 +∫ x
x0
f(t, y(t)) dt, x ∈ I,
gilt. D.h. genau dann, wenn die Funktion y ∈ C(I) die Integralgleichung(5.16), die eine abstrakte Fixpunktgleichung ”y = ϕ(y)“ im Raum C(I)darstellt, erfullt.
Beweis: ⇒: Ist y Losung des AWP’s, so folgt durch Integration (5.16).⇐: Es ist f(t, y(t)) ∈ C(I). Daher ist wegen (5.16) y ∈ C1(I) mit y′(x) =f(x, y(x)), x ∈ I, nach dem Fundamentalsatz. 2
Wir wollen uns jetzt mit dem Verlauf einer maximalen Losung befassen.Lauft sie von ”Rand zu Rand“, des Gebietes D?
Definition 5.2 Sei A ⊂ Rm. ∂A bezeichnet den Rand von A. x ∈ ∂A genaudann, wenn in jeder ε-Umgebung von x, Bε(x), mindestens ein Punkt ausA und mindestens ein Punkt, der nicht zu A gehort, liegt.
D.h. x ∈ ∂A⇔ fur jedes ε > 0 ist Bε(x)∩A 6= ∅und Bε(x) ∩ CA 6= ∅. Da x ∈ A genau dann,wenn Bε(x) ∩A 6= ∅ fur jedes ε > 0 gilt, ist
(5.17) ∂A = A ∩ CA;
oder wegen x ∈ CA genau dann, wenn x keininnerer Punkt, d.h. x /∈
A, hat man auch ∂A =
A \A.
A
x• ← Bε(x)
x
y
D
∂D
∂D
x
y
D
∂D
x
y
∂D
← ∂D
D
5 Gewohnliche und partielle DGL, numerische Methoden 71
Wir betrachten das folgende AWP mit einer unstetigen rechten Seite f :
(5.18) f(x, y) =
1, x ∈ R, 0 < y < 70, x ∈ R, y ≤ 0
, D = (x, y) |x ∈ R, y < 7
7
y0
y
x0x0−y0 x
•
D
∂D
Wir suchen die maximaleLosung mit dem Anfangswert(x0, y0) ∈ D.
Fall 1: 0 < y0 < 7: Es gilt y′(x) = 1, y(x0) = y0, also y(x) = y0 + (x− x0)solange y(x) > 0 ist; denn nur dann ist y′(x) = f(x, y(x)) = 1, alsonur fur x0 − y0 < x < x1 (y0 + (x1 − x0) = 7).
Diese lokale Losung y ∈ C1(I), I = (x0 − y0, x1) stoßt bei x1 gegenden Rand und ist aber nicht nach links fortsetzbar. Denn sie erfullt inx0 − y0 nicht die Differentialgleichung, day′(x0 − y0) = 1 6= f(x0 − y0, y(x0 − y0)) = f(x0 − y0, 0) = 0.
Unser Losungsbegriff verlangt, dass die Dgl. punktweise erfullt ist. Da-her ist das obige Geradenstuck die maximale Losung zum Anfangswert(x0, y0) ∈ D mit y0 > 0. Sie endet in D!
Fall 2: y0 ≤ 0: Dann ist offenbar y(x) = y0, x ∈ R die maximale Losungmit I = R. Sie lauft in diesem Sinne von ”Rand zu Rand“.
Ist f stetig in D, so werden wir sehen, dass die maximale Losung nicht inD enden kann. Wegen y ∈ C1(I) konnen Singularitaten der Losung nur inden Randpunkten von I auftreten:
(5.19)y′(x) = y2(x)y(0) = 1
D = (x, y) |x ∈ R, y > 0
Es ist y(x) = 11−x , x ∈ I =
(−∞, 1) die maximale Losung.Die Losung lauft von ”Rand zuRand“; es gilt limx↑1 y(x) =∞.
1 x
y
1D
∂D
5 Gewohnliche und partielle DGL, numerische Methoden 72
Es kann also vorkommen, dass bei einer maximalen Losung y ∈ C1(I),I = (a, b), limx↑b ‖y(x)‖ = ∞ gilt (dazu muss naturlich D unbeschranktsein).
Ist eine maximale Losung y ∈ C1(a, b) beschrankt, d.h. supx∈(a,b) ‖y(x)‖ <∞, so muss leider limx↑b y(x) nicht notwendig existierten, wie das folgendeBeispiel zeigt.
(5.20)y′(x) = − 1
x2(1− y2(x))
12
y(1) = sin(1)D = (0, 2)× (−1, 1).
Es ist namlich y(x) = sin 1x , x ∈ (0, 2) die maximale Losung, |y(x)| ≤ 1.
Aber limx↓0 y(x) existiert nicht!Es gibt namlich zu jedem y ∈ [−1, 1] eine Folge xk∞k=1, xk > 0 mitlimk→∞ xk = 0 und y = limk→∞ sin 1
xk, da wk = 1
xk→∞, sin 1
xk= y.
−1
1y
sinw
wwk
• • • • • •
Die Losung lauft von ”Rand zu Rand“ in dem Sinn, dass
limk→∞
(xk
y(xk)
)=(
0y
)∈ ∂D.
x
y
D
← ∂D
1
•
Wir formulieren jetzt ohne Beweis den Hauptsatz uber das Randverhalteneiner maximalen Losung, wenn die rechte Seite f stetig ist.
5 Gewohnliche und partielle DGL, numerische Methoden 73
Satz 5.1 Gegeben sei das AWP aus Definition 5.1 mit einem stetigen f . Esexistiere die maximale Losung y ∈ C1(I) mit einem (notwendigerweise(‡))offenen Intervall I. Es sei b ∈ R rechter Randpunkt von I. Sei ε > 0, sodass [b− ε, b) ⊂ I ist.
Dann gilt entweder
(a) ‖y(x)‖ ist auf [b−ε, b) unbeschrankt (limk→b
‖y(x)‖ =∞) d.h. es existierteine Folge xk∞k=1 aus I mit limk→∞ xk = b und limk→∞ ‖y(xk)‖ =∞(dazu muss D unbeschrankt sein)
oder
(b) ‖y(x)‖ ist auf [b− ε, b) beschrankt (limk→b
‖y(x)‖ <∞).Existiert dann limx→b y(x) = z, so ist (b, z) ∈ ∂D.Existiert limx→b y(x) nicht (s. (5.20)), so gibt es eine Folge xk∞k=1
aus I mit limk→∞ xk = b, limk→∞ y(xk) = z und (b, z) ∈ ∂D.
Entsprechendes Verhalten liegt naturlich auch in einem linken Randpunkta ∈ R vor, so dass die Losung in diesem Sinne von ”Rand zu Rand“ lauft.Dies meint man auch im Fall I = (a,∞) (a ∈ R), I = (−∞, b) (b ∈ R), I =(−∞,∞).
Wir kommen nun zu der Frage der Existenz und Eindeutigkeit von lokalenLosungen, was uns dann zur maximalen Losung fuhren wird. Der beruhmteExistenzsatz von Peano besagt:
Das AWP aus Definition 5.1 besitzt eine lokale Losung, wenn fstetig ist.
Der Beweis (s. Literatur) beruht auf einer Folge von ”Naherungslosungen“,die man sich mit den einfachsten numerischen Verfahren (Euler-Cauchy-Verfahren) erzeugt und dann zeigt, dass eine Teilfolge gegen eine lokaleLosung des AWP’s konvergiert.
Die Stetigkeit impliziert aber nicht die eindeutige Losbarkeit! Diese mussman fur die numerische Approximation naturlich haben.
(5.21) Beispiel:
y′(x) = |y(x)|α mit 0 < α < 1y(0) = 0
(‡)s. Satz 5.4
5 Gewohnliche und partielle DGL, numerische Methoden 74
und D = R × R, f(x, y) = |y|α. Es ist f stetig in D, und y(x) = 0 ist einelokale Losung mit I = R. Dies ist keine maximale Losung, denn fur b > 0ist z.B.
y(x) =
0, −∞ < x < b
(1− α)1
1−α (x− b)1
1−α , b ≤ x <∞
auch eine lokale Losung (die von der ersten abzweigt).
xy
b
y
Bemerkung: Man erhalt keine von y = 0 abzweigende Losung im Fallα ≥ 1.
Ist f stetig und hangt f nicht von y ab, so ist das AWP naturlich durchy0 +
∫ xx0f(t) dt eindeutig losbar. Eindeutigkeitsbedingungen sind daher Be-
dingungen an y in f(x, y). Im Beispiel (5.21) hat f = |y|α eine starke Sin-gularitat der Ableitung im Nullpunkt fur 0 < α < 1,
fy =
α|y|α−1, y > 0−α|y|α−1, y < 0
,
die fur α ≥ 1 nicht auftaucht.
Die folgende Lipschitzbedingung verlangt die Beschranktheit der Differen-zenquotienten bzgl. y und wird uns Eindeutigkeit bescheren.
Definition 5.3 Gegeben sei das AWP aus Definition 5.1.
(a) f genugt auf D einer Lipschitzbedingung, wenn eine Lipschitzkon-stante L > 0 existiert, so dass fur (x, y), (x, y∗) ∈ D gilt
(5.22) ‖f(x, y)− f(x, y∗)‖ ≤ L‖y − y∗‖.
(b) f genugt auf D einer lokalen Lipschitzbedingung, wenn es zu je-dem (x, y) ∈ D eine Umgebung Br(x, y) ⊂ D gibt, so dass f aufBr(x, y) die Lipschitzbedingung (5.22) mit einem L = L(Br(x, y))erfullt.
Mit dem Mittelwertsatz der Differentialrechnung zeigt man sofort (Br istkonvex):
5 Gewohnliche und partielle DGL, numerische Methoden 75
Lemma 5.3 Gegeben sei das AWP aus Definition 5.1. f besitze in D stetigepartielle Ableitungen ∂fi
∂yk(x, y), i, k = 1, . . . , n. Dann genugt f auf D einer
lokalen Lipschitzbedingung.
Wir werden nur die lokale Lipschitzbedingungen fordern, die naturlich we-sentlich schwacher als die globale ist. Z.B. erfullt f(x, y) = 1 + y2 keineLipschitzbedingung in D = R× R, aber eine lokale Lipschitzbedingung.
Ein wichtiges Hilfsmittel nicht nur zum Beweis der Eindeutigkeit ist dasfolgende Gronwall-Lemma:
Lemma 5.4 Es seien u, v ∈ C[a, b] mit u, v ≥ 0. Gilt mit einer Konstantenc ≥ 0 die Ungleichung
(5.23) u(x) ≤ c+∫ x
au(t)v(t) dt, x ∈ [a, b],
so folgt
(5.24) u(x) ≤ c eR xa v(t) dt, x ∈ [a, b].
Ist insbesondere c = 0 in (5.23), so folgt u = 0.
Beweis: Sei h(x) := c+∫ xa u(t)v(t) dt.
Ist c > 0, so gilt h(x) > 0, und es folgt h′(x) = u(x)v(x) ≤ h(x)v(x), bzw.ddx log h(x) ≤ v(x).
Also ist ∫ x
a
d
dtlog h(t) dt = log
h(x)h(a)
≤∫ x
av(t) dt,
woraus wegen h(a) = c (5.24) fur c > 0 folgt.
Ist jetzt c = 0 in (5.23), so gilt (5.23) fur beliebige c > 0 und somit auch(5.24) fur beliebige c > 0; woraus u = 0 folgt. 2
Satz 5.2 (Eindeutigkeit der Losung) Gegeben sei das AWP aus Defini-tion 5.1 mit einem stetigen f , das einer lokalen Lipschitzbedingung genugt.
(a) Sind yi ∈ C1(Ii), i = 1, 2, zwei lokale Losungen des AWP (x0 ∈I i,
i = 1, 2), dann gibt es ein offenes Intervall I ⊂ I1 ∩ I2 mit x0 ∈I und
y1(x) = y2(x), x ∈ I.
(b) Losungen zu verschiedenen Anfangswerten konnen sich nicht in Dschneiden. Insbesondere gilt in Teil (i) sogar y1(x) = y2(x), x ∈ I1∩I2.
5 Gewohnliche und partielle DGL, numerische Methoden 76
Bemerkung: Ist y ∈ Rn, n ≥ 2, so konnen sich in (ii) einzelne Komponen-ten von Losungen zu verschiedenen Anfangswerten durchaus schneiden (nurnicht alle in einem Punkt).
Beweis von Satz 5.2: : Wir brauche nur (i) zu zeigen. Nach Lemma 5.2gilt fur die lokalen Losungen yi ∈ C1(Ii),
(1) yi(x) = y0 +∫ x
xo
f(t, yi(t)) dt, x ∈ Ii, i = 1, 2.
Es existiert eine Kugel Br(x0, y0) ⊂ D, so dass f dort eine Lipschitzbedin-gung erfullt.Ferner gibt es wegen der Stetigkeit der yi(x) in x0, limx→x0 yi(x) = y0, einoffenes Intervall I ⊂ I1 ∩ I2, x0 ∈ I mit
(2) (xi, yi(x)) ∈ Br(x0, y0), x ∈ I, i = 1, 2.
Sei L die Lipschitzkonstante von f bzgl. Br(x0, y0). Fur x ∈ I, x ≥ x0 giltmit (1)
‖y1(x)− y2(x)‖ =∥∥∥∥∫ x
x0
(f(t, y1(t))− f(t, y2(t))) dt∥∥∥∥
≤∫ x
x0
‖f(t, y1(t))− f(t, y2(t))‖ dt
(2)
≤ L
∫ x
x0
‖y1(t)− y2(t)‖ dt.
Die letzte Ungleichung ergibt nun mit dem Gronwall Lemma (fur c = 0) dieAussage y1(x) = y2(x), x ≥ x0, x ∈ I. 2
Unter den Voraussetzungen von Satz 5.2 werden wir jetzt die Existenz einerlokalen Losung zeigen. Wir gehen dabei von der Fixpunktgleichung (5.16)aus, und zeigen die Existenz wie bei dem Banachschen Fixpunktsatz (Satz3.1).
Satz 5.3 (Picard-Lindelof ∼ 1893) Gegeben sei das AWP aus Definiti-on 5.1 mit einem stetigen f , das auf D eine lokale Lipschitzbedingung erfullt.Dann existiert ein Intervall Iα = [x0 − α, x0 + α], α > 0, so dass das AWPgenau eine Losung auf Iα besitzt.
Beweis (nur zur Information): Da D offen ist, gibt es ein Rechteck
(1) R = (x, y) | |x− x0| ≤ a, ‖y − y0‖ ≤ b ⊂ D,
so dass f dort eine Lipschitzbedingung mit der Konstanten L erfullt.Sei ferner
(2) M = max(x,y)∈R
‖f(x, y)‖
5 Gewohnliche und partielle DGL, numerische Methoden 77
und
(3) α = mina,
b
M,
17L
.
Wir zeigen, dass ein y ∈ C(Iα) existiert mit
(4) y(x) = y0 +∫ x
x0
f(t, y(t)) dt, x ∈ Iα.
Sei fur y ∈ C(Iα), ϕ(y;x) := y0 +∫ xx0f(t, y(t)) dt. Dann ist auch ϕ(x; y) ∈
C(Iα); und wir suchen mit (4) einen Fixpunkt von ϕ in C(Iα). Dazu gehenwir genauso vor wie im Banachschen Fixpunktsatz, Satz 3.1.
(5) G = y ∈ C(Iα) | ‖y(x)− y0‖ ≤ b, x ∈ Iα
ϕ bildet G in sich ab: Sei y ∈ G, dann ist
‖ϕ(y;x)− y0‖ =∥∥∥∥∫ x
x0
f(t, y(t)) dt∥∥∥∥ ≤ ∣∣∣∣∫ x
x0
‖f(t, y(t))‖ dt∣∣∣∣
≤ M |x− x0| ≤ Mα(3)
≤ Mb
M= b
fur x ∈ Iα (da (x, y(x)) ∈ R, x ∈ Iα).
ϕ ist kontrahierend auf G: Es seien y1, y2 ∈ G und x ∈ Iα. Dann gilt furx ≥ x0
‖ϕ(y1;x)− ϕ(y2;x)‖ =∥∥∥∥∫ x
x0
(f(t, y1(t))− f(t, y2(t)) dt∥∥∥∥
≤∫ x
x0
‖(f(t, y1(t))− f(t, y2(t))‖ dt
≤∫ x
x0
L‖y1(t)− y2(t)‖ dt
≤ L maxt∈Iα‖y1(t)− y2(t)‖ |x− x0|
= L‖y1 − y2‖∞;Iα |x− x0|
≤ Lα‖y1 − y2‖∞;Iα =: K‖y1 − y2‖∞;Iα
mit K ≤ 17 wegen (3).
Also gilt fur y1, y2 ∈ G
(6) ‖ϕ(y1)− ϕ(y2)‖∞;Iα ≤ K‖y1 − y2‖∞;Iα .
Wie in (3.8),(3.9) erhalten wir daher, dass die Folge yk(x)∞k=0 definiertdurch
(7) yk+1(x) = ϕ(yk;x), k = 0, 1, . . .
5 Gewohnliche und partielle DGL, numerische Methoden 78
das Cauchy-Kriterium fur glm. Konvergenz auf Iα erfullt: Zu jedem ε > 0existiert ein N(ε) mit
‖yn − ym‖∞;Iα < ε, m > n ≥ N(ε).
Also konvergiert diese Folge glm. gegen eine stetige Funktion y ∈ C(Iα).Wegen yk(x) ∈ G, ‖yk(x)−y0‖ ≤ b, x ∈ Iα, gilt auch y ∈ G, bzw. (x, y(x)) ∈R, x ∈ Iα.
Wegen‖f(t, yk(t))− f(t, y(t))‖ ≤ L‖yk(t)− y(t)‖, t ∈ Iα
erhalt man f(t, yk(t))→→ f(t, y(t)) auf Iα.
Daher folgt fur jedes feste x ∈ Iα
yk+1(x) = y0 +∫ x
x0
f(t, yk(t)) dt→ y0 +∫ x
x0
f(t, y(t)) dt,
und dahery(x) = y0 +
∫ x
x0
f(t, y(t)) dt, x ∈ Iα;
woraus mit Lemma 5.2 die Behauptung folgt. 2
Wir kommen nun zum Hauptsatz uber die Existenz einer maximalen Losung,die man im Prinzip durch ”Hintereinandersetzen“ von lokalen Losungenerhalt.
Sei y1 ∈ C1[x0, x1] (x1 > x0) die lokale Losung nach Satz 5.3 und y2 ∈C1[x1, x2] die lokale Losung nach Satz 5.3 des AWP’s
y′2(x) = f(x, y2(x))y2(x1) = y1(x1).
Dann ist
y(x) =y1(x), x0 ≤ x ≤ x1
y2(x), x1 ≤ x ≤ x2
x0 x1 x2
y1y2
•
eine lokale Losung des ursprunglichen Problems (die y1 fortsetzt). Denn yist (automatisch) in x1 differenzierbar, da
y′+(x1) = limh→0h>0
y(x1 + h)− y(x1)h
= f(x1, y2(x1) = y′−(x1) = f(x1, y1(x1)).
Ferner ist y′(x) = f(x, y(x)) stetig auf [x0, x2].
5 Gewohnliche und partielle DGL, numerische Methoden 79
Satz 5.4 (Hauptsatz uber die Existenz einer maximalen Losung)Gegeben sei das AWP aus Definition 5.1 mit einem stetigen f , das auf Deine lokale Lipschitzbedingung erfullt. Dann besitzt das AWP (5.15) einemaximale Losung y ∈ C1(I) mit einem offenen Intervall I(§).
In den Beispielen (5.19), (5.20) ist jeweils die lokale Lipschitzbedingungerfullt, so dass die angegebenen Losungen maximale Losungen sind.
BeweisSkizze von Satz 5.4 (nur zur Information): Sei L die Mengealler lokalen Losungen durch (x0, y0). Nach Satz 5.3 ist L 6= ∅. Das Definiti-onsintervall von y ∈ L sei I(y). Stets ist x0 ∈
I (y), y ∈ L.
SeiI =
⋃y∈L
I(y) (was sonst?).
Man zeigt: I ist auch ein Intervall, dass offen sein muss.
Wegen der eindeutigen Losbarkeit der Anfangswertprobleme mit Anfangs-werten in D ist durch die folgende Vorschrift eine Funktion Y (x), x ∈ Idefiniert:
Zu x ∈ I wahle (irgendein) I(y) mit x ∈ I(y), und setzeY (x) = y(x).
Offenbar ist Y Fortsetzung von jedem y ∈ L.
Man zeigt: Y erfullt die Differentialgleichung in jedem x ∈ I.
Y ist daher eine lokale Losung, die jede andere lokale fortsetzt. Y ist daherdie maximale Losung. 2
Einen Sonderfall bilden lineare Systeme der Form
(5.25)y′(x) = A(x)y(x) + b(x)y(x0) = y0.
Hier sind A(x) ∈ Rn×n, b(x) ∈ Rn auf dem offenen Intervall I stetige Funk-tionen; x0 ∈ I. Es ist also D = I×Rn und man sieht leicht, dass die Funktionf(x, y) = A(x)y + b(x) auf D eine lokale Lipschitzbedingung erfullt.
Sei I = (a, b). Nach Satz 5.4 existiert daher die maximale Losung des AWP’s(5.25) y ∈ C1(Imax) mit Imax = (a′, b′) ⊂ I. Wir behaupten Imax = I.
(§)Ware z.B. I = (a, b], so konnte man y uber b hinaus fortsetzen.
5 Gewohnliche und partielle DGL, numerische Methoden 80
Sei etwa b′ < b. Nach dem Randver-halten einer maximalen Losung, Satz5.1, kann nur limx↑b′‖y(x)‖ = ∞ gel-ten. Denn der Rand kann nicht erreichtwerden (x→ b′ < b). a bb′ x
D
← ∂D← ∂D
Rn
Nun ist fur x0 ≤ x < b′
y(x) = y0 +∫ x
x0
A(t)y(t) dt+∫ x
x0
b(t) dt,
bzw.
‖y(x)‖ ≤
(‖y0‖+
∫ b′
x0
‖b(t)‖ dt
)+∫ x
x0
‖A(t)‖ ‖y(t)‖ dt,
woraus mit dem Gronwall-Lemma der Widerspruch limx↑b′‖y(x)‖ <∞ folgt.Also ist b′ = b. 2
Satz 5.5 (Hauptsatz fur lineare Differentialgleichungen) Gegeben seimit einem offenen Intervall I das Anfangswertproblem (5.25) mit stetigenDaten A(x), b(x). Dann existiert die maximale Losung dieses AWP’s auf I.
Losungen zu verschiedenen Anfangswerten des AWP’s (5.25) existieren dem-nach stets auf ganz I, und sie hangen auf jedem kompakten Intervall I1 ⊂I, x0 ∈
I 1, lipschitzstetig von den Anfangswerten ab:
Sind y(y0;x), y2(y∗0;x) zwei Losungen zum Anfangswert y0 bzw. y∗0, so gilt
y(y0;x)− y(y∗0;x) = y0 − y∗0 +∫ x
x0
A(t)(y(y0; t)− y(y∗0; t)) dt, x ∈ I,
bzw. mit dem Gronwall-Lemma
‖y(y0;x)− y(y∗0;x)‖ ≤ eR xx0‖A(t)‖ dt‖y0 − y∗0‖, x ∈ I.
Daher gilt fur ein kompaktes Intervall I1 ⊂ I, x0 ∈I 1
(5.26) ‖y(y0)− y(y∗0)‖∞;I1 ≤ eRI1‖A(t)‖ dt‖y0 − y∗0‖.
Im allgemeinen nichtlinearen Fall kann die Aussage (5.26) so nicht gelten,da maximale Losungen mit ”benachbarten“ Anfangswerten nicht dasselbeExistenzintervall besitzen mussen:
5 Gewohnliche und partielle DGL, numerische Methoden 81
(5.27) Beispiel: Sei D = R× y | y > 0.
y′ = xy(y − 2)y(0) = y0
Fur y0 = 2 lautet die maximale Losung y(2;x) = 2, x ∈ R.Fur y0 > 2 lautet sie
y(x) =2y0
y0 + (2− y0)ex2 mit x ∈ (−xs, xs), xs = log(
y0
y0 − 2
).
Es ist limx→±xs y(y0;x) =∞.Die obige Formel gilt auch fur 0 < y0 < 2; es ist dann limx→±∞ y(y0;x) = 0.
Wegen xs → 0, y0 →∞ ist keine ”globale“ stetige Abhangigkeit zu erwarten.
y
y0
2
y0
−xs → 0 ← xs x
y(2;x)
← y(y0;x)
y(y0;x)
•
•
•
•
Es gilt:
Satz 5.6 Unter den Voraussetzungen des Hauptsatzes (Satz 5.4) sei x0×Br(y0) ⊂ D, BR(y0) = y |‖y − y0‖ ≤ r, und fur ein geeignetes ε > 0 sei(o.E.) auch R = [x0 − ε, x0 + ε]× Br+ε(y0) ⊂ D.
L sei die Lipschitzkonstante von f bzgl. des Rechtecks R. Dann existiertein Intervall Iα = [x0 − α, x0 + α], α > 0, so dass jede maximale Losungy(y∗0;x) mit Anfangswert y∗0 ∈ Br(y0) auch auf Iα definiert ist. Ferner giltfur y∗0, y
∗1 ∈ Br(y0)
(5.28) ‖y(y∗0;x)− y(y∗1;x)‖ ≤ eL|x−x0|‖y∗0 − y∗1‖, x ∈ Iα.
Beweis: Sei M = max(x,y)∈R ‖f(x, y)‖.Fur y∗0 ∈ Br(y0) sei Ry∗0 = (x, y) | |x − x0| ≤ ε, ‖y − y∗0‖ ≤ ε. Dann istRy∗0 ⊂ R ⊂ D, denn fur (x, y) ∈ Ry∗0 ist ‖y−y0‖ ≤ ‖y−y∗0‖+‖y∗0−y0‖ ≤ ε+r.
5 Gewohnliche und partielle DGL, numerische Methoden 82
Also existiert nach dem Beweis von Satz 5.3 (fur das Rechteck Ry∗0 ) einelokale Losung (durch y∗0) auf dem Intervall Iα mit α aus (3) des Beweises:
α = minε,
ε
M,
17L
.
Die Abschatzung (5.28) ergibt sich wieder mit dem Gronwall-Lemma. 2
Die stetige Abhangigkeit von den Anfangswerten ist naturlich eine wichti-ge Voraussetzung fur die numerische Approximierbarkeit der Losung. Dennjedes Verfahren liest den Anfangswert ein.
Wir befassen uns jetzt zunachst mit der numerischen Approximation ge-wohnlicher Anfangswertprobleme und behandeln dann exemplarisch die ”Dis-kretisierung“ einiger partieller Differentialgleichungen (s. Beispiel (5.6)).
Unter den Voraussetzungen des Hauptsatzes wollen wir die maximale Losungdes AWP’s y ∈ C1(I) auf einem beschrankten Teilintervall [a, b] ⊂ I, x0 = a,naherungsweise berechnen. Wir suchen die Naherungswerte in einem hinrei-chend großen Schlauch um die Losung
(5.29) Sδ(y) = (x, z) | ‖y(x)− z‖ ≤ δ, x ∈ [a, b] ⊂ D.
a xb
z
y
← Sδ(y)
Wir konnen daher annehmen, dass f auf Sδ(y) eine Lipschitzbedingungerfullt (Ubung).
Wir beginnen mit dem einfachsten Verfahren, dem Euler-Cauchy Verfahren.Sei t1 = t0 + h0, t0 = a (h0 > 0). Um eine Naherung fur y(t1) zu erhalten,geht man entlang der Tangente der Losung im Punkt (t0, y0):
t0 t0+h0 x
y
y0
y1
y(x)
← y0 + y′(t0)(x− t0) = y0 + f(t0, y0)(x− t0)
•
•
Wir setzen y1 = y0 + h0f(t0, y0). Um eine Naherung fur y(t2) = y(t1 + h1)
5 Gewohnliche und partielle DGL, numerische Methoden 83
zu erhalten, geht man entlang der Tangente der lokalen Losung mit dem(falschen) Anfangswert y1 und erhalt y2 = y1 + h1f(t1, y1); usw.
Auf diese Weise erzeugen wir uns auf dem Gitter Th := tj | ti+1 = ti +hi, i = 0, . . . , N − 1, t0 = 0, tN = b mit dem Schrittweitenvektor h =(h0, . . . , hN−1)t Naherungswerte yi fur y(ti) durch die Rekursion
(5.30) yi+1 = yi + hif(ti, yi), i = 0, . . . , N − 1,
falls dies moglich ist.
Die Naherungslosung ist eine Funktion auf dem Gitter Th, eine Gitterfunk-tion yh(t), t ∈ Th, definiert durch
(5.31) yh(tj) = yj , j = 0, . . . , N − 1.
(Durch Interpolation der yi erhalt man eine kontinuierliche Naherungslosung;
”dense output“)
Uns interessiert naturlich eine Aussage uber den Fehler
(5.32) maxi=0,...,N
‖y(ti)− yi‖ = maxt∈Th‖y(t)− yh(t)‖ = ‖y − yh‖h,∞,
bzw. auch die Frage der Konvergenz : limh→0 ‖y− yh‖h,∞ = 0, d.h. zu jedemε > 0 existiert ein δ(ε) > 0, so dass fur jeden Schrittweitenvektor h mitmaxi=0,...,N−1 hi < δ, ‖y − yh‖h,∞ < ε folgt.
Fur die reine Konvergenzanalyse ist es zweckmaßig, den folgenden lokalenAbbruchfehler Th(t) = T (h, t) einzufuhren. Fur 0 < h < b − a und t ∈[a, b− h] ist der lokale Abbruchfehler des Euler-Cauchy-Verfahrens definiertdurch
(5.33) Th(t) =1h
(y(t+ h)− y(t))− f(t, y(t)),
wobei y die Losung des AWP’s ist.
Oder anschaulich: Startet man das Verfahren im Punkt t ∈ [a, b − h] mitdem exakten Wert y(t), so ist der Fehler nach einem Schritt mit Schrittweiteh
hTh(t) = y(t+ h)− (y(t) + hf(t, y(t))).
Th(t) kann man i.a. nicht ausrechnen, da y unbekannt ist; man kann jedochdie Konvergenzordnung fur h→ 0 bestimmen.
(5.34) Th(t) =O(1), y ∈ C1[a, b]O(h), y ∈ Cm[a, b], m ≥ 2 (Taylor)
5 Gewohnliche und partielle DGL, numerische Methoden 84
Man nennt daher das Verfahren konsistent bzw. konsistent von 1. Ordnung(fur glatte Losungen y ∈ Cm,m ≥ 2). Wir werden bei allen Einschrittver-fahren (allein aus yj ergibt sich yj+1) des Typs (5.30) sehen, dass aus derKonsistenz die Konvergenz folgt (was sonst in der Numerik nicht der Fallist).
Wir fordern die Durchfuhrbarkeitsbedingungen:
(5.35)Fur Schrittweitenvektoren h mit maxhi ≤ H0 existiere die Nahe-rungslosung yh(t), t ∈ Th, und es sei (t, yh(t)) ∈ Sδ(y) ⊂ D miteinem δ > 0.
Bei Anwendung ist (5.35) meistens erfullt; und sie gilt, wie wir sehen werden,wenigstens asymptotisch fur ”H0 klein genug“.
Es ist
y(ti+1) = y(ti) + hif(ti, y(ti)) + hiThi(ti) undyi+1 = yi + hif(ti, yi), i = 0, 1, . . . , N − 1.
Unter der Voraussetzung (5.35) erhalt man mit den Abkurzungen ∆i =‖y(ti)− yi‖, bi = hi‖Thi(ti)‖ unmittelbar die Abschatzung
∆i+1 ≤ (1 + hiL)∆i + bi ≤ ehiL∆i + bi, i = 0, . . . , N − 1,
wobei L die Lipschitzkonstante bzgl. Sδ(y) ist.
Durch Induktion (oder Inspektion) findet man
∆i ≤ eL(h0+...+hi−1)i−1∑k=0
bk, i = 1, 2, . . . , N.
Daher gilt fur i = 0, 1, . . . , N
(5.36) ‖y(ti)− yi‖ ≤ eL(ti−a)i−1∑k=0
hk‖Thk(tk)‖,
woraus naturlich die Konvergenz des Verfahrens folgt.
Satz 5.7 Unter der Voraussetzung (5.35) gilt fur das Euler-Cauchy-Verfah-ren die ”Fehlerabschatzung“
(5.37)‖y(ti)− yi‖ ≤ eL(ti−a)
i−1∑k=0
hk‖Thk(tk)‖
≤ eL(ti−a)(ti − a) maxk=0,1,...,i−1
‖Thk(tk)‖.
Ferner ist die Bedingung (5.35) fur H0 klein genug erfullt.
5 Gewohnliche und partielle DGL, numerische Methoden 85
Beweis: Wir beweisen den Zusatz. Sei ein Sδ(y) ⊂ D gegeben. Wahle dannH0 so klein, dass
eL(b−a)N−1∑k=0
hk‖Thk(tk)‖ < δ
gilt. Durch Induktion zeigt man dann, dass die yi wohldefiniert sind (in Sδliegen), und dass (5.37) gilt. 2
Dieses Verfahren ist bestenfalls nur von erster Ordnung konvergent. Mit(5.37) kann man aber durch Steuerung der Schrittweiten hk den Fehler (theo-retisch) kontrollieren. Wir kommen darauf noch zuruck.
Ein Verfahren 2. Ordnung kann man z.B. so konstruieren:Es ist (siehe (4.46), k = 0)
y(t+ h)− y(t)h
= y′(t+h
2) +O(h2),
falls y ∈ C3.Also ist wegen y′(t+ h
2 ) = f(t+ h2 , y(t+ h
2 ))
Th(t) =y(t+ h)− y(t)
h− f(t+
h
2, y(t+
h
2)) = O(h2).
y(t + h2 ) ist nicht zu gebrauchen (dies ergabe ein 2-Schritt-Verfahren: yj ,
yj+ 12→ yj+1). Wir approximieren bis auf O(h2)-Terme y(t + h
2 ) durchbrauchbare Großen:
y(t+ h2 ) = y(t) + h
2y′(t) +O(h2) = y(t) + h
2f(t, y(t)) +O(h2).
Also ist (Taylor)
f(t+ h2 , y(t+ h
2 )) = f(t+ h2 , y(t) + h
2f(t, y(t))) +O(h2).
Damit haben wir
(5.38) Th :=y(t+ h)− y(t)
h− f(t+ h
2 , y(t) + h2f(t, y(t))) = O(h2).
Das zugehorige Verfahren ist das verbesserte Euler-Cauchy-Verfahren
(5.39)k1 = f(tj , yj)k2 = f(tj + hj
2 , yj + hj2 k1)
yj+1 = yj + hjk2, j = 0, 1, . . . , N − 1.
Man beachte, dass die hohere Konsistenzordnung 2 nur realisiert wird, wennf ∈ C2 ist! Ist f ∈ C1, f /∈ C2, so ist i.a. Th(t) = O(h), und man wurde(5.39) nicht verwenden, da man pro Schritt zwei f -Werte berechnen muss.
5 Gewohnliche und partielle DGL, numerische Methoden 86
Ist das Verfahren (5.39) durchfuhrbar im Sinne der Bedingung (5.35), dannist die Verfahrensfunktion
(5.40) fh(t, y) = f(t+ h2 , y(t) + h
2f(t, y(t)))
fur h ≤ H0 und (t, y) ∈ Sδ, t ∈ [a, b− h2 ] wohldefiniert, da fur (t, y) ∈ Sδ, t ∈
[a, b− h2 ] (t+ h
2 , y + h2f(t, y)) ∈ Sδ′ ⊂ D mit Sδ ⊂ Sδ′ .
Ist daher L die Lipschitzkonstante bzgl. Sδ′ , so gilt fur (t, y), (t, y∗) ∈ Sδ, t ∈[a, b− h
2 ]
‖fh(t, y)− fh(t, y∗)‖ = ‖f(t+ h2 , y + h
2f(t, y))−f(t+ h2 , y∗ + h
2f(t, y∗))‖
≤ L‖(y + h2f(t, y))− (y∗ + h
2f(t, y∗))‖
≤ L‖y − y∗‖+ L2H02 ‖y − y
∗‖
= L(1 + L2H0)‖y − y∗‖
= L0‖y − y∗‖.
Es isty(ti+1) = y(ti) + hifhi(ti, y(ti)) + hiThi(ti)
yi+1 = yi + hifhi(ti, yi), i = 0, . . . , N − 1.
Wir erhalten daher wie beim Euler-Verfahren analog zu (5.37)
(5.41) ‖y(ti)− yi‖ ≤ eL0(ti−a)i−1∑k=0
hk‖Thk(tk)‖, i = 0, . . . , N.
Satz 5.8 Das Verfahren (5.39) sei durchfuhrbar im Sinne von (5.35). Danngilt die Abschatzung (5.41). Ferner ist die Bedingung (5.35) fur H0 kleingenug stets erfullt.
Beweis: Wir beweisen den Zusatz ausfuhrlich (zur Information).Sei Sδ(y) ⊂ D gegeben. M = max(t,y)∈Sδ/2 ‖f(t, y)‖.Es sei L die Lipschitzkonstante von f bzgl. Sδ, und
(1) L0 = L(1 +L
2(b− a)).
Es sei maxi=0,...,N−1 ≤ H0 und H0 ≤ (b− a) so gewahlt, dass
(2) eL0(b−a)N−1∑k=0
hk‖Thk(tk)‖ <δ
2und H0 ≤
δ
2M
5 Gewohnliche und partielle DGL, numerische Methoden 87
gilt. Dann existiert die Naherungslosung yh(t), t ∈ Th, (t, yh(t)) ∈ Sδ/2,t ∈ Th, und es gilt die Abschatzung (5.41).
Wir zeigen dies durch Induktion uber i.Fur i = 0 ist (5.41) richtig.(5.41) sei richtig fur i. Wegen (2) gilt dann ‖y(ti) − yi‖ ≤ δ
2 , d.h. (ti, yi) ∈Sδ/2.
Es ist ‖y(ti + hi2 )− y(ti)‖ ≤
∫ ti+hi/2ti
‖f(t, y(t))‖ dt ≤ H02 M , und daher∥∥∥y(ti + hi
2 )− (yi + hi2 f(ti, yi))
∥∥∥ ≤ H0
2M +
∥∥∥y(ti)− yi − hi2 f(ti, yi)
∥∥∥≤ H0
2M +
δ
2+H0
2M
= H0M +δ
2
≤ δ wegen (2).
Daher ist (ti + hi2 , yi + hi
2 f(ti, yi)) ∈ Sδ, und es existiert
yi+1 = yi + hif(ti + hi2 , yi + hi
2 f(ti, yi)).
y(ti+1) = y(ti) + hif(ti + hi2 , y(ti) + hi
2 f(ti, y(ti) + hiThi(ti))
Es folgt
‖y(ti+1)− yi+1‖ ≤ (1 + hiL0)‖y(ti)− yi‖+ hi‖Thi(ti)‖
≤ ehiL0‖y(ti)− yi‖+ hi‖Thi(ti)‖
(5.41)
≤ ehiL0eL0(ti−a)i−1∑k=0
hk‖Thk(tk)‖+ hi‖Thi(ti)‖
≤ eL0(ti+1−a)i∑
k=0
hk‖Thk(tk)‖ ≤δ
2
2
Wozu braucht man Verfahren hoherer Ordnung?Bei einem Verfahren hoherer Ordnung kann man sicher die Schrittweitengroßer wahlen (¶), um eine Genauigkeit ε zu erzielen. Dafur ist aber auchder Aufwand pro Schritt hoher.
(¶) Wegen der Rundungsfehler kann man die Schrittweiten nicht beliebig klein machen(Ubung)!
5 Gewohnliche und partielle DGL, numerische Methoden 88
Wir machen das folgende Experiment : Gegeben sei eine Menge von Verfahrender Ordnung p = 1, 2, . . . , 50 und ein hinreichend glattes AWP, das dieseOrdnung auch realisiert. Durch Schrittweitensteuerung (s.u.) sei es moglich,jeweils die Genauigkeit ε okonomisch zu erreichen. Fur die Genauigkeit ε =10−r, r = 1, 2, . . . , 14 bestimmen wir den Aufwand des Verfahrens mit derOrdnung p, A(p, r). Fur festes r ist der Aufwand minimal fur das Verfahrenmit der Ordnung pr, A(pr, r) ≤ A(p, r), p = 1, . . . , 50. Man stellt dann fest:
(5.42) pr ∼ c0r + c1 (c0, c1 > 0).
Die folgende heuristische Uberlegung ergibt (5.42): Das Verfahren der Ord-nung p erreicht mit Schrittweitensteuerung
ε '‖y − yh‖h,∞‖y‖h,∞
∼N−1∑i=0
hihpi = (b− a)hp
mit einer mittleren Schrittweite h.
(Es ist namlich min[0,b−a] xp ≤
PN−1i=0 hih
piPN−1
i=0 hi≤ max[0,b−a] x
p, also existiert (Zwi-
schenwertsatz) h ∈ (0, b− a) mit hp =PN−1i=0 hih
pi
b−a ,minhi ≤ h ≤ maxhi).
Also ist ε ∼ κhp. Der Aufwand pro Schritt wird ∼ cp sein, und die Zahlder Schritte ∼ b−a
h . Daher wird fur ε = 10−r gelten: A(p, r) ∼ b−ah p ∼
p(κε )1/p = p(10rκ)1/p. Die Funktion g(p) = p(10rκ)1/p, p > 0 ist konvex. Esist g′(pr) = 0, pr ∼ log(10rκ) = r log 10 + log κ (κ ≥ 1).
Die beiden vorgestellten Einschrittverfahren sind Spezialfalle der allgemei-nen Runge-Kutta-Verfahren, die Runge 1895 und Kutta 1901 begrundethaben.Ausgehend von einem Quadraturverfahren auf [0, 1]
(5.43)∫ 1
0g(t) dt =
m∑i=1
γig(αi) +R0(g), g ∈ C[0, 1]
mit Knoten 0 ≤ α1 ≤ α2 ≤ . . . ≤ αm und Gewichten γi mit
(5.44)m∑i=1
γi = 1 (⇔ R0(1) = 0)
ist wegen y(t+ h) = y(t) +∫ t+ht f(u, y(u)) du
(5.45) y(t+ h) = y(t) + hm∑i=1
γif(t+ αih, y(t+ αih)) + hR0(f(. . .)).
5 Gewohnliche und partielle DGL, numerische Methoden 89
Die y(t+αih) konnen wir zur Konstruktion eines 1-Schritt Verfahrens nichtgebrauchen.
Es seien Quadraturformeln auf [0, αi] (falls αi > 0) gegeben der Form
(5.46)∫ αi
0g(t) dt =
m∑l=1
βilg(αl) +Ri(g), i = 1, 2, . . . ,m.
Wegen y(t+ αih) = y(t) +∫ t+αiht f(u, y(u)) du erhalten wir
(5.47)y(t+ αih) = y(t) + h
m∑l=1
βilf(t+ αlh, y(t+ αlh)) + hRi(f(. . .)),
i = 1, 2, . . . ,m.
Das numerische Verfahren entsteht naturlich durch Streichen der RestgliederR0, R1, . . . , Rm in (5.45), (5.47).Man erhalt dann mit t = ti, h = hi, y(t)→ yj , y(t+αih)→ zi die Vorschrift
(5.45’) yj+1 = yj + hj
m∑i=1
γif(t+ αihj , zi),
wobei die zi (die) Losung von
(5.47’) zi = yj + hj
m∑l=1
βilf(t+ αlhj , zl), i = 1, 2, . . . ,m
sein sollen.Man braucht in (5.45’) nur die K-Werte
(5.48) ki = f(t+ αihj , zi),
und wegen (5.47’) sind diese (die) Losung von
(5.49) ki = f(tj + αihj , yj + hj
m∑l=1
βilkl), i = 1, 2, . . . ,m.
Fur die Durchfuhrbarkeit des Verfahrens werden wir naturlich voraussetzen,dass fur aktuelle Schrittweiten hi ≤ H0 (5.49) eindeutig losbar ist. Diesgilt zumindest, wenn H0 ”klein genug“ ist: Fur hj = 0 ist (5.49) losbar;kj(0) = f(tj , yj). Sei Fi(k1, . . . , km, h) = ki − f(tj + αih, yj + h
∑ml=1 βilkl)
und F = (F1, . . . , Fm)t, so ist F (k1(0), . . . , km(0), 0) = 0 und falls f ∈C1, ∂F
∂k (k1(0), . . . , km(0), 0) = 1. Nach dem Satz uber implizite Funktionenexistiert daher fur hj klein genug genau eine Losung von (5.49).
5 Gewohnliche und partielle DGL, numerische Methoden 90
Definition 5.4 Sei m ∈ N. Ein m-stufiges Runge-Kutta-Verfahren ist de-finiert durch
yj+1 = yj + hj
m∑i=1
γiki, j = 0, 1, . . . , N − 1
mit den K-Werten ki = ki(hj , tj , yj) als Losung von (5.49), wobei die αi ∈[0, 1], γi mit (5.44) und βil ∈ R die das Verfahren definierenden Parametersind.
Schema:
α1 β11 · · · β1m...
......
αm βm1 · · · βmmγ1 · · · γm
Das Verfahren heißt explizit, falls βil = 0 fur l ≥ i. Dann ist namlichk1 = f(tj + α1hj , yj) und ki = f(tj + αihj , yj + hj
∑i−1l=1 βilkl).
Das Verfahren heißt halbimplizit, falls βil = 0 fur l > i. Dann hat man furjedes ki nur eine implizite Gleichung. Sonst heißt das Verfahren vollimplizit.
Ein vollimplizites Verfahren ist naturlich sehr aufwendig, und man verwen-det dieses auch nur in besonderen Situationen (s.u.).
Den lokalen Abbruchfehler Th(t), t ∈ [a, b − h], h ≤ H0 erhalt man wieder,indem man ausgehend von der exakten Losung im Punkt (t, y(t)) einenSchritt mit der Schrittweite h rechnet und die Differenz zum exakten Wertbildet:
(5.50) y(t+ h) = y(t) + h
m∑i=1
γiki(h, t, y(t)) + hTh(t)
Mit Hilfe des Banachschen Fixpunktsatzes zeigt man im vollimpliziten Fall,dass fur (t, y) ∈ Sδ, t ∈ [a, b− h], h ≤ H0 klein genug
(5.49’) ki = f(t+ αih, y + h
m∑l=1
βilkl), i = 1, 2, . . . ,m
genau eine Losung ki(h, t, y) in Br(k0), k0 = (f(t, y), . . . , f(t, y))t besitzt.Insbesondere gilt dann wegen der gleichmaßigen Stetigkeit von f in S′δ
(5.51) max(t,y)∈Sδt∈[a,b−h]
‖ki(h, t, y)− f(t, y))‖ → 0, h→ 0.
Wegen∑m
i=1 γi = 1 erhalten wir auch
(5.52) maxt∈[a,b−h]
‖Th(t)‖ → 0, h→ 0,
5 Gewohnliche und partielle DGL, numerische Methoden 91
so dass m-stufige Runge-Kutta-Verfahren konsistent sind.
Wegen (5.49’) existiert die Verfahrensfunktion
(5.53) fh(t, y) =m∑i=1
γiki(h, t, y) fur t ∈ [a, b− h],
h klein genug, und (t, y) ∈ Sδ.Sie genugt auch dort einer Lipschitzbedingung. Denn fur (t, y), (t, y∗) ∈ Sδist
‖ki(h, t, y)− ki(h, t, y∗)‖
= ‖f(t+ αih, y + hm∑l=1
βilkl(h, t, y)− f(t+ αih, y∗ + h
m∑l=1
βilkl(h, t, y∗)‖
≤ L(Sδ′)‖y − y∗‖+ h maxi=1,...,m
m∑l=1
|βil| maxl=1,...,n
‖kl(h, t, y)− kl(h, t, y∗)‖
= L(Sδ′)‖y − y∗‖+ hC maxl=1,...,m
‖kl(h, t, y)− kl(h, t, y∗)‖;
woraus fur h klein genug , z.B. hC ≤ 12 ,
(5.54) maxi=1,...,m
‖ki(h, t, y)− ki(h, t, y∗)‖ ≤ 2L(Sδ′)‖y − y∗‖
folgt.
Damit zeigt man dann wie gehabt:
Satz 5.9 Das m-stufige Runge-Kutta-Verfahren sei durchfuhrbar im Sinnevon (5.35), wobei im vollimpliziten Fall
(5.55)
∥∥∥∥∥m∑i=1
γi(ki(h, tj , yj)− ki(h, tj , y(tj)))
∥∥∥∥∥ ≤ L0‖yj − y(tj)‖
noch gelten soll.
Dann gilt die Abschatzung (5.41). Fur H0 klein genug sind die Bedingungen(5.35) und (5.55) erfullt.
5 Gewohnliche und partielle DGL, numerische Methoden 92
Wir geben noch weitere Beispiele von Runge-Kutta-Verfahren an:
Klassisches Runge-Kutta-Verfahren (m = 4)(5.56)
0 k1 = f(tj , yj)12
12 k2 = f(tj + 1/2hj , yj + 1/2hjk1)
12 0 1
2 k3 = f(tj + 1/2hj , yj + 1/2hjk2)
1 0 0 1 0 k4 = f(tj + hj , yj + hjk3)16
13
13
16
yj+1 = yj + hj16
(k1 + 2(k2 + k3) + k4).
Es ist nicht ganz einfach, Th(t) = O(h4) fur f ∈ C4 zu zeigen.
Trapezregel (m = 2)(5.57)
0 k1 = f(tj , yj)
1 12
12 k2 = f(tj+1, yj + 1
2hj(k1 + k2))
12
12
yj+1 = yj +12hj(f(tj , yj) + f(tj+1, yj+1)).
Es ist Th(t) = O(h2), f ∈ C2.
ruckwartiger Euler (m = 1)(5.58)
1 1 k1 = f(tj + hj , yj + hjk1)1
yj+1 = yj + hjk1(⇒ k1 = f(tj+1, yj+1)), oderyj+1 = yj + hjf(tj+1, yj+1).
Es ist Th(t) = O(h), f ∈ C1.
5 Gewohnliche und partielle DGL, numerische Methoden 93
Runge-Kutta-Gauß (m = 2)(5.59)
16(3−
√3) 1
4112(3− 2
√3)
16(3 +
√3) 1
12(3 + 2√
3) 14
12
12
k1 = f(tj +16
(3−√
3)hj , yj +112hj(3k1 + (3− 2
√3)k2))
k2 = f(tj +16
(3 +√
3)hj , yj +112hj((3 + 2
√3)k1 + 3k2))
yj+1 = yj +12hj(k1 + k2).
Es ist Th(t) = O(h4), f ∈ C4.
Auf die Verwendung impliziter Verfahren gehen wir spater ein.
Wir kommen jetzt zur Schrittweitensteuerung . Um diese richtig zu begrun-den, benutzen wir eine etwas andere Fehlerabschatzung als in (5.41). DerAbbruchfehler Th(t) ist eigentlich nicht sauber berechenbar. Wir setzen vor-aus, dass das Verfahren durchfuhrbar ist im Sinne von (5.35). Das Ziel ist,die Schrittweiten hi,maxhi ≤ H0, moglichst okonomisch zu wahlen, um eineFehlertoleranz ε in der Naherungslosung zu garantieren.
Dabei setzen wir voraus, dass bei gegebener Naherungslosung yj , j = 0, 1,. . . , N (mit maxhi ≤ H0) die Anfangswertprobleme
(5.60)y′(x) = f(x, y(x))y(ti) = yi fur x ∈ [ti, ti+1]
auf [ti, ti+1] (eindeutig) losbar sind mit Losung y(yi;x), x ∈ [ti, ti+1], i =0, 1, . . . , N − 1.(Fur hinreichend kleines H0 ist (5.60) erfullt, s. Beweis von Satz 5.5)
Es ist dann
y(ti+1)− yi+1 = y(ti+1)− y(yi; ti+1) + y(yi; ti+1)− yi+1.
Fur 0 < h ≤ hi ist der die Naherungslosung begleitende Abbruchfehler Thdefiniert durch
(5.61) hTi(h) = y(yi; ti + h)− yi − hfh(ti, yi), i = 0, 1, . . . , N − 1.
Mit y(x) = y(yi;x) ist offenbar
Ti(h) =1h
(y(ti + h)− y(ti)− hfh(ti, y(ti))),
5 Gewohnliche und partielle DGL, numerische Methoden 94
so dass sich Ti(h) fur h → 0 genauso abschatzen laßt wie Th(ti) (unter ent-sprechenden Glattheitsvoraussetzungen). Also erhalten wir mit Ti := Ti(hi)
(5.62) y(ti+1)− yi+1 = y(ti+1)− y(yi; ti+1) + hiTi.
Es ist fur ti ≤ t ≤ ti+1
y(t)− y(yi; t) = y(ti)− yi +∫ t
ti
(f(y, y(u))− f(u, y(yi;u))) du,
und mit der Lipschitzkonstanten L = L(Sδ) erhalt man
‖y(t)− y(yi; t)‖ ≤ ‖y(ti)− yi‖+ L
∫ t
ti
‖y(u)− y(yi;u)‖ du.
Das Gronwall-Lemma ergibt daher die Abschatzung
(5.63) ‖y(ti+1)− y(yi; ti+1)‖ ≤ ‖y(ti)− yi‖eLhi .
Setzt man dies in (5.62) ein, so erhalt man die Rekursion
‖y(ti+1)− yi+1‖ ≤ eLhi‖y(ti)− yi‖+ hi‖Ti‖, i = 0, 1, . . . , N − 1,
woraus schließlich die Abschatzung
(5.64) ‖y(ti)− yi‖ ≤ eL(ti−a)i−1∑k=0
hk‖Tk‖, i = 0, . . . , N,
folgt.
Satz 5.10 Das m-stufige Runge-Kutta-Verfahren sei durchfuhrbar im Sin-ne von (5.35), und es seien die Bedingungen (5.60) erfullt. Dann gilt dieAbschatzung (5.64) mit der Lipschitzkonstanten L = L(Sδ).Fur H0 klein genug gilt (5.64).
Die Schrittweitensteuerung versucht die hk so zu wahlen, dass hk‖Tk‖ ' tolgilt. Wegen (5.64) wird dann der Gesamtfehler auch diese Großenordnungbesitzen.Mit einem Verfahren der Ordnung p haben wir bis zum Punkt ti gerechnetund suchen eine geeignete Schrittweite h∗i .
ti ti+h ti+2h
yi
y(yi;x)
•
•
5 Gewohnliche und partielle DGL, numerische Methoden 95
Mit einer Testschrittweite h rechnen wir zwei Schritte zum Punkt ti + 2h,d.h. wir berechnen yh(ti + h) = yi + hfh(ti, yi) und yh(ti + 2h) = yh(ti +h) + . . .. Dann machen wir einen Kontrollschritt mit der Schrittweite 2h,d.h. y2h(ti + 2h) = yi + 2hf2h(ti, yi).
Fur Verfahren auf einem aquidistanten Gitter kann man die folgende asym-ptotische Entwicklung zeigen:Es gibt Funktionen ei(x), gi(x) ∈ C1[ti, b] mit ei(ti) = gi(ti) = 0, so dass furdie auf dem Gitter x = ti, ti + h, ti + 2h, . . . ≤ b mit dem Anfangswert yiberechnete Naherungslosung yh(x) gilt (s. Literatur):
(5.65) yh(x)− y(yi;x) = ei(x)hp + gi(x)hp+1 +O(hp+2).
Es ist yh(ti + h)− y(yi, ti + h) = −hTi(h) und daher
(5.66) hTi(h) = −e′i(ti)hp+1 +O(hp+2).
Wegen (5.65) haben wir
yh(ti + 2h)− y(yi; ti + 2h) = ei(ti + 2h)hp + gi(ti + 2h)hp +O(hp+2)= e′i(ti)2hh
p +O(hp+2);
yh(ti + 2h)− y(yi; ti + 2h) = e′i(ti)2hp+1 +O(hp+2),
y2h(ti + 2h)− y(yi; ti + 2h) = e′i(ti)(2h)p+1 +O(hp+2),
woraus
1(2p − 1)
(y2h(ti + 2h)− yh(ti + 2h))(5.67)
= 2e′i(ti)hp+1 +O(hp+2)
(5.66)= −2hTi(h) +O(hp+2)
folgt.
Wir konnen daher hTi bis auf Fehler hoherer Ordnung durch den berechen-baren Ausdruck in (5.67) darstellen. Wir setzen als lokalen Schatzer von2h‖Ti(h)‖
(5.68) err(h) =1
2p − 1‖y2h(ti + 2h)− yh(ti + 2h)‖.
Die gewunschte Schrittweite h∗ erfullt
(5.69) 2h∗‖Ti(h∗)‖ ' tol ∼ ‖l′(ti)‖2(h∗)p+1.
Wegen (s. (5.67)) err(h) ∼ ‖l′(ti)‖2hp+1 erhalten wir
(5.70)h∗
h'(
tolerr(h)
) 1p+1
.
5 Gewohnliche und partielle DGL, numerische Methoden 96
Bemerkung: Wurde man in (5.69) ‖Ti(h∗)‖ ' tol fordern, so erhielte man
statt (5.70) h∗
h '(
htolerr(h)
) 1p ; was naturlich auch geht.
Man kann die folgende Strategie benutzen:
(1) Ist err(h) ≤ tol, dann werden beide Werte yh(ti + h), yh(ti + 2h) ak-zeptiert, und man startet im Punkt ti + 2h mit der Schrittweite hneu
nach der Formel (5.71).
(2) Ist err(h) > tol, dann startet man die Rechnung im Punkt ti erneutmit der (kleineren) Schrittweite hneu.
(5.71) hneu = hmin
qmax,max
qmin, s
(tol
err(h)
) 1p+1
.
Hier ist qmax ∈ [1, 5] (z.B), qmin < qmax und der tuning-Faktor s ∼ 0.8, 0.9,der eine zu schnelle Anderung der Schrittweite vermeiden kann.
Insbesondere ist fur err(h) tol
hneu = hmin
qmax, s
(tol
err(h)
) 1p+1
∼ h qmax
und fur err(h) tol hneu = hqmin.
Wir befassen uns noch kurz mit der Verwendung impliziter Verfahren. Sei
(5.72) y′(x) =(
998 1998−999 −1999
)y(x), x ≥ 0, y(0) =
(10
).
Die Matrix A in (5.72) besitzt die Eigenwerte λ1 = −1, λ2 = −1000 mit den
(l.u.) Eigenvektoren c1 =(
2−1
), c2 =
(−11
).
Offenbar sind c1ex, c2e
−1000x und jede Linearkombination auch Losung vony′ = Ay. Wegen y(0) = α1c1 +α2c2 mit eindeutigen α1, α2 lautet die Losungdes AWP’s (5.72)
(5.73) y(x) =(
2−1
)e−x +
(−11
)e−1000x =
(y1(x)y2(x)
).
Da sich die Losung aus sehr unterschiedlich schnell fallenden Anteilen zu-sammensetzt, nennt man das System steif .
5 Gewohnliche und partielle DGL, numerische Methoden 97
Ebenso fuhrt die Linienmethode fur die Warmeleitungsgleichung (5.6) aufein System (s. (5.11))
y′ = Ay + b mit A = − 1h2
2 −1−1 2 −1
. . .−1 2 −1
−1 2
,
das ebenfalls steif ist. Denn die Eigenwerte von A lauten λj = − 4h2 sin2 jπ
2N ,
j = 1, . . . , N − 1, und es ist max |λj |min |λj | = |λN−1|
|λ1| '4π2N
2 1, wenn h kleinwird.
Lost man jetzt (5.72) mit einem expliziten Verfahren, z.B. Euler, so erhaltman mit der Transformation
yj = Tzj , T = (c1, c2), zj+1 = (1 + h diag(λ1, λ2))zj ,
oder
(5.74)y1j = 2(1− h)j − (1− 1000h)j
y2j = −(1− h)j + (1− 1000h)j , j = 0, 1, 2, . . . .
Man sieht, dass Konvergenz auf einem festen Intervall nicht alles ist: Bei(5.74) muss wahrend der gesamten Rechnung |1−1000h| < 1 sein, damit dienumerische Losung nicht exponentiell explodiert! Es muss also h < 0.002sein, obwohl der Term, der das bewirkt, fur große j keinen Beitrag zurLosung bringt.
Mit expliziten Verfahren kann man daher steife AWP’s auf großen Inter-vallen nicht losen. Offenbar liegt dies daran, dass ejhλ = (ehλ)j durch ein(Polynom(λh))j approximiert wird, und |Polynom(λh)| < 1 immer eine Be-schrankung fur h bedeutet.
Bei geeigneten, teuren, impliziten Verfahren ist das nicht der Fall. Dortklingt unabhangig von h auch die numerische Losung (fur j → ∞) ab, sodass die Schrittweite nur anfangs aus Genauigkeitsgrunden kleiner zu wahlenist.
Das ruckwartige Euler-Verfahren (5.58) ergibt z.B.
(5.75)y1j = 2 1
(1+h)j− 1
(1+1000h)j
y2j = − 1
(1+h)j+ 1
(1+1000h)j, j = 0, 1, . . .
5 Gewohnliche und partielle DGL, numerische Methoden 98
Hier wird (ehλ)j approximiert durch (rat. Funktion(hλ))j , und rationaleFunktion konnen fur die obigen hλ beschrankt sein.
Man erhalt eine vollstandige Diskretisierung des sog. parabolischen Anfangs-randwertproblems (5.6), indem man auf die Semidiskretisierung (5.11) eingeeignetes implizites Verfahren anwendet. Wir nehmen ein Verfahren 2. Ord-nung, die Trapezregel (5.57) und erhalten das Crank-Nicolson Verfahren.Wir formulieren dies nur auf einem aquidistanten Gitter. Sei k die Schritt-weite in t Richtung.
(5.76) k =T
M, tj = jk, j = 0, . . . ,M.
Mit der Matrix A aus (5.11) und mit
(5.77) b(t) =1h2
ϕ1(t)
0...0
ϕ2(t)
erhalt man fur j = 0, 1, . . . ,M − 1
yj+1 = yj +k
2(Ayj + b(tj) +Ayj+1 + b(tj+1))
bzw. das Gleichungssystem
(5.78)(1− k
2A)yj+1 = (1 + k
2A)yj + k2 (b(tj) + b(tj+1)),
j = 0, . . . ,M − 1,y0 = y(0).
Man kann zeigen, dass hier unter entsprechenden Glattheitsvoraussetzungenmit u(tj) = (u(h, jk), . . . , u((N − 1)h, jk))t die Abschatzung
(5.79) maxj=0,...,M
‖yj − u(tj)‖ = C(u)(h2 + k2)
gilt(‖); und zwar ohne Bedingung an die Schrittweiten (z.B. λ = kh2 ≤ 1),
die man sonst bei parabolischen Problemen hat. Das System (5.78) lost manfur jedes j = 0, 1, . . . ,M − 1 naturlich iterativ.
(‖)M. Lees: Approximate Solutions of Parabolic Equations, I. Soc. Indust. Appl. Math.7 (1959), 167-183.
6 Optimierung 99
6 Optimierung
Die Optimierung befasst sich mit der Minimierung gewisser (Kosten-) Funk-tionale unter Nebenbedingungen. Naturgemaß gibt es daher sehr viele An-wendungen (Ubung).
Definition 6.1 Sei f : M ⊂ Rn → R. x ∈ M heißt ein lokales Minimumvon f auf M , wenn eine Kugel Br(x) existiert mit
f(x) ≤ f(x), x ∈ Br(x) ∩M ;
und x heißt ein globales Minimum, wenn
f(x) ≤ f(x), x ∈M ;
gilt. Entsprechend lautet die Definition fur ein striktes (lokales) Minimum.
Aus der Analysis sind hinreichende Bedingungen bekannt, wann ein Mini-mum existiert. Es ist der Traum jedes Optimierers, einen Algorithmus zurBerechnung eines globalen Minimums zu entwerfen!
Die Grundlage fur Optimalitatskriterien ist
Satz 6.1 Sei M ⊂ D ⊂ Rn, D offen und f : D → Rn sei stetig differen-zierbar; f ∈ C1(D).
(i) Ist x ∈M ein lokales Minimum von f auf M , so gilt fur jede Richtungn (‖n‖ = 1) mit x+ tn ∈M , 0 ≤ t ≤ δ (δ > 0)
(6.1) Dnf(x) = (gradf(x), n) ≥ 0.
(ii) Ist x ∈M ein lokales Minimum von f auf M , so ist
(6.2) f ′(x) = 0 (kritischer Punkt).
Ferner ist fur f ∈ C2(D)
(6.3) f ′′(x) positiv semidefinit (p.s.d.),
d.h. (x, f ′′(x)x) ≥ 0, x ∈ Rn.
6 Optimierung 100
Beweis:(i) Es ist F (t) := f(x + tn) ∈ C1[0, δ] und F (t)−F (0)
t ≥ 0, t ∈ (0, δ]. Alsoist F ′+(0) ≥ 0; aber F ′+(t) = (gradf(x+ tn), n).
(ii) Ist x ∈M , so gilt (6.1) fur jede Richtung n. Daher ist x ein kritischer
Punkt.Sei n beliebig fest. Dann besitzt die Funktion F (t) = f(x + tn), t ∈(−δ, δ) mit einem geeigneten δ > 0 in t = 0 ein Minimum. Also istF ′′(0) = (n, f ′′(x)n) ≥ 0. 2
Die Bedingung (6.3) (mit (6.2)) ist nicht hinreichend fur das Vorliegen einesMinimums (z.B. f(x) = x3
1, x = 0). Es konnen z.B. Sattelpunkte(∗) (d.h.f ′′(x) ist indefinit) auftreten.
Satz 6.2 Unter den Voraussetzungen von Satz 6.1 sei f ∈ C2(D) und x ∈M
ein kritischer Punkt. Existiert eine Kugel Br(x) ⊂M , so dass
(6.3’) f ′′(x) p.s.d. fur x ∈ Br(x) ist,
dann ist x ein lokaler Minimalpunkt.
Beweis: Fur x ∈ Br(x) ist nach der Taylorformel
f(x) = f(x) + f ′(x)(x− x) +12
(x− x, f ′′(θx+ (1− θ)x)(x− x))
mit einem 0 < θ < 1 und θx+ (1− θ)x ∈ Br(x).Also folgt f(x) ≥ f(x), x ∈ Br(x). 2
Bemerkung: Hinreichend fur (6.3’) ist naturlich die Bedingung: f ′′(x) po-sitiv definit. Dann ist x auch ein strikt lokales Minimum (Ubung).
(6.4) Beispiel:
(i) Die Rosenbrock -Funktion
f(x1, x2) = 100(x2 − x21)2 + (1− x1)2
besitzt den einzigen kritischen Punkt x =(
11
)und f ′′(x) ist s.p.d.
(ii) Die Funktion
f(x1, x2) = (x21 + x2 − 11)2 + (x1 + x2
2 − 7)2
besitzt vier globale Minimalstellen mit Wert Null und vier Sattelpunk-te.
(∗) z.B. f(x1, x2) = x1x2, x = 0, oder f(x1, x2) = x21 − x2
2, x = 0.
6 Optimierung 101
(iii) Die Funktion
f(x1, x2) = (x1 − 2)4 + (x1 − 2x2)2
besitzt ein globales Minimum bei x =(
21
), f ′′(x) ist singular und
f ′′(x), x 6= x, ist s.p.d..
Es gibt eine große und wichtige Klasse von Problemen, fur die die notwendi-gen Optimalitatsbedingungen oft auch hinreichend sind: konvexe Probleme.Hier ist die Menge M konvex, d.h. bekanntlich mit x, y ∈ M ist auch dieVerbindungsstrecke λx + (1 − λ)y, λ ∈ [0, 1], in M . Und das Zielfunktionalf ist konvex.
Definition 6.2 Sei K ⊂ Rn konvex, und f : K → R. Die Funktion f heißtkonvex, wenn fur x, y ∈ K und 0 ≤ λ ≤ 1 gilt
(6.5) f(λx+ (1− λ)y) ≤ λf(x) + (1− λ)f(y).
x y
f
•
•
λf(x)+(1−λ)f(y)
f(x)+f ′(x)(y−x)
Satz 6.3 Sei f : K ⊂ D ⊂ Rn → R, K konvex, D offen und f ∈ C1(D).Dann ist f konvex genau dann, wenn
(6.6) f(y) ≥ f(x) + f ′(x)(y − x), x, y ∈ K
gilt.
Beweis: ⇒: Sei f konvex. Es ist
f ′(x)(y − x) =d
dtf(x+ t(y − x))
∣∣t=0
= limt↓0
1t(f(x+ t(y − x))− f(x)).
f(x+ t(y − x))− f(x) = f((1− t)x+ ty)− f(x)≤ (1− t)f(x) + tf(y)− f(x)= t(f(y)− f(x)) fur 0 ≤ t ≤ 1.
6 Optimierung 102
Daher gilt (6.6).
⇐: Sei z ∈ K. Dann gilt wegen (6.6) fur x, y ∈ K fur 0 ≤ λ ≤ 1
λf(y) ≥ λf(z) + λf ′(z)(y − z)(1− λ)f(x) ≥ (1− λ)f(z) + (1− λ)f ′(z)(x− z),
woraus durch Addition
λf(y) + (1− λ)f(x) ≥ f(z) + f ′(z)(λy + (1− λ)x− z)
folgt.
Wahlt man hier z = λy + (1− λ)x, so ergibt sich die Konvexitat von f . 2
Nutzlich ist haufig die Beschreibung der Konvexitat uber die zweite Ablei-tung.
Satz 6.4 Sei K ⊂ D ⊂ Rn, K konvex, D offen und f : D → Rn sei ausC2(D).
(i) Ist f ′′(x), x ∈ K, p.s.d., dann ist f konvex auf K.
(ii) Ist f konvex auf K, so ist f ′′(x) fur jedes x ∈K p.s.d..
Beweis:
(i) Es gilt nach der Taylorformel fur x, y ∈ K
(1) f(y) = f(x) + f ′(x)(y − x) +12
(y − x, f ′′(θx+ (1− θ)y)(y − x))
mit einem 0 < θ < 1. Daher gilt f(y) ≥ f(x) + f ′(x)(y − x), was mitSatz 6.3 die Behauptung ist.
(ii) Ist x ∈K und z ∈ Rn beliebig, so gilt y = x+ tz ∈ K fur |t| < δ. Mit
Satz 6.3 folgt daher aus (1) fur t 6= 0 (z, f ′′(θx+ (1− θ)y)z) ≥ 0, bzw.fur t→ 0 (z, f ′′(x)z) ≥ 0. 2
(6.7) Beispiel: Fur A ∈ Rm×n, b ∈ Rn ist f(x) = 12‖Ax − b‖
22 konvex auf
Rn. Denn f ′′(x) = AtA ist p.s.d..
Der entscheidende Satz fur konvexe Probleme ist
Satz 6.5 Sei K ⊂ D ⊂ Rn, K konvex, D offen und f ∈ C1(D), f : D → R,konvex auf K. Dann sind folgende Aussagen aquivalent:
6 Optimierung 103
(i) x ∈ K ist ein globaler Minimalpunkt von f auf K.
(ii) x ∈ K ist ein lokales Minimum von f auf K.
(iii) f ′(x)(x− x) ≥ 0, x ∈ K.
Beweis: (i)⇒(ii): klar
(ii)⇒(iii): Sei x ∈ K und o.E. x 6= x, n = 1‖x−x‖(x− x).
x+ tn = (1− t
‖x− x‖)x+
t
‖x− x‖x ∈ K fur 0 ≤ t ≤ ‖x− x‖.
Daher ergibt (6.1) die Aussage (iii).
(iii)⇒(i): Da f konvex ist, folgt aus (6.6) fur x ∈ K
f(x) ≥ f(x) + f ′(x)(x− x) ≥ f(x).
2
Insbesondere gilt:
Satz 6.6 Unter den Voraussetzungen von Satz 6.5 sind folgende Aussagenaquivalent:
(i) x ∈K ist globaler Minimalpunkt,
(ii) x ∈K ist lokaler Minimalpunkt,
(iii) f ′(x) = 0, d.h. x ist kritischer Punkt.
Im Beispiel (6.7) ist K = Rn und x daher globaler Minimalpunkt genaudann, wenn f ′(x) = 0, d.h. AtAx = Atb (s. Kapitel 2!).
Nichtlineare Optimierungsprobleme kann man i.a. nicht geschlossen losen.Ein numerisches Verfahren wird iterativ eine Folge xk∞k=0 zulassiger xk ∈M,k = 0, 1, . . . erzeugen, die hoffentlich gegen ein lokales Minimum x von fkonvergiert.Bei den sog. Abstiegsverfahren wahlt man, um Maxima zu vermeiden, diexk so, dass
f(xk+1) < f(xk), k = 0, 1, . . .
gilt. Oft ist man mit einem suboptimalen xk zufrieden.
Man unterscheidet zwei Typen von Verfahren:Bei den Verfahren mit Schrittweitensteuerung geht man von einer Abstiegs-richtung dk (s.u.) aus;
f(xk + sdk) < f(xk), 0 < s ≤ sk,
6 Optimierung 104
und wahlt dann die Schrittweite σk ∈ (0, sk] so, dass der Abstieg moglichstgroß ausfallt, und setzt
xk+1 = xk + σkdk.
Bei den Trust-Region-Verfahren (die hochaktuell sind) wird f(xk+d) durcheine lokales Modell fk(d) auf der Trust-Region ‖d‖ ≤ ρk approximiert. Eswird eine Abstiegsrichtung dk des Modells ermittelt, so dass der Abstieg vonf(xk + d) nahezu mit dem Abstieg von fk(d) ubereinstimmt, d.h. f(xk) −f(xk + dk) ' f(xk)− fk(dk); und man setzt xk+1 = xk + dk. Das Levenberg-Marquardt-Verfahren ist ein Verfahren dieser Art.
Im folgenden befassen wir uns mit dem unrestringierten Optimierungspro-blem
(6.8) minx∈Rn
f(x), f : Rn → R,
und wir suchen ein lokales Minimum x von f .Dazu ist naturlich
(6.9) f ′(x) = 0
notwendig.
Man konnte z.B. das Newton-Verfahren auf (6.9) anwenden. Wir suchenaber unter den Nullstellen von f ′ ein lokales Minimum von f . Daher istes sinnvoller, eine Abstiegsvariante zu verwenden, namlich das gedampfteNewton-Verfahren; dann wird die Konvergenz gegen ein lokales Maximumvermieden (Sattelpunkte konnen aber auftauchen).
Die Grundlage fur Abstiegsverfahren ist das einfache
Lemma 6.1 Sei f : Rn → R differenzierbar in x. Sei d ∈ Rn mit
(6.10) (gradf(x), d) < 0.
Dann existiert zu jedem 0 < α < 1 ein s0(x, d, α) > 0, so dass fur 0 ≤ s ≤ s0
gilt
(6.11) f(x+ sd)− f(x) ≤ sα(gradf(x), d).
Beweis: ddsf(x+ sd)
∣∣t=0
= (gradf(x), d) < 0.
Also existiert fur 0 < β < 1 zu −β(gradf(x), d) ein s0 > 0 mit
f(x+ sd)− f(x)s
≤ (1− β)(gradf(x), d) fur 0 < |s| ≤ s0.
Mit α := 1− β folgt daher (6.11). 2
6 Optimierung 105
Definition 6.3 Ein Vektor d mit (6.10) heißt Abstiegsrichtung von f in x.
Ist A ∈ Rn×n eine s.p.d. Matrix und gradf(x) 6= 0, so ist
(6.12) d = −A gradf(x)
eine Abstiegsrichtung. Denn (Ax, x) ≥ c‖x‖2(c > 0), und daher
(gradf(x), A gradf(x)) ≥ c‖gradf(x)‖2 bzw.
(gradf(x),−A gradf(x)) ≤ −c‖gradf(x)‖2 < 0.
Das Newton-Verfahren angewendet auf F (x) = gradf(x) = 0 lautet(F ′(x) = f ′′(x))
(6.13)f ′′(xk)∆xk = −gradf(xk)
xk+1 = xk + ∆xk = xk + dkN
mit der Newtonrichtung
(6.14) dkN = −f ′′(xk)−1gradf(xk).
Ist daher f ′′(xk) positiv definit, so ist dkN eine Abstiegsrichtung, aber esmuss nicht f(xk+1) < f(xk) sein. Dies gelingt erst durch ”Dampfung“ miteiner richtigen Schrittweite 0 < σk ≤ 1, die wir anschließend begrunden.Dies ergibt ein gedampftes Newton-Verfahren.
Sei (gradf(x0), d) < 0. Das Ziel ist, in (6.11) ein moglichst großes Inter-vall [0, s0] zu erreichen. Dazu machen wir die (bei Optimierern ubliche undBeweise erleichternde) Voraussetzung
(6.15) N(x0) = x ∈ Rn | f(x) ≤ f(x0) sei kompakt. (!)
Fur festes 0 < α < 1 betrachten wir fur f ∈ C1 die Funktion
(6.16) g(t) = f(x0 + td)− f(x0)− tα(gradf(x0), d), t ≥ 0.
Wegen Lemma 6.1 ist die Menge
(6.17) A = T > 0 | g(t) < 0, 0 < t < T nicht leer.
Offenbar ist fur T ∈ A x0 + td ∈ N(x0) bzw. x0 + Td ∈ N(x0). Also ist Abeschrankt.Sei
(6.18) σ∗ = supA.
6 Optimierung 106
Es ist dann g(t) < 0, 0 < t < σ∗ und
(6.19) g(σ∗) = 0
(denn es ist g(σ∗) ≤ 0, und g(σ∗) < 0 ergabe einen Widerspruch).
Neben (6.15) setzten wir noch (fast o.E.)
(6.20) ‖f ′(x)− f ′(y)‖ ≤ L‖x− y‖, x, y ∈ N(x0)
voraus.
Fur 0 ≤ t ≤ σ∗ ist
g(t) = g(0) + g′(0)t+∫ t
0(g′(u)− g′(0)) du.
g(0) = 0, g′(t) = (gradf(x0 + td), d)− α(gradf(x0), d).
|g′(u)− g′(0)| ≤ Lu‖d‖2, 0 ≤ u ≤ t ≤ σ∗.
Und wegen (6.19) erhalten wir
0 = g(σ∗) ≤ 0 + (1− α)(gradf(x0), d)σ∗ +12σ∗2L‖d‖2;
oder
(6.21) σ∗ ≥ −2(1− α)L
(gradf(x0), d)‖d‖2
.
Wegen g(σ∗) = 0 gilt auch
(6.22) f(x0 + σ∗d)− f(x0) ≤ σ∗α(gradf(x0), d)
und daher
(6.23) f(x0 + σ∗d)− f(x0) ≤ −2α(1− α)L
((gradf(x0), d)
‖d‖
)2
.
Das Verfahren von Armijo berechnet Schrittweiten σk, die Ungleichungender Form (6.21) (mit einer gewissen Konstanten, da man L nicht kennt) und(6.22) erfullen.
Ein erster Plan des Verfahrens lautet etwa so:
Sei (gradf(xk), dkN ) < 0, und 0 < α < 1, γ > 0, 0 < β1 < β2 < 1fest gewahlt.
6 Optimierung 107
1. Wahle eine Ausgangsschrittweite σ0 = σ0k mit
σ0 ≥ −γ(gradf(xk), dkN ) und setze j = 0.
2. Gilt dann
f(xk + σjdkN ) ≤ f(xk) + ασj(gradf(xk), dkN ),
so setze σk = σj , xk+1 = xk + σkdkN .
Sonst wahle σj+1 ∈ [β1σj , β2σ
j ] (σj+1 ≤ β2σj),
setze j = j + 1 und gehe zu 2.
Man kann zeigen, dass unter den Voraussetzungen (6.15),(6.20) die Schritt-weite σk nach endlich vielen Schritten gefunden wird (dies ist fur (6.22)klar). Unter diesen Voraussetzungen ist naturlich f(xk)∞k=0 konvergentund xk∞k=0 enthalt eine konvergente Teilfolge xkj → x. Wegen (6.23) gilt
Dnkf(xk) =(
gradf(xk),dkN‖dkN‖
)→ 0
(nk =
dkN‖dkN‖
).
Setzt man auch voraus, dass f ′′(x) positiv definit ist fur x ∈ N(x0), so giltlimk→∞ gradf(xk) = gradf(x) = 0. Ist dann noch x das einzige Minimum inN(x0), so konvergiert xk → x. Und man kann zeigen (s. Literatur), dass dasunter diesen (starken) Voraussetzungen das gedampfte Newton-Verfahrenfur große k σk = 1 ergibt; also quadratische Konvergenz.
Ein Nachteil dieses Verfahrens ist, dass man i.a. keine Abstiegsrichtungerhalt, wenn f ′′(xk) nicht positiv definit ist. Bei den Trust-Region-Verfahrenist auch dann noch ein Abstieg moglich. Zum Abschluss behandeln wir eineTrust-Region-Variante des Newton-Verfahrens.
Es ist fur f ∈ C2
f(xk + d) = f(xk) + f ′(xk)d+12
(d, f ′′(xk), d) + O(‖d‖2)
und daher fur kleine ‖d‖ der quadratische Ausdruck
(6.24) fk(d) = f(xk) + f ′(xk)d+12
(d, f ′′(xk)d)
eine gute Approximation von f(xk + d) auf der Trust-Region Bρk(0) =d | ‖d‖ ≤ ρk.
Wir suchen nun eine Abstiegsrichtung dk als Losung des quadratischen Op-timierungsproblems mit Nebenbedingung
(6.25) min‖d‖2≤ρk
fk(d) = fk(dk).
6 Optimierung 108
Dieses Problem ist naturlich stets losbar, wohingegen fk(d) → min i.A.unlosbar ist, und im Fall f ′′(xk) positiv definit die Losung dkN besitzt!
Es ist stets
(6.26) fk(dk) ≤ fk(0) = f(xk),
so dass durch Steuern von ρk mit einem Abstieg zu rechnen ist.
Ist fk(dk) = f(xk), so folgt
f(xk) ≤ f(xk) + f ′(xk)d+12
(d, f ′′(xk)d), ‖d‖ ≤ ρk,
oder0 ≤ f ′(xk)d+
12
(d, f ′′(xk)d), ‖d‖ ≤ ρk.
Also ist 0 ein lokales Minimum der Funktion rechts, und somit f ′(xk) = 0;(d, f ′′(xk)d) ≥ 0, d ∈ Rn. Daher verwendet man die Abfrage ”fk(d
k) =f(xk)?“ als Stoppkriterium.
Die Crux des Verfahrens ist naturlich die Losung von (6.25).
Wir formulieren zunachst das(6.27) Verfahren
Sei 0 < δ1 < δ2 < 1, σ1 ∈ (0, 1), σ2 > 1 und ρ0 > 0 gewahlt.
1. Wahle x0 ∈ Rn,berechne b0 = gradf(x0), A0 = f ′′(x0),und setze k = 0.
2. Berechne eine Losung dk des Problems
min‖d‖≤ρk
fk(d), fk(d) = f(xk) + (bk, d) +12
(d,Akd).
Ist f(xk) = fk(dk), so stoppt das Verfahren, und xk ist einKandidat fur ein lokales Minimum.
3. Ist f(xk) > fk(dk) (s. (6.26)), berechne
rk =f(xk)− f(xk + dk)f(xk)− fk(dk)
.
Ist rk ≥ δ1, so liegt ein Abstieg vor, setze
xk+1 = xk + dk
bk+1 = gradf(xk+1)Ak+1 = f ′′(xk+1).
6 Optimierung 109
Ist rk ≤ δ2 (d.h. rk ∈ [δ1, δ2]), so lasst man die Trust-Region(fast) unverandert, und das Modell fk(d) beschreibt dortf(xk + d) gut;man wahlt ρk+1 ∈ [σ1ρk, ρk].Ist sogar rk > δ2, dann kann man die Trust-Region ver-großern und wahlt ρk+1 ∈ [ρk, σ2ρk].Setze k = k + 1 und gehe zu 2.
4. Ist rk < δ1, dann wird nichts akzeptiert, die Trust-Regionist zu groß. Man wahlt ρk+1 ∈ (0, σ1ρk) und setzt
xk+1 = xk
bk+1 = bk
Ak+1 = Ak
k = k + 1
und geht zu 2.
Ist die Niveaumenge N(x0) kompakt, D offen mit N(x0) ⊂ D und f ∈C2(D), so kann man zeigen (s. Literatur), dass limk→∞ gradf(xk) = 0 gilt.Daher ist jeder Haufungspunkt der Folge xk∞k=0 ein stationarer Punkt.Ferner fallt das Verfahren (unter bestimmten Bedingungen) fur große k mitdem Newton-Verfahren zusammen.
Wir befassen uns noch mit der Losung des Problems (6.25).Sei also A ∈ Rn×n, At = A, b ∈ Rn und
(6.28) F (x) =12
(x,Ax) + (b, x); (†)
und zu losen ist das Problem
(6.29) min‖x‖2≤r
F (x).
Sei
(6.30) g(x) =12
(r2 − (x, x)).
Es gilt die folgende Charakterisierung eines Minimalpunktes x von F aufBr(0); dies sind die Karush-Kuhn-Tucker-Bedingungen aus der Optimie-rung, die in diesem Fall auch hinreichend sind.
Satz 6.7 x ist eine Losung des Problems (6.29) genau dann, wenn ein (ein-deutiges) µ ≥ 0 existiert, mit
(†)O.E. kann man A = At annehmen, da A = 12(A + At) + 1
2(A−At).
6 Optimierung 110
(i) (A+ µ)x+ b = 0,
(ii) µ(r2 − ‖x‖2) = 0, ‖x‖ ≤ r,
(iii) A+ µ ist p.s.d..
Beweis: ⇒: Ist ‖x‖2 < r, so ist gradF (x) = Ax + b = 0 und F ′′(x) = Ap.s.d. (Satz 6.1, 6.2). Also sind die Bedingungen (i)-(iii) fur µ = 0 erfullt.
Sei ‖x‖2 = r. Dann gilt nach Satz 6.1:
(1) (x, n) < 0⇒ (gradF (x), n) ≥ 0.
Ist (y, x) = 0, y 6= 0, so gilt furnk = y − 1
k x, (nk, x) = − 1k‖x‖
2 < 0,und daher (gradF (x), nk) ≥ 0, woraus(gradF (x), y) ≥ 0 folgt.Da dies auch fur −y gilt, haben wir
(y, x) = 0⇒ (gradF (x), y) = 0; d.h.
gradF (x) ∈ x | (x, x) = 0⊥ ≡< x > .
n
x
x
0
← z|(z, x) = 0
x− 2(x, x)x
Also existiert ein α ∈ R mit
(2) Ax+ b = αx.
Sei (x, n) < 0. Dann ist wegen (1) (Ax + b, n) = α(x, n) ≥ 0. Also istα = −µ, µ ≥ 0; und (i) ist bewiesen.
Multipliziert man (i) mit x, so folgt (Ax + b, x) = −µ‖x‖2, oder µ =− 1r2
(Ax+ b, x). µ ist eindeutig bestimmt.
(iii) erhalt man so: Sei ‖x‖2 = 1 und (x, x) 6= 0. Dann hat der gespiegelteVektor (s. Householder Verfahren)
(3) x− 2(x, x)x = x+ tx, t := −2(x, x)
dieselbe Lange wie x. Also gilt
0 ≤ F (x+ tx)− F (x) = t(Ax+ b, x) +12t2(x,Ax)
(i)= −tµ(x, x) +
12t2(x,Ax)
(3)=
12t2µ+
12t2(x,Ax)
=12t2(µ(x, x) + (x,Ax)) =
12t2((A+ µ)x, x).
Optimierung 111
⇐: Sei µ > 0. Dann ist ‖x‖ = r, d.h. g(x) = 0. Wir betrachten die Lagrange-Funktion zu x, d.h.
(4) L(x;x) = F (x)− µg(x) = −12µr2 + (b, x) +
12
(x, (A+ µ)x).
Da A+ µ p.s.d. ist L(x;x) konvex in x.Wegen gradxL(x;x)x=x = (A+ µ)x+ b = 0 ist x ein globaler Minimalpunktvon L(x;x) (s. Satz 6.6). also ist fur x ∈ R
(5) F (x) = L(x; x) ≤ L(x;x) = F (x)− µg(x).
Fur ‖x‖2 ≤ r, d.h. g(x) ≥ 0, erhalten wir aus (5) dann F (x) ≤ F (x), ‖x‖2 ≤r.
Sei µ = 0. Nach Voraussetzung ist A p.s.d.. Also ist F (x) konvex. x ist einglobaler Minimalpunkt von F in Rn. Also gilt F (x) ≤ F (x), x ∈ Rn. 2
Man lasse sich zur Losung von (i)-(iii) aus Satz 6.7 etwas einfallen.
A Lineare Differentialgleichungen 112
A Lineare Differentialgleichungen
Wir befassen uns jetzt noch weiter mit linearen Differentialgleichungen
(A.1) y′(x) = A(x)y(x) + b(x), x ∈ I,
wobei I ein offenes Intervall und A(x), b(x) ∈ IRn stetig sind.Ist b = 0, so heißt die Differentialgleichung homogen, sonst inhomogen.Die Menge der Losungen der homogenen Gleichung (A.1), d.h. y ∈ C1(I) |y′(x) = A(x)y(x), x ∈ I, bilden einen Vektorraum: Sind y1, y2 zwei Losun-gen der homogenen Gleichung, so ist auch α y1(x) + βy2(x) (α, β ∈ IR) eineLosung der homogenen Gleichung.Wir werden sehen, dass dieser Vektorraum endlich dimensional ist.Doch zunachst
Definition A.1 Die auf einem Intervall I definierten (vektorwertigen) Funk-tionen yi, i = 1, 2, . . . , p heißen linear unabhangig (l.u.), wenn aus
αi ∈ IR, i = 1, 2, . . . , p,p∑i=1
αiyi(x) = 0, x ∈ I, stets αi = 0, i = 1, 2, . . . , p,
folgt.Andernfalls heißen sie linear abhangig (l.a.).
Zum Beispiel sind 1, x, x2, . . . , xn linear unabhangig auf [a, b].
Lemma A.1 Es seien yi, i = 1, 2, . . . ,m Losungen der homogenen Glei-chung ((A.1) mit (b ≡ 0)). Ist m > n, so sind die yi, i = 1, 2, . . . ,m, linearabhangig.
Beweis: yi(x), i = 1, 2, . . . ,m, x ∈ I, seien m Losungen der homogenenGleichung. Fur jedes feste x0 ∈ I sind die Vektoren y1(x0), . . . , ym(x0) ∈ IRnlinear abhangig.
Also existiert ein (α1, . . . , αm)t 6= 0 mit 0 =m∑i=1
αiyi(x0). Die Funktion
y(x) =m∑i=1
αiyi(x) ist Losung der homogenen Gleichung, und es gilt y(x0) =
0. Wegen der Eindeutigkeit ist dann y(x) = 0, x ∈ I. 2
Lemma A.2 Es gibt (hochstens) n linear unabhangige Losungen der ho-mogenen Gleichung ((A.1) mit b ≡ 0).Ist insbesondere det(y1, . . . , y
n) 6= 0, so sind die n Losungen yi, i = 1, 2, . . . , n,
definiert durch y′i(x) = A(x) yi(x), x ∈ I, yi(x0) = yi linear unabhangig.
A Lineare Differentialgleichungen 113
Beweis: Sein∑i=1
αiyi(x) = 0, x ∈ I. Dann ist fur x = x0
n∑i=1
αi yi(x0) = 0 =n∑i=1
αi yi , woraus αi = 0, i = 1, 2, . . . , n folgt. 2
Definition A.2 Ein Fundamentalsystem (FS) zur homogenen Gleichung((A.1) mit b ≡ 0) besteht aus n linear unabhangigen Losungen y1, . . . , yn. Diezugehorige n×n Matrix W (x) = (y1(x), . . . , yn(x)) heißt Wronski-Matrix.
Lemma A.3 Es sei y1, . . . , yn ein Fundamentalsystem. Dann ist y eineLosung der homogenen Gleichung ((A.1) mit b ≡ 0) genau dann, wenn dieDarstellung
(A.2) y(x) =n∑i=1
ciyi(x) = W (x) c, x ∈ I
mit c = (c1, . . . , cn)t gilt.
Beweis: Ist y eine Losung der homogenen Gleichung, so sind nach Lem-ma A.1 die Funktionen y(x), y1(x), . . . , yn(x) linear abhangig auf I. Daherexistiert (α, α1, . . . , αn)t 6= 0 mit
α y(x) +n∑i=1
αi yi(x) = 0, x ∈ I.
Es ist α 6= 0. Ware namlich α = 0, so folgtn∑i=1
αi yi(x) = 0, x ∈ I, woraus
wegen der linearen Unabhangigkeit der yi folgt, dass ai = 0 gilt, was einWiderspruch zur Annahme ware. Also ist α 6= 0, und daher
y(x) = −n∑i=1
αiαyi(x), x ∈ I.
2
Kennt man eine spezielle Losung ys der inhomogenen Gleichung ((A.1),b 6≡ 0), so kennt man alle Losungen:
Satz A.1 Es sei ys eine spezielle Losung der inhomogenen Gleichung ( (A.1),b 6≡ 0), z.B. die Losung des AWP’s, und y1, . . . , yn sei ein Fundamentalsy-stem.Dann besitzt jede Losung y der inhomogenen Gleichung die Darstellung(A.3)
y(x) =n∑i=1
ciyi(x)+ys(x) = W (x) c+ys(x), x ∈ I mit c = (c1, . . . , cn)T .
A Lineare Differentialgleichungen 114
Beweis: Da die Funktion Z(x) = y(x) − ys(x) die homogene Gleichung((A.1), b 6≡ 0) erfullt, folgt mit Lemma A.4 die Behauptung: 2
Satz A.2 Sei W Wronski-Matrix zu einem Fundamentalsystem der Diffe-rentialgleichung (A.1).Dann gilt:
(i) W ′(x) = A(x)W (x), x ∈ I,
(ii) detW (x) 6= 0, x ∈ I, sgn det W (x) = const, x ∈ I
(iii) W (x), x ∈ I ist eine Wronski-Matrix genau dann, wenn eine nicht-singulare Matrix B existiert mit
(A.4) W (x) = W (x)B, x ∈ I.
Beweis:
(i) Trivial.
(ii) Existiert etwa ein x1 ∈ I mit det W (x1) = 0, so sind die Vektoreny1(x1), . . . , yn(x1) linear abhangig. Also existiert (α1, . . . , αn)t 6= 0 mit
0 =n∑i=1
αiyi(x1). Die Funktion y(x) =n∑i=1
αiyi(x) erfullt y′ = Ay und
y(x1) = 0; also ist y(x) ≡ 0, woraus αi = 0 fur i = 1, 2, . . . , n folgt; dasist ein Widerspruch. Also besitzt detW (x) keine Nullstellen in I. DadetW (x) stetig ist, kann detW (x) in I nicht das Vorzeichen wechseln(Zwischenwertsatz!).
(iii) Sei W eine beliebige Wronski-Matrix. Setzt man C(x) = W (x) −W (x)B mitB = W−1(x0) W (x0) (x0 ∈ I), so gilt C ′(x) = A(x)C(x), x ∈I, C(x0) = 0, woraus C(x) ≡ 0 folgt.Umgekehrt ist jede Matrix der Form (A.4) (mit detB 6= 0) eineWronski-Matrix. 2
Die in (A.3) benutzte spezielle Losung ergibt sich auch mit Hilfe des Funda-mentalsystems, und zwar durch die Methode der ”Variation der Konstan-ten“. Hierbei macht man fur die spezielle Losung ys den Ansatz ys(x) =n∑i=1
ci(x)yi(x) = W (x) c(x) (Es wird der “Vektor c variiert“).
Es folgt
y′s(x) = W ′(x) c(x)+W (x) c′(x) != A(x) ys(x)+b(x) = A(x)W (x) c(x)+b(x),
A Lineare Differentialgleichungen 115
oder wegenW ′ = AW ist dannW (x) c′(x) = b(x), oder c′(x) = W−1(x) b(x).Setzt man daher
c(x) =∫ x
x0
W−1(t) b(t) dt, d.h. c(x0) = 0,
so istys(x) = W (x)
∫ x
x0
W−1(t) b(t) dt
eine spezielle Losung mit ys(x0) = 0.
Satz A.3 Sei W Wronski-Matrix. Dann besitzt jede Losung von y′(x) =A(x)y(x) + b(x), x ∈ I die Darstellung
(A.5) y(x) = W (x) c+∫ x
x0
W (x)W−1(t) b(t) dt, x ∈ I
mit einem Vektor c ∈ IRn.Die Losung des AWP’s lautet
(A.6) y(x) = W (x)W−1(x0) y0 +∫ x
x0
W (x)W−1(t) b(t) dt.
Beweis: (A.5) ist klar (Satz A.1), und (A.6) folgt daraus. 2
Beispiel:
(A.7) x′′(t) + w2x(t) = f(t), t ≥ 0
Diese Differentialgleichung 2. Ordnung (Oszillator) ist aquivalent zu demSystem erster Ordnung
y′(t) = A(t)y(t) +(
0f(t)
), t ≥ 0,
mit der Wronski-Matrix
W (t) =(
sin(wt) cos(wt)w cos(wt) −w sin(wt)
).
Es ist W−1(t) = − 1w
(−w sin(wt) − cos(wt)−w cos(wt) sin(wt)
),
A Lineare Differentialgleichungen 116
woraus
ys(t) =
t∫0
(cosw(t− u) 1
w sinw(t− u)−w sinw(t− u) cosw(t− u)
)(0
f(u)
)du =
(xs(u)x′s(u)
)
folgt. Daher ist eine spezielle Losung von (A.7) gegeben durch
xs(t) =∫ t
0
1w
sinw(t− u)f(u) du,
und die allgemeine Losung ( y(t) =(
x(t)x′(t)
)= W (t)c+ ys(t)) ergibt daher
die Darstellung
x(t) = c1 sinwt+ c2 coswt+∫ t
0
1w
sinw(t− u)f(u) du
mit beliebigen Konstanten c1, c2.
Fur eine lineare Anfangswertaufgabe n-ter Ordnung (n ≥ 2)
(A.8)y(n)(x) +
n−1∑k=0
Ak(x) y(k)(x) = f(x), x ∈ I, Ak, f ∈ C(I)
y(k)(x0) = γk, k = 0, 1, . . . , n− 1
erhalt man durch Umwandlung in ein System erster Ordnung (setze z(x) =(y(x), . . . , y(n−1)(x)
), s. Lemma 5.1) analoge Resultate:
(A.9)y(x)
...
...y(n−1)(x)
′
=
0 1
. . .1
−A0 . . . . . . −An−1(x)
y(x)......
y(n−1)(x)
+
0...0
f(x)
Zum Schluss dieses Abschnitts befassen wir uns noch mit linearen Differen-tialgleichungen mit konstanten Koeffizienten
(A.10) y′(t) = Ay(t) + b(t), t ∈ IR,
wobei A eine n×n-Matrix ist und b ∈ C(IR), bzw.
(A.11) y(n)(t) +n−1∑k=0
Ak y(k)(t) = f(t), t ∈ IR,
A Lineare Differentialgleichungen 117
oder aquivalent
( A.11)’ z′(t) =
0 1
. . .1
−A0 · · · · · · −An−1
z(t)+
0......
f(t)
mit reellen Zahlen Ak und f ∈ C (IR) ist. Die Matrix
B =
0 1
. . .1
−A0 · · · · · · −An−1
heißt Begleitmatrix (zu dem Polynom xn +
n−1∑k=0
Akxk).
Unser Ziel ist es nun, aussagekraftige Fundamentalsysteme fur (A.11), (A.11)’zu konstruieren.Diese werden sich aus den Eigenwerten der Matrix A bzw. aus denen der Be-gleitmatrix B ergeben. Nun sind aber die Eigenwerte von B die Nullstellendes charakteristischen Polynoms
(A.12) pn(x) = xn +n−1∑k=0
Akxk,
denn es gilt
(A.13) det(B − λI) = ±pn(λ),
was sich z.B. durch Entwicklung nach der letzten Zeile in det(B−λI) ergibt.
Sei y′i(t) = Ayi(t), i = 1, 2, . . . , n, x ∈ IR.Nach Lemma A.3 bilden diese yi(t) ein Fundamentalsystem genau dann,wenn
(A.14) det(y1(0), . . . , yn(0)) 6= 0 ist.
Macht man fur die Losung der homogenen Gleichung den Ansatz y(t) =ceλ·t (c 6= 0), so gilt y′(t) = Ay(t) ⇔ Ac = λc. D.h., λ ist Eigenwert vonA, und c ist ein Eigenvektor dazu. (Diese Großen konnen naturlich komplexsein!) Da y(0) = c ist, gilt daher
A Lineare Differentialgleichungen 118
Lemma A.4 Besitzt die Matrix A ∈ IRn×n n linear unabhangige Eigenvek-toren c1, . . . , cn (∈ ICn) zu Eigenwerten λ1, . . . , λn, so bilden die Funktionen
(A.15) yi(t) = cieλit
ein Fundamentalsystem von (A.10).
Matrizen, die die Voraussetzungen von Lemma A.4 erfullen, heißen diago-nalisierbar. Ist namlich T die Matrix der Eigenvektoren von A
(A.16) T = (c1, . . . , cn),
so gilt A · T = T · diag(λ1, . . . , λn), oder
(A.17) T−1AT = diag(λ1, . . . , λn).
Setzt man daher y(t) = Tz(t), so geht y′ = Ay uber in das entkoppelteSystem
(A.18) z′(t) = diag(λ1, . . . , λn)z(t)
mit dem Fundamentalsystem zi(t) = eieλit, ei = (0, . . . , i, . . . , 0)t, worausdurch Rucktransformation wieder (A.15) folgt.Eine große Klasse von diagonalisierbaren Matrizen sind die normalen Ma-trizen: Wir nennen eine Matrix A normal, falls AtA = AAt gilt. Hier istsogar (ci, ck) = δik, d.h. T ist unitar.Fur symmetrisches A sind λi ∈ IR, ci ∈ IRn, und T ist eine orthogonaleMatrix (T tT = 1).
Um den allgemeinen Fall zu behandeln, benotigen wir weitere Tatsachen ausder Algebra. Sei λ ein EW von A, dann sind die Nullraume
(A.19) N((A− λ)j) = x | (A− λ)jx = 0, j = 1, 2, . . .
nicht leer, und es gilt N((A− λ)j) ⊂ N((A− λ)j+1), j = 1, 2, . . . .Wegen dim N((A− λ)j) ≤ n existiert ein j0 ∈ IN mit
(A.20) N((A− λ)j0) = N((A− λ)j), j ≥ j0.
Das kleinste j0 mit (A.20) heißt die Nullkettenlange von λ, α(λ) (≡Elementarteilerexponent!) Fur α(λ) ≥ 1 existiert daher ein x mit (A −
A Lineare Differentialgleichungen 119
λ)α−1x 6= 0, und (A − λ)αx = 0. (Denn sonst ware N((A − λ)α−1) =N((A − λ)α) = N((A − λ)l), l ≥ α.) Ein solches x ist ein Hauptvektor(HV) α-ter Stufe. Ein Hauptvektor x l-ter Stufe (l ≥ 1) ist auch aquivalentdefiniert durch
(A.21)(A− λ)jx 6= 0, j = 0, 1, . . . , l − 1
(A− λ)lx = 0.
In der Algebra zeigt man:
(A.22) dimN((A− λ)α(λ)) = m(λ),
wobei m die algebraische Vielfachheit von λ ist (d.h. λ ist eine m-facheNullstelle des charakteristischen Polynoms det(A− xI)), und
(A.23) IRn =s⊕i=1
N(
(A− λi)α(λi)),
wobei λi die verschiedenen EW von A sind. D.h., jedes x ∈ IRn ist eindeu-tig darstellbar in der Form
x =s∑i=1
zi mit zi ∈ N(
(A− λi)α(λi)).
Manchmal zeigt man in der Algebra sogar:
(A.24) N((A− λ)α(λ)) besitzt eine Jordan-Basis:
Es sei dim N(A−λ) = ρ = maximale Zahl linear unabhangiger Eigenvekto-ren zum Eigenwert λ. ρ ist die geometrische Vielfachheit von λ, ρ ≤ m.Es gibt dann ρ linear unabhangige Eigenvektoren x1, . . . , xρ von λ mit denfolgenden Eigenschaften:Fur jedes i = 1, 2, . . . ρ gibt es Vektoren yli mit(A.25)y1i = xi , (A− λ)yl+1
i = yli, l = 1, 2, . . . (Mi − 1) (Jordan-Kette)
Dann ist jedes yli, l = 1, 2, . . . ,Mi, ein Hauptvektor der Stufe l. Es istmax
i=1,2,...,ρMi = α(λ), und die Gesamtheit dieser Vektoren bilden eine Ba-
sis fur N((A− λ)α(λ)
), m(λ) =
ρ∑i=1
Mi.
Wegen (A.23) bilden daher die Haupt- und Eigenvektoren der Matrix A eineBasis fur den IRn. Fasst man diese in einer Matrix T zusammen (s. (A.16)),
A Lineare Differentialgleichungen 120
so gibt bekanntlich die Ahnlichkeitstransformation (A.17) die JordanscheNormalform.
Warnung (J. Grande): Die Jordanketten mussen nicht fur jede Basis vonN(A− λ) existieren! Zum Beispiel fur
A =
1 1 00 1 00 0 1
ist λ = 1 dreifacher Eigenwert vonA und die Eigenvektoren
101
,
10−1
bilden eine Basis von N(A− 1). Das System
(A− 1)y =
10±1
ist aber nicht losbar. Das System (A−1)y = α
101
+β
10−1
muss fur
gewisse α, β losbar sein, namlich α = β; d.h. fur den Eigenvektor
100
.
Man erhalt y =
010
.
Damit konnen wir den allgemeinen Fall behandeln:
Satz A.4 A ∈ IRn×n besitze s paarweise verschiedene Eigenwerte λi, i =
1, 2, . . . , s, mit der algebraischen Vielfachheit m(λi) (det(A−x) = ±s∏i=1
(x−
λi)m(λi)) und der geometrischen Vielfachheit ρ(λi) ≤ m(λi). Dann bildendie folgenden n Funktionen ein Fundamentalsystem fur y′(t) = Ay(t) :Fur λ ∈ λ1, . . . , λs mit den Vielfachheiten m, ρ bilde
(A.26) yi(t) = xieλt, i = 1, 2, . . . , ρ,
wobei xi die Eigenvektoren zum Eigenwert λ sind. Und fur jeden Hauptvek-tor c der Stufe l ≥ 2 aus den Jordan-Ketten (A.25), c = yli bilde
(A.27) yc(t) = eλtl−1∑k=0
1k!
(A− λ)kctk.
A Lineare Differentialgleichungen 121
Beweis: Wegen yi(0) = xi, yc(0) = c, und da die Gesamtheit dieser Vekto-ren eine Basis fur den IRn bilden, brauchen wir nur zu zeigen, dass (A.27)die homogene Gleichung y′ = Ay erfullt:
y′c(t) = λeλtl−1∑k=0
1k!(A− λ)kctk + eλt
l−1∑k=1
1k!(A− λ)kcktk−1
= λeλtl−1∑k=0
1k!(A− λ)kctk + eλt
l−2∑k=0
1k!(A− λ)k+1ctk.
Ayc(t) = eλtl−1∑k=0
1k!(A− λ+ λ)(A− λ)kctk
= λeλtl−1∑k=0
1k!(A− λ)kctk + cλt
l−1∑k=0
1k!(A− λ)k+1ctk
= λeλtl−1∑k=0
1k!(A− λ)kctk + cλt
l−2∑k=0
1k!(A− λ)k+1ctk,
da (A− λ)lc = 0. 2
Bemerkung: Formal ist
yc(t) = eλtl−1∑k=0
1k!
(A− λ)kctk = eλt∞∑k=0
1k!
(A− λ)kctk
= eλte(A−λ)t · c = eAt · c;
d.h. yc(t) = eAt · c
Es gilt (wegen Akxi = λkxi): yi(t) = eAtxi. Man kann dies auch exaktbegrunden, und zeigen, dass eAt (uber die unendliche Reihe definiert) eineWronski-Matrix ist.Wir schließen das Kapitel mit einem Beispiel:
(A.28) y′(t) = Ay(t) + b(t), A =
1 0 07 1 00 0 2
p(x) = det(A− x) = −(x− 1)2(x− 2).
Also ist λ1 = 1, m(λ1) = 2, λ2 = 2, m(λ2) = 1.Zum Eigenwert λ1 = 1 gibt es nur einen Eigenvektor x1, da
(A− 1)x1 = 0 =
07x1
1
x13
⇒ x11 = x1
3 = 0, x1 =
010
.
A Lineare Differentialgleichungen 122
(A− 1)y21 = x1 ergibt
y21 =
1760
als Hauptvektor 2. Stufe,
und der Eigenvektor zum Eigenwert λ2 = 2 ist z1 =
001
.
Das Fundamentalsystem besteht daher aus
y1(t) =
010
et, y3(t) =
001
e2t
und
y2(t) = et(y21 + (A− λ)y2
1t) = et(y21 + x1 · t) =
17
6 + t0
et.
Mit T = (x1, y21, z
1) findet man wegen
AT = (x1, x1 + y21, 2z1) = (x1, y2
1, z1)
1 1 00 1 00 0 2
die Jordansche Normalform
J = T−1AT =
1 1 00 1 00 0 2
.
Literatur 123
Literatur
[1] Amann, H.: Gewohnliche Differentialgleichungen. Walter de Gruyter,Berlin, 1983.
[2] Arnold, V.I.: Gewohnliche Differentialgleichungen. Springer Verlag,Berlin, 1980.
[3] de Boor, C.: A Practical Guide to Splines. Springer Verlag, New York,1978.
[4] Braun, M.: Differential Equations and their Applications. SpringerVerlag 1983.
[5] Coddington, E. A. und N. Levinson: Theory of Ordinary Differen-tial Equations. McGraw–Hill, New York, 1955.
[6] Collatz, L.: Funktionalanalysis und Numerische Mathematik. Sprin-ger Verlag, Berlin 1964.
[7] Cronin, J.: Differential Equations. Dekker, New York, 1980.
[8] Dahmen, W. und A. Reusken: Numerik fur Ingenieure und Natur-wissenschaftler. Springer, Berlin, 2008.
[9] Deuflhard, P. und A. Hohmann: Numerische Mathematik I. Walterde Gruyter, Berlin, 1993.
[10] Deuflhard, P. und F. Bornemann: Numerische Mathematik II.Walter de Gruyter, Berlin, 1994
[11] Engeln-Mullges, E. und F. Reutter: Formelsammlung zur Nume-rischen Mathematik. 6. Auflage, BI-Verlag, Mannheim, 1988.
[12] Erwe, F.: Gewohnliche Differentialgleichungen. BI-Verlag, Mannheim,1964.
[13] Franklin, J.N.: Matrix Theory. Practice Hall, Englewood Cliffs, NJ,1968.
[14] Golub, G.H. und C.F. van Loan: Matrix Computations. 2. Auflage,The Johns Hopkins University Press, Baltimore, MD, 1989.
[15] Golub, G.H. und J.M. Ortega: Scientific Computation and Diffe-rential Equations. Academic Press, Boston, MA, 1992.
[16] Grigorieff, R.D.: Numerik Gewohnlicher Differentialgleichungen 1.Teubner Verlag, Stuttgart, 1987.
Literatur 124
[17] Hammerlin, G. und K.-H. Hoffmann: Numerische Mathematik.Springer Verlag, Berlin, 1989.
[18] Hairer, E.,S.P. Nørsett und G. Wanner: Solving Ordinary Diffe-rential Equations I. Springer Verlag, Berlin, 1987.
[19] Hale, J.: Ordinary Differential Equations. John Wiley and Sons,New York, 1969.
[20] Hartmann, P.: Ordinary Differential Equations. John Wiley and Sons,New York, 1964.
[21] Hestenes, M.R. und E. Stiefel: Methods of conjugate gradients forsolving linear systems, Jour. Res. National Bureau of Standards, 29(1952), pp 409-439.
[22] Hille, E.: Lectures on Ordinary Differential Equations. Addision–Wesley, Reading, MA, 1968.
[23] Hirsch, M. und S. Smale: Differential Equations, Dynamic Systems,and Linear Algebra. Academic Press, New York, 1974.
[24] Kamke, E.: Differentialgleichungen: Losungsmethoden und LosungenI. Teubner Verlag, Stuttgart, 1977.
[25] Kamke, E. : Differentialgleichungen I, Gewohnliche Differentialglei-chungen. Akademische Verlagsgesellschaft Geest u. Portig, Leipzig,1964.
[26] Knobloch, H.W. und F. Kappel: Gewohnliche Differentialgleichun-gen . Teubner Verlag, Stuttgart, 1974.
[27] Locher, F.: Numerische Mathematik fur Informatiker. Springer Ver-lag, Berlin, 1991.
[28] Luenberger, D.: Introduction to Linear and Nonlinear Programming.Addison–Wesley, Reading, MA, 1973.
[29] Meinardus, G. und G. Merz: Praktische Mathematik I. Bibliogra-phisches Institut, Mannheim, 1979.
[30] Natanson, I.P.: Konstruktive Funktionentheorie I, II. Akademie Ver-lag, Berlin, 1955.
[31] Niederdrenk, K.: Die endliche Fourier– und Walsh–Transformationmit einer Einfuhrung in die Bildverarbeitung. Vieweg Verlag, Braun-schweig, 1982.
[32] Peyerimhoff, A.: Gewohnliche Differentialgleichungen I, II. Studien-text, Akademische Verlagsgesellschaft 1970.
Literatur 125
[33] Press, H. et al.: Numerical Recipes, The Art of Scientific Computation.Cambridge University Press, New York, 1988.
[34] Reimer, M.: Grundlagen der Numerischen Mathematik I. Akademi-sche Verlagsgesellschaft, Wiesbaden, 1980.
[35] Rutishauser, H.: Vorlesungen uber Numerische Mathematik.Birkhauser Verlag, Basel, 1976.
[36] Schabak, R.: Numerische Mathematik II. 4. Auflage, Springer Verlag,Berlin, 1991.
[37] Schafke, F.W. und D. Schmidt: Gewohnliche Differentialgleichun-gen. Springer Verlag, Berlin, 1973.
[38] Schmeißer, G. und H. Schirmeier: Praktische Mathematik. Walterde Gruyter, Berlin, 1976.
[39] Schumaker, L.L.: Spline Functions: Basic Theory. John Wiley andSons, New York, 1981.
[40] Schwarz, H.R.: Numerische Mathematik. Teubner Verlag, Stuttgart,1988.
[41] Stoer, J.: Numerische Mathematik I, II, V bzw.III. 3. Auflage, Sprin-ger Verlag, Berlin, 1990.
[42] Stummel, F. und K. Hainer: Praktische Mathematik. Teubner Stu-dienbucher, Teubner Verlag, Stuttgart, 1971.
[43] Wendroff, B.: Theoretical Numerical Analysis. Academic Press,New York, 1966.
[44] Wilkinson, J.H.: Rundefehler. Springer Verlag, New York, 1969.
[45] Wilkinson, J.H.: The Algebraic Eigenvalue Problem. Clarendon Press,Oxford, 1965.
Index
Abbruchfehlerlokaler, 83
Absorptionsfunktion, 2Abstiegsrichtung, 105Abstiegsverfahren, 103Aitken-Neville-Algorithmus, 53algebraische Vielfachheit, 119Algorithmus
Aitken-Neville-, 53Gauß-, 5Stabilitat, 18
Anfangswertproblem, 65, 69lokale Losung, 69maximale Losung, 69n-ter Ordnung, 116
AufwandAitken-Neville, 54Cholesky-Zerlegung, 14Einschrittverfahren, 88Gauß-Algorithmus, 8Gesamtschrittverfahren, 23Householder-Verfahren, 37
Ausloschung, 37
B-Splines, 56Basis
Orthonormal-, 32, 37Begleitmatrix, 117Bestapproximation, 26
charakteristisches Polynom, 117Cholesky-Verfahren, 15Cholesky-Zerlegung, 14Cramersche Regel, 5
DarstellungLagrange-, 52Newton-, 54
Diagonalmatrix, 41Differentialgleichung, 65
Fundamentalsystem, 113homogene, 112
inhomogene, 112lineare Systeme, 79lineare, Hauptsatz, 80m-ter Ordnung, 68partielle, 66steife, 96System, 66
Differentiation, numerische, 61Diskretisierung, 67dividierte Differenz, 55Drehung, 35Dreiecksmatrix
obere, 7, 32untere, 8
Einschrittverfahren, 84Euler-Cauchy-, 82Runge-Kutta-, 88verbessertes Euler-Cauchy-, 85zweiter Ordnung, 85
Einzelschrittverfahren, 24Euler-Cauchy-Verfahren, 82
Fehlerabschatzung, 84Konvergenz, 83, 85lokaler Abbruchfehler, 83verbessertes, 85
Euler-McLaurin-Summenformel, 61
Fehlerrelativer, 16
Fehlerabschatzunga posteriori, 46a priori, 46
Fehlerdarstellung, 52Fehlerquadratsumme, 26finiter Ausdruck, 61Fixpunkt, 43
-gleichung, 43-iteration, 43-satz, Banachscher, 45
Fundamentalsystem, 113Funktion
126
Index 127
Hut-, 57konvexe, 101Rosenbrock-, 100
FunktionenLineare Abhangigkeit, 112
Gauß-Algorithmus, 5Gesamtschrittverfahren, 23
Konvergenz, 24Gewicht, 59Gitter, 83
aquidistandes, 95-funktion, 83
Gleichungssystemlineares, 4nichtlineares, 43
Gleitkommadarstellung, 15Gram-Schmidtsches Orthogonalisierungs-
verfahren, 32Gronwall-Lemma, 75
Hauptvektor, 119Horner-Schema, 54Householder-Verfahren, 32Hutfunktion, 57
Instabilitatnumerische, 17
Integralgleichung, 70Interpolation, 50Iterationsverfahren, 21, 22
Jordan-Basis, 119Jordansche Normalform, 120
Karush-Kuhn-Tucker-Bedingungen, 109klassisches Runge-Kutta-Verfahren, 92Kondition, 16, 29
Normalgleichungen, 30Konsistenz, 84Kontraktion, 44Konvergenz
Iterationsverfahren, 21Banachscher Fixpunktsatz, 45Einzelschrittverfahren, 25
Gesamtschrittverfahren, 24, 25quadratische, 47
Konvergenzordnung, 46konvexes Problem, 101Konvexitat, 101, 102kritischer Punkt, 27, 99, 103
Lagrange-Darstellung, 52-Funktion, 111-Grundpolynome, 61-Polynom, 51
Levenberg-Marquardt-Verfahren, 104lineares Gleichungssystem, 4Linienmethode, 67Lipschitzbedingung, 74
lokale, 74lokale Losung, 69
Eindeutigkeit, 75Existenz, 73, 76
lokaler Abbruchfehler, 83LR-Zerlegung, 8, 9
Maschinengenauigkeit, 15Maschinenzahlen, 15Matrix
dunnbesetzte, 20diagonaldominante, 18diagonale, 41diagonalisierbare, 118normale, 118orthogonale, 11, 32, 35, 41Permutations-, 11, 12spd, 13–15symmetrische, 35Vandermonde-, 51Wronski-, 113
maximale Losung, 69Existenz, 79, 80Randverhalten, 70, 72
Minimalpunktlokaler, 27
Minimumglobales, 99, 103
Index 128
lokales, 99, 103striktes, 99
Multigridverfahren, 25
Newton-Darstellung, 54Newton-Verfahren, 47
gedampftes, 104, 105, 107vereinfachtes, 47
NormEuklidische, 22Operator-, 22Spaltensummen-, 22Zeilensummen-, 22
Normalgleichung, 27, 28, 38Nullkettenlange, 118Nullstelle
Einfachheit der, 48Nullstellenproblem, 43numerische Differentiation, 58, 61numerische Instabilitat, 17numerische Integration, 58
Obere Dreiecksmatrix, 7, 32Optimalitatskriterium, 99Optimierung, 99orthogonale Matrix, 11orthogonale Projektion, 28Orthonormalbasis, 37Oszillator, 115
Peano-Existenzsatz, 73Permutationsmatrix, 11, 12Picard-Lindelof, Satz von, 76Pivotelement, 17polynomiales Interpolationsproblem,
50Polynominterpolation, 50Projektion
orthogonale, 28Pseudoinverse, 40
QR-Zerlegung, 32Quadraturfehler, 60Quadraturformel
interpolatorische, 59
zusammengesetzte, 60Quadraturverfahren, 58
ruckwartiges Euler-Verfahren, 92Ruckwartseinsetzen, 8, 10Radon-Transformation, 3Rand, 70relativer Fehler, 16Residuum, 26Riemann-Integrierbarkeit, 69Rosenbrock-Funktion, 100Runge-Kutta-Gauß-Verfahren, 92Runge-Kutta-Verfahren, 88
Durchfuhrbarkeit, 89explizites, 90halbimplizites, 90klassisches , 92Konsistenz, 91lokaler Abbruchfehler, 90m-stufiges, 90, 94Verfahrensfunktion, 91vollimplizites, 90
Sattelpunkt, 100schlecht konditioniert, 17Schrittweitensteuerung, 88, 93, 103Schwachungskoeffizient, 2Singularwert, 40
Zerlegung-, 41Singularwertzerlegung (SVD), 41Spaltenpivotsuche, 18Spaltensummenkriterium, 24Spektralradius, 21Spiegelung, 35Spline, 50
-Funktionen, 55B-, 56
Storung der Daten, 29Stutzstelle, 59Stutzstellenpolynom, 51System erster Ordnung, 67
Trager, 57Transformation
Radon-, 3
Index 129
Trapezregel, 60, 92Fehlerdarstellung, 60summierte, 61
Trust-Region-Verfahren, 104, 107tuning-Faktor, 96
Untere Dreiecksmatrix, 8
Vandermonde-Matrix, 51Variation der Konstanten, 114verbessertes Euler-Cauchy-Verfahren,
85Verfahren
Abstiegs-, 103Cholesky-, 15Einzelschritt-, 24Euler-Cauchy-, 82gedampftes Newton-, 104, 105Householder-, 32Levenberg-Marquardt-, 104Multigrid-, 25Newton-, 47Orthogonalisierungs-
Gram-Schmidtsches, 32Quadratur-, 58ruckwartiges Euler-, 92Runge-Kutta-, 88Runge-Kutta-Gauß-, 92Trust-Region-, 104, 107von Armijo, 106
Verfahrensfunktion, 86Vorwartseinsetzen, 10
Warmeleitungsgleichung, 66, 97Wachstumsrate, 65Wronski-Matrix, 113
Zeilenaquilibrierung, 17Zeilensummenkriterium, 24Zeilenvertauschung, 11Zerlegung
Cholesky-, 14LR-, 8, 9QR-, 32, 37Singularwert-, 41