Numerische Mathematik I

234 Pages·2008·0.877 MB·German

by  Peter Knabner

Numerische Mathematik I
Peter Knabner
Sommersemester 2008
¨
Uberarbeitete Fassung: 30.04.2008
1
Inhaltsverzeichnis
1 Direkte Verfahren fu¨r lineare Gleichungssysteme 4
1.1 Das Eliminationsverfahren von Gauß . . . . . . . . . . . . . . . . . . . 4
1.2 Dreieckszerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4 Weitere Formen von Eliminationsverfahren und das Cholesky Verfahren 23
2 Fehleranalyse und St¨orungsrechnung 28
2.1 Darstellung von Zahlen und Maschinenzahlen . . . . . . . . . . . . . . 28
2.2 Rundungsfehler und Gleitpunktarithmetik . . . . . . . . . . . . . . . . 32
2.3 Exkurs: Normen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.4 Kondition eines Problems – Normweise Konditionsanalyse . . . . . . . 48
2.5 Normweise Konditionsanalyse bei der L¨osung linearer Gleichungssysteme 52
2.6 Komponentenweise Konditionsanalyse . . . . . . . . . . . . . . . . . . . 57
2.7 Schlechtgestellte Probleme . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.8 Die Kondition (Stabilit¨at) von Algorithmen . . . . . . . . . . . . . . . 63
3 Orthogonalisierungsverfahren 71
3.1 QR–Zerlegung durch Househoulder–Transformationen . . . . . . . . . . 72
3.2 Lineare Ausgleichsrechnung . . . . . . . . . . . . . . . . . . . . . . . . 78
3.3 Pseudoinverse, Singul¨arwertzerlegung . . . . . . . . . . . . . . . . . . . 86
4 Nichtlineare Gleichungen und Gleichungssysteme 95
4.1 Fixpunktiteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.2 Das Newtonverfahren und Varianten . . . . . . . . . . . . . . . . . . . 102
4.3 Nichtlineare Ausgleichsprobleme: Das Gauß–Newton Verfahren . . . . . 111
4.4 Lineare Iterationsverfahren fu¨r lineare Gleichungssysteme . . . . . . . . 115
4.5 Das Newtonverfahren Teil 2: Modifikationen . . . . . . . . . . . . . . . 124
2
5 Interpolation 131
5.1 Lagrange’sche Polynominterpolation . . . . . . . . . . . . . . . . . . . . 133
5.2 Hermite–Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
5.3 Spline–Funktionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.4 Trigonometrische Interpolation . . . . . . . . . . . . . . . . . . . . . . . 159
6 Numerische Integration (Quadratur) 167
6.1 Newton–Cotes–Formeln . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
6.2 Die Euler–MacLaurin’sche Summenformel . . . . . . . . . . . . . . . . 174
6.3 Extrapolation. Das Romberg–Verfahren . . . . . . . . . . . . . . . . . . 177
6.4 Gauß–Quadratur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
7 Nichtlineare Optimierung ohne Nebenbedingungen 191
7.1 Nichtlineare Optimierung (mit linearen Gleichungsbedingungen) . . . . 191
7.2 Eindimensionale Optimierung . . . . . . . . . . . . . . . . . . . . . . . 194
7.3 Zusammenhang Optimierung – Gleichungssysteme, Newtonverfahren
und Varianten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197
7.3.1 Verfahren fu¨r station¨are Punkte . . . . . . . . . . . . . . . . . . 198
7.4 Konvexe und quadratische Funktionale . . . . . . . . . . . . . . . . . . 199
7.5 Gradientenverfahren und
Methode der konjugierten Gradienten fu¨r quadratische Optimierungs-
probleme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
7.6 Exkurs: Grosse du¨nnbesetzte lineare Gleichungssysteme aus der Dis-
kretisierung partieller Differentialgleichungen und die Effizienz von
L¨osungsverfahren dafu¨r . . . . . . . . . . . . . . . . . . . . . . . . . . . 213
7.7 Vorkonditionierte CG-Verfahren . . . . . . . . . . . . . . . . . . . . . . 218
7.8 Krylov-Unterraum-Methoden fu¨r nichtsymmetrische Gleichungssysteme 225
7.9 Verfahren der konjugierten Gradienten fu¨r nichtquadratische Optimierung231
Literaturverzeichnis 234
3
1 Direkte Verfahren fu¨r lineare Gleichungssysteme
1.1 Das Eliminationsverfahren von Gauß
Notation:
Rm,n := R–Vektorraum der m×n–Matrizen (m Zeilen, n Spalten)
A ∈ Rm,n : A = (a ) = (a )
ij ij i=1,...,m
j=1,...,n
AT := (a ) ∈ Rn,m transponierte Matrix
ji
Rn := Rn,1 := R–Vektorraum der Vektoren mit n Komponenten
(beachte: Vektoren = Spalten)
x ∈ Rn : x = (x ,...,x )T bzw. xT = (x ,...,x ) .
1 n 1 n
Komponenten von Matrizen und Vektoren tragen also tiefgestellte Indizes.
Indizierung von Vektoren durch hochgestellte Indizes, insbesondere bezeichne
ei den i-ten Einheitsvektor, d.h.
(cid:40)
1 fu¨r i = j
ei := δ := .
j ij
0 fu¨r i (cid:54)= j
Fu¨r u,v ∈ Rn ist das Skalarprodukt
uTv ∈ R
definiert durch
(cid:88)n
uTv = u v .
i i
i=1
Fu¨r u ∈ Rm, v ∈ Rn ist das dyadische oder ¨außere Produkt
uvT ∈ Rm,n
definiert durch
uvT = (u v ) .
i j
A ∈ Rm,n : A = (a1,...,an), ai ∈ Rm Spalten von A
(cid:88)n
E = E := (e1,...,en) = ekekT ∈ Rn,n Einheitsmatrix
n
k=1
4
A ∈ Rm,n : Aej : j-te Spalte von A
eiTA : i-te Zeile von A
eiTAej= a =: (A) .
ij ij
A kann als Summe der Spalten dargestellt werden:
(cid:88)n (cid:88)n
A = AE = A ekekT = (Aek)ekT = (a1,0,...,0)+···+(0,...,0,an) (1.1)
k=1 k=1
Bemerkung: Statt R als Grundk¨orper ist auch C m¨oglich.
Das Skalarprodukt ist dann definiert durch
(cid:88)n
uTv¯ = u v¯ .
i i
i=1
Problem: Gegeben sei A ∈ Rn,n, nichtsingul¨ar (d. h. invertierbar), b ∈ Rn.
Gesucht ist x ∈ Rn mit
Ax = b . (1.2)
x ∈ Rn existiert eindeutig fu¨r beliebige b ∈ Rn (bei exakter Arithmetik in Rn).
(Rechnen in endlicher Teilmenge von Rn mit inexakter Arithmetik → s.u.)
Spezialf¨alle:
1. A sei Diagonalmatrix, d.h. a = 0 fu¨r i (cid:54)= j, A := diag(a ).
ij ii
A nichtsingul¨ar ⇔ a (cid:54)= 0 fu¨r alle i =⇒
ii
x = b /a fu¨r alle i . (1.3)
i i ii
Also sind n Divisionen n¨otig zur Berechnung der L¨osung.
DieKomplexit¨atvon(1.3)wirdinElementaroperationengemessen,d.h.nurMul-
tiplikationen und Divisionen werden gez¨ahlt:
Diese Anzahl ist O(n) fu¨r n → ∞.
Die genaue Definition der Landausymbole erfolgt unten, hier reicht:
f(n) sei Polynom in n k–ten Grades, dann f(n) = O(nk)
(wesentlicher Wachstumsanteil fu¨r n → ∞).
(1.3) ist also ein “optimaler” Algorithmus, da n Unbekannte und n Eintr¨age in
A vorliegen, wobei a (cid:54)= 0 Eintrag in A heißt.
ij
5
2. A sei obere Dreiecksmatrix, d.h. a = 0 fu¨r i > j.
ij
Dann liegt ein gestaffeltes Gleichungssystem vor
a x + ··· ······ + a x = b
11 1 1n n 1
a x + ······ + a x = b
22 2 2n n 2
... ... ...
(1.4)
... ... ...
a x + a x = b
n−1,n−1 n−1 n−1,n n n−1
a x = b
n,n n n
und a (cid:54)= 0 fu¨r alle i, da fu¨r nichtsingul¨are A gilt
ii
0 (cid:54)= detA = a ···a .
11 nn
Ru¨ckw¨artssubstitution fu¨r (1.4):
x := b /a
n n nn
x := (b −a x )/a
n−1 n−1 n−1,n n n−1,n−1
.
.
. (1.5)
.
.
.
x := (b −a x −···−a x )/a
1 1 12 2 1n n 11
Es werden ben¨otigt:
n Divisionen,
(cid:88)n (cid:88)n−1
(n−1)n
(i−1) = i = Multiplikationen,
2
i=2 i=1
n(n−1+2)
d.h. Divisionen und Multiplikationen.
2
Die Komplexit¨at von (1.5) ist also O(n2) fu¨r n → ∞.
(n+1)n
Diesistwiederein“optimaler”Algorithmus,da Eintr¨ageinAvorliegen.
2
Algorithmus 1.1 (Ru¨ckw¨artssubstitution:)
for i=n:-1:1
x(i)=b(i);
for j=i+1:n
x(i)=x(i)-x(j)*A(i,j);
end
x(i)=x(i)/A(i,i);
end
(Die innere Schleife u¨ber j ist fu¨r i = n leer.)
6
3. Analog: A heißt untere Dreiecksmatrix, wenn a = 0 fu¨r i < j.
ij
Wieder gilt: A nichtsingul¨ar ⇔ a (cid:54)= 0 fu¨r alle i.
ii
x = A−1b wird bestimmt durch
Algorithmus 1.2 (Vorw¨artssubstitution:)
for i=1:n
x(i)=b(i);
for j=1:i-1
x(i)=x(i)-x(j)*A(i,j);
end
x(i)=x(i)/A(i,i);
end
Ist A eine nichtsingul¨are obere (untere) Dreiecksmatrix, dann ist (wie immer)
A−1 = (z1,...,zn), wobei Azi = ei
und zi, i = 1,...,n, ist mittels Algorithmus 1.1 (1.2) bestimmbar, insbesondere
zi = 1/a und zi = 0 fu¨r j = i+1,...,n (j = 1,...,i−1) .
i ii j
Es gilt also
Satz 1.3 Die Inverse einer nichtsingul¨aren oberen (unteren) Dreiecksmatrix A
ist eine obere (untere) Dreiecksmatrix A−1 und
(A−1) = 1/(A) fu¨r alle i = 1,...,n .
ii ii
4. Sei A ∈ Rn,n eine beliebige nichtsingul¨are Matrix:
Die Strategie des Gauß’schen Eliminationsverfahrens ist die Ru¨ckfu¨hrung auf
Spezialfall 2 durch sukzessive Elimination von Unbekannten:
Elimination von x aus den Gleichungen 2,...,n mittels Gleichung 1,
1
dann (wenn m¨oglich) Elimination von x aus den Gleichungen 3,...,n mittels
2
Gleichung 2 (modifiziert), usw..
Vorl¨aufige Zusatzvoraussetzung sei: a (cid:54)= 0.
11
Elimination von x aus Gleichung i, i = 2,...,n, erfolgt durch Addition von
1
−a /a ∗ Gleichung 1 zu Gleichung i. Das ist ¨aquivalent mit dem U¨bergang
i1 11
zum Gleichungssystem
A(2)x = b(2).
7
Genauer:
(cid:161) (cid:162) (cid:161) (cid:162)
A(1),b(1) A(2),b(2)
= =
   
a a ··· a b a a ··· a b
 11 12 1n 1   11 12 1n 1 
 a ... ... b   0 a(2) ··· a(2) b(2) 
 21 2   22 2n 2 
  −→  
 ... ... ... ...   ... ... ... ... ... 
   
a ··· ··· a b 0 a(2) ··· a(2) b(2)
n1 nn n n2 nn n
Dies ist darstellbar als:
(A(2),b(2)) := L(1)(A(1),b(1)), wobei A(1) := A, b(1) := b (1.6)
und  
1 0 ··· ··· 0
 
 a 
 21 
− 1 0 ··· 0
 
 a11 
L(1) :=  ... 0 ... ... ...  ∈ Rn,n . (1.7)
 
 ... ... ... ... 0 
 
 
a
n1
− 0 ··· 0 1
a
11
L(1) ist eine untere Dreiecksmatrix und l¨aßt sich darstellen als
L(1) = E −m(1)e1T ,
(cid:181) (cid:182)
a a T
wobei m(1) = 0, 21,..., n1 . (1.8)
a a
11 11
Dann gilt:
A(2)e1 = a e1,
11
d.h. in der ersten Spalte ist nur a(2) (cid:54)= 0 und
11
e1T(A(2),b(2)) = e1T(A(1),b(1)),
d.h. die erste Zeile bleibt unver¨andert.
8
Ist a(2) (cid:54)= 0 (weitere vorl¨aufige Voraussetzung), kann mit dem reduzierten Glei-
22
chungssystem fu¨r x ,...,x analog verfahren werden, um nun x aus den Glei-
2 n 2
chungen 3,...,n zu eliminieren. Das ist darstellbar als
(A(3),b(3)) := L(2)(A(2),b(2))
mit L(2) := E −m(2)e2T , (1.9)
(cid:195) (cid:33)
T
a(2) a(2)
wobei m(2) := 0,0, 32 ,···, n2
a(2) a(2)
22 22
denn:
L(2)A(2)e1 = L(2)a e1 = a (E −m(2)e2T)e1 = a e1 ,
11 11 11
d.h. die erste Spalte von A(2) bleibt unver¨andert, und
(cid:195) (cid:33)
e1TL(2)A˜(2) = e1T(E −m(2)e2T)A˜(2) = e1T −e1Tm(2)e2T A˜(2) = e1TA˜(2) .
(cid:124) (cid:123)(cid:122) (cid:125)
=0
(wobei A˜(2) = (A(2),b(2))), d.h. die erste Zeile von A˜(2) bleibt unver¨andert.
Allgemein gilt
Satz 1.4 Betrachte (1.2) mit nichtsingul¨arem A. Ist der folgende Algorithmus
(Eliminationsverfahren nach Gauß) durchfu¨hrbar, d.h. sind a(i) (cid:54)= 0 fu¨r alle i =
ii
1,...,n−1, dann formt er (1.2) in ein ¨aquivalentes Gleichungssystem mit oberer
Dreiecksmatrix um:
(cid:179) (cid:180)
A(1) := a(1) := A, b(1) := b, i := 1
ij
Fu¨r i = 1,...,n−1 :
Ist a(i) (cid:54)= 0, setze
ii
(cid:195) (cid:33)
a(i) a(i) T
m(i) := 0,...,0, i+1,i,..., n,i ,
a(i) a(i)
ii ii
L(i) := E −m(i)eiT, (1.10)
(cid:161) (cid:162) (cid:161) (cid:162)
A(i+1),b(i+1) := L(i) A(i),b(i) .
Ist a(i) = 0, breche ab: Der Algorithmus ist dann nicht durchfu¨hrbar.
ii
9
Beweis: Es genu¨gt, durch Induktion u¨ber k zu zeigen, daß die A(k) erfu¨llen:
eiTA(k)ej = 0 fu¨r 2 ≤ k ≤ n, 1 ≤ j < k, j < i ≤ n. (1.11)
Fu¨r k = 2 ist (1.11) schon gezeigt.
Es gelte (1.11) fu¨r k < n. Zu zeigen ist dann, daß (1.11) auch fu¨r k +1 gilt.
Seien i,j ∈ {1,...,n} mit 1 ≤ j < k +1, j < i ≤ n :
(cid:179) (cid:180)
eiTA(k+1)ej = eiTL(k)A(k)ej = eiT E −m(k)ekT A(k)ej
= eiTA(k)ej −eiTm(k)(ekTA(k)ej).
Hieraus folgt eiTA(k+1)ej = 0 nach Induktionsannahme fu¨r j < k, i > j.
Fu¨r j = k, i > j ist
eiTA(k)ek
eiTA(k+1)ej = eiTA(k)ek − ekTA(k)ek = 0
ekTA(k)ek
wegen
a(k) eiTA(k)ek
eiTm(k) = ik = .
a(k) ekTA(k)ek
kk
(cid:50)
Speichert man die Multiplikatoren, d.h. die Eintr¨age in den m(i) auf den jeweils
frei werdenden Pl¨atzen von A in der i-ten Spalte ab Zeile i + 1, so lautet der
Algorithmus
Algorithmus 1.5 (Gaußelimination ohne Pivotisierung)
for k=1:n-1
d=1/A(k,k);
for i=k+1:n
A(i,k)=A(i,k)*d;
b(i)=b(i)-A(i,k)*b(k);
for j=k+1:n
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
Das obige Vorgehen erfordert 1(n3−n)+ 1(n2−n) = O(n3/3) Multiplikationen
3 2
bzw. Divisionen (bei i.a. n2 Eintr¨agen in A).
10