Einführung in die Numerische Mathematik für Studierende der Fachrichtungen Informatik und Ingenieurwesen
Vorlesung Sommersemester 2010
Nicolas Neuÿ Institut für Angewandte und Numerische Mathematik Karlsruher Institut für Technologie
[email protected]
c
2008
Dieses Skript wird unter der GNU Free Documentation License, Version 1.3 zur Verfügung gestellt.
Diese Version wurde am 9. April 2011 erstellt.
Inhaltsverzeichnis 1
Einführung
1
1.1
Was ist Numerische Mathematik? . . . . . . . . . . . . . . . . . . . . . . . . .
1.2
Beispiel
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
Problemstellung
1.2.2
Modellierung
1.2.3
Mathematisches Modell
1.2.4
Numerisches Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.2.5
Theorie
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2.6
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2.7
Programm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2.8
Beobachtungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.2.9
Sinnvolle Problemstellungen . . . . . . . . . . . . . . . . . . . . . . . . .
4
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2
2.3
3
2 2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 6 7
Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.1
7
Beispiel
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gleitkomma-Arithmetik
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.1
Ganzzahl-Arithmetik . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.2
Wissenschaftliche Darstellung . . . . . . . . . . . . . . . . . . . . . . . .
8
2.2.3
IEEE-Gleitkommazahlen . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.4
Arithmetik
9
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
Kondition und Stabilität . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.3.1
Qualitative De nition der Kondition
. . . . . . . . . . . . . . . . . . . .
10
2.3.2
Absolute Kondition in 1D . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.3.3
Relative Kondition in 1D
11
2.3.4
Wiederholung: Mehrdimensionale Dierentiation
2.3.5
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
Absolute Kondition allgemeiner Funktionen
. . . . . . . . . . . . . . . .
13
2.3.6
Relative Kondition allgemeiner Funktionen
. . . . . . . . . . . . . . . .
14
2.3.7
Vektornormen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
2.3.8
Matrixnormen/Operatornormen . . . . . . . . . . . . . . . . . . . . . . .
15
2.3.9
Anwendung für Konditionsabschätzungen
. . . . . . . . . . . . . . . . .
16
. . . . . . . . . . . . . . . . . . . . . . . . .
17
Übungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2.3.10 Stabilität von Algorithmen 2.4
1
Ziel dieses Kurses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Grundlagen 2.1
1
1.2.1
1.2.10 Etwas Scilab 1.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
Lineare Gleichungssysteme
19
3.1
Problemde nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
3.2
Die Kondition der Matrixmultiplikation
20
i
. . . . . . . . . . . . . . . . . . . . . .
3.3
3.4
4
Abschätzung des absoluten Fehlers . . . . . . . . . . . . . . . . . . . . .
20
3.2.2
Abschätzung des relativen Fehlers
. . . . . . . . . . . . . . . . . . . . .
21
Die Kondition der Lösung linearer Gleichungssysteme . . . . . . . . . . . . . . .
21
3.3.1
Variation der rechten Seite . . . . . . . . . . . . . . . . . . . . . . . . . .
22
3.3.2
Variation der Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Das Gauÿsche Eliminationsverfahren
22
. . . . . . . . . . . . . . . . . . . . . . . .
23
3.4.1
Klassisches Vorgehen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.4.2
Motivation der (P)LR-Zerlegung
. . . . . . . . . . . . . . . . . . . . . .
24
3.4.3
Die (P)LR-Zerlegung anhand eines Beispiels . . . . . . . . . . . . . . . .
25
3.4.4
Verwendung der (P)LR-Zerlegung . . . . . . . . . . . . . . . . . . . . . .
26
3.4.5
Theorie der LR-Zerlegung ohne Zeilentausch . . . . . . . . . . . . . . . .
27
3.4.6
Die LR-Zerlegung für symmetrisch positiv de nite Matrizen . . . . . . .
27
3.4.7
Theorie der LR-Zerlegung mit Zeilentausch
. . . . . . . . . . . . . . . .
28
3.4.8
Implementation der LR-Zerlegung
. . . . . . . . . . . . . . . . . . . . .
28
3.4.9
Aufwandsbetrachtungen . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
3.4.10 Die Stabilität der LR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . .
31
Lineare Ausgleichsprobleme
32
4.1
Ein eindimensionales Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
4.2
Mehrdimensionale Probleme . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
4.2.1
Problem der kleinsten Fehlerquadrate
. . . . . . . . . . . . . . . . . . .
33
4.2.2
Ein mehrdimensionales Beispiel . . . . . . . . . . . . . . . . . . . . . . .
33
4.3 4.4
5
3.2.1
Lösung der Normalgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.3.1
. . . . . . . . . . . . . . . . . . . . . . . . . . .
34
Lösung mittels Orthogonaltransformationen . . . . . . . . . . . . . . . . . . . .
34
Kondition des Problems
4.4.1
Orthogonale Matrizen
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.4.2
QR-Zerlegung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.4.3
Anwendung der QR-Zerlegung
. . . . . . . . . . . . . . . . . . . . . . .
35
4.4.4
Anwendung auf Gleichungssysteme . . . . . . . . . . . . . . . . . . . . .
36
Nichtlineare Gleichungssysteme
37
5.1
37
5.2 5.3
Der eindimensionale Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1
Halbierungsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
5.1.2
Das Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
Der mehrdimensionale Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
5.2.1
40
Beispiel
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Theorie des Newton-Verfahrens
. . . . . . . . . . . . . . . . . . . . . . . . . . .
41
5.3.1
Kondition des Problems
. . . . . . . . . . . . . . . . . . . . . . . . . . .
41
5.3.2
Fixpunktiterationen
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
5.3.3
Kontraktionen und Fixpunktsatz
. . . . . . . . . . . . . . . . . . . . . .
44
5.3.4
Anwendung auf das Newton-Verfahren . . . . . . . . . . . . . . . . . . .
46
5.4
Das globale Newton-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
5.5
Iterative Lösung linearer Gleichungssysteme . . . . . . . . . . . . . . . . . . . .
48
5.5.1
Einführung
48
5.5.2
Konvergenztheorie
5.5.3
Die Beziehung zur Kondition
. . . . . . . . . . . . . . . . . . . . . . . .
50
5.5.4
Anwendung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ii
50
5.5.5 6
9
52 54
6.1
Einführung
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
6.2
Lagrange-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
6.2.1
Konstruktion
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
6.2.2
Gute Basiswahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
6.2.3
Interpolationsfehler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
6.2.4
Kondition der Polynominterpolation
58
. . . . . . . . . . . . . . . . . . . .
Stückweise polynomiale Approximationen (Splines) . . . . . . . . . . . . . . . .
60
6.3.1
Lineare Splines
60
6.3.2
Kubische Splines
6.3.3
Berechnung eines interpolierenden kubischen Splines
. . . . . . . . . . .
64
6.3.4
Fehlerabschätzung der kubischen Splines . . . . . . . . . . . . . . . . . .
65
6.3.5
Anwendung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
Numerische Quadratur
68
7.1
Einführung
68
7.2
Quadraturformeln
7.3
Konstruktion durch Polynominterpolation 7.3.1
8
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Interpolation
6.3
7
Optimale Dämpfung
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
68
. . . . . . . . . . . . . . . . . . . . .
69
Newton-Cotes-Formeln . . . . . . . . . . . . . . . . . . . . . . . . . . . .
70
7.4
Zusammengesetzte Formeln
7.5
Ausblicke
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
Numerische Lösung gewöhnlicher Dierentialgleichungen
75
8.1
Grundlagen
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
8.2
Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
8.2.1
Pendelgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
8.2.2
Räuber-Beute-Modell
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
76
8.3
Geometrische Interpretation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
8.4
Numerische Lösung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
8.4.1
Das explizite Euler-Verfahren . . . . . . . . . . . . . . . . . . . . . . . .
78
8.4.2
Verfahren höherer Ordnung
79
8.4.3
Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
8.4.4
Anwendung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81
8.4.5
Ausblicke
82
. . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Übersetzungstabelle
83
iii
1 Einführung 1.1 Was ist Numerische Mathematik? De nition: Numerische Mathematik beschäftigt sich mit der Entwicklung und Analyse von Methoden zur Lösung mathematischer Problemstellungen. Schema: Problem des praktischen Lebens
⇒ Mathematisches Modell & Messung von Parametern
⇒ . . . mit einem guten numerischen Verfahren . . . Brauchbare Lösung
1.2 Beispiel Wir wollen die Schwingung eines Stabpendels berechnen.
1.2.1 Problemstellung
• L:
Länge des Pendelstabs
• m:
Masse des Pendels
• g:
Erdbeschleunigung
• G = mg : • θ:
L
Gewichtskraft
Auslenkungswinkel
• FR = G sin θ:
G .
Rückstellkraft
FR
1
θ
1.2.2 Modellierung Wissen: Aus der Mechanik wissen wir
•
Auslenkungswinkel
•
Winkelgeschwindigkeit
•
Geschwindigkeit des Pendelgewichts ist
•
Beschleunigung des Pendelgewichts
•
Rückstellkraft
•
Andererseits:
θ
ist eine Funktion
˙ , ω(t) = θ(t)
FR (t) = −mg sin θ(t)
θ(t)
wobei
der Zeit
˙ := θ(t)
t
d θ(t). dt
v(t) = Lω(t)
a(t) = v(t) ˙ = Lω(t) ˙ = L¨θ(t) (Vorzeichen!)
FR (t) = ma(t) = mL¨θ(t)
(Newtonsches Kraftgesetz)
1.2.3 Mathematisches Modell Zusammensetzen liefert die Dierentialgleichung
¨θ(t) = − g sin θ(t) L zu deren Lösung man noch Anfangsbedingungen braucht:
θ(0) = θ0
˙ θ(0) = ω(0) = ω0 .
Bemerkung: Für diese Dierentialgleichung kann man die Lösung
θ(t)
nicht mittels ele-
mentarer Funktionen ausdrücken (was für viele Dierentialgleichungen der Fall ist). Man muss sie daher numerisch lösen.
1.2.4 Numerisches Verfahren
•
Approximiere
˙ θ(t)
und
ω(t) ˙
durch die Dierenzenquotienten (mit kleinem
˙ ≈ θ(t + h) − θ(t) , θ(t) h •
Andererseits wissen wir:
•
Man kann nun aus
˙ = ω(t) , θ(t)
ω(t) ˙ ≈
h > 0)
ω(t + h) − ω(t) . h
ω(t) ˙ =−
g sin θ(t) . L
(θ(t), ω(t)) eine Approximation zu (θ(t+h), ω(t+h)) berechnen: θ(t + h) ≈ θ(t) + hω(t) g ω(t + h) ≈ ω(t) + h(− sin θ(t)) L
2
1.2.5 Theorie Satz: (Analysis) Die Dierentialgleichung
¨θ(t) = − g sin θ(t) , L besitzt eine eindeutige Lösung
θ:R→R
˙ θ(0) = θ0 , θ(0) = ω0 für alle Anfangsbedingungen
(θ0 , ω0 ).
Satz: (Numerik) Das numerische Verfahren (es heiÿt explizites Euler-Verfahren)
θh,0 ωh,0
=
θ0 ω0
,
t>0
konvergiert, d.h. für jedes
lim
t N →∞, h= N →0
θh,k+1 ωh,k+1
=
θh,k ωh,k
ωh,k +h − Lg sin θh,k
für
k≥0
gilt
θh,N → θ(t) ,
lim
t N →∞, h= N →0
˙ . ωh,N → ω(t) = θ(t)
1.2.6 Implementation Eine für die Numerik sehr gut geeignete Programmiersprache ist das kommerzielle Matlab, oder aber freie Derivate davon, insbesondere Octave und Scilab. Vorteile:
•
Starke Sprache zur Matrizenmanipulation
•
Interaktives Arbeiten
•
Die Sprache ist recht konventionell, d.h. die meisten Sprachelemente kennt man aus anderen Computersprachen in sehr ähnlicher Form.
•
Viele Bibliotheken
Bemerkung: Wegen der freien Verfügbarkeit werden wir in diesem Kurs häu g Scilab benutzen. Da am KIT auch eine Matlab-Campuslizenz verfügbar ist, werden auch Beispiele in Matlab verwendet.
1.2.7 Programm Programm: (pendel.sci) g =
9.81;
function
L =
1.0;
[ theta , time ]
= pendel ( theta0 ,
h = T/N ; time ( 1 ) =
0.0;
3
omega0 ,
T,
N)
theta (1) = theta0 ; for
omega ( 1 ) = omega0 ;
i =1:N
t i m e ( i +1) = t i m e ( i )+h ; t h e t a ( i +1) = t h e t a ( i )+h ∗ omega ( i ) ; omega ( i +1) = omega ( i )−h ∗ g /L ∗ s i n ( t h e t a ( i ) ) ; end endfunction theta0 = T =
10.0;
0.0;
omega0 =
6.23;
N = 10000;
[ theta , time ]
= pendel ( theta0 ,
omega0 ,
T,
N) ;
p l o t ( time , t h e t a ) //
or :
p l o t ( t i m e , [ t h e t a , omega ] )
//
...
test
also
N=10000
and
values
of
omega0
around
6.26
1.2.8 Beobachtungen
•
Die Lösungskurve approximiert die korrekte Pendelbewegung für
•
Die Approximation ist leider nicht sehr genau (Fehler
•
Wenn man auch
ω
N → ∞.
O(h)).
plottet, sieht man den Austausch zwischen potentieller und
kinetischer Energie.
•
Das Verfahren ist nicht exakt energieerhaltend, sondern die Energie nimmt im Verlauf der Iteration zu.
1.2.9 Sinnvolle Problemstellungen Problem: Gegeben sei
L = 1m, θ0 = π , ω0 = 0 ms .
Zu dieser instabilen Lage betrachten
wir drei Fragestellungen: (I) Fällt das Pendel rechts herum, links herum, oder bleibt es oben stehen? (II) Wie ist die Lage des Pendels nach einer Sekunde? (III) Wie ist die Lage des Pendels nach einer Minute?
Frage: Welche der obigen Fragestellungen ist für die Praxis brauchbar? Antwort:
4
•
Frage (I) ist in der Praxis unbrauchbar, weil die Antwort Das Pendel bleibt für immer bei
θ=π
stehen durch beliebig kleine Störung der Anfangsdaten (welche
in der Praxis unvermeidbar sind, z.B. durch Luftzug oder Erschütterungen) falsch wird.
•
Frage (II) ist sinnvoll, weil die Antwort Das Pendel be ndet sich nach einer Sekunde bei
θ = π
auch bei kleinen Störungen der Anfangsdaten immer noch eine
nützliche Approximation ist.
•
Auch Frage (III) ist in der Praxis unbrauchbar. Eine genauere Betrachtung zeigt
√ ε mit der Zeit t um den Faktor et g/L 80 Faktor etwa 4 · 10 .
wächst. Für
θ=π
beobachten
nämlich, dass ein Anfangsfehler
L = 1m
und
t = 60s
ist dieser
Bemerkung: Falls man experimentell dennoch ein Stehenbleiben bei
sollte, liegt das wahrscheinlich an Reibungseekten in der Aufhängung, die im mathematischen Modell der Pendelgleichung nicht berücksichtigt werden. De nition: Man nennt Probleme/Fragestellungen, bei denen sich die Antwort stark ändert, obwohl die Eingangsdaten nur wenig variiert werden, schlecht gestellt oder schlecht konditioniert. Bemerkung: In der Praxis sind leider viele Probleme recht schlecht konditioniert. Manchmal ist eine geänderte Fragestellung die beste Abhilfe, ansonsten sollte aber der Numeriker danach trachten, den Fehler nicht auch noch durch ungeeignete Verfahren zu vergröÿern.
1.2.10 Etwas Scilab
•
zeros(m,n) Nullmatrix
•
A' Transposition
•
eye(m,n) Einheitsmatrix
•
det(A) Determinante
•
[1,2;3,4] Matrix aus Elementen
•
spec(A) Eigenwerte/Spektrum
•
A(i,j) Zugri auf Matrixelement
•
dot(x,y) Skalarprodukt
•
i:j Integerbereich
•
if condition . . . (Codezeilen). . . end
•
•
if condition . . . else . . . end
i:d:j Bereich mit Schrittweite d
•
for i=1:n . . . end Schleife
•
A*B Matrizenprodukt
•
function y = f (x) . . . endfunction
•
A\b Lösung LGS
Funktionsde nition
5
1.3 Ziel dieses Kurses
•
Sie lernen die Arbeitsweise
•
von numerischen Verfahren kennen,
•
die für viele wissenschaftliche Anwendungen wichtig sind.
•
Dieses Hintergrundwissen ist notwendig,
•
um Standardsoftware
•
korrekt und sicher anwenden zu können.
6
2 Grundlagen 2.1 Motivation 2.1.1 Beispiel
Aufgabe: Man berechne
ex = exp(x) =
∞ X xk k=0
Algorithmus: Berechne
N X k=0
xk k!
wobei
N
k!
auf eine Genauigkeit
so gewählt ist, dass
Satz: Bei exakter Arithmetik berechnet dieser Algorithmus
|
ε > 0.
x 1 |≤ N 2
und
|
exp(x) auf einen Fehler ε > 0
genau. Beweis:
∞ ∞ l ∞ X X xk xN X xl 1 |≤| || |≤ε | =ε k! N ! (N + 1) · · · (N + l) 2 l=1 l=1 k=N +1
Programm: (exp_series.sce) e p s = 1 . 0 e −15; function
summe = e x p _ s e r i e s
(x)
summe = 1 ; k = 1; summand = x / k ; while
( a b s ( summand )>e p s )
|
( abs ( x/k ) >0.5)
summe = summe + summand ; k = k +1; summand = summand ∗ x / k ; end endfunction for
x =40: − 10: − 40
d i s p ( [ x , exp ( x ) , e x p _ s e r i e s ( x ) ] ) ; end
Beobachtung: Wir erhalten folgende Tabelle:
7
do
xN | ≤ ε. N!
x
exp(x)
exp_series(x)
40
2.354E+17
2.354E+17
30
1.069E+13
1.069E+13
20
4.852E+08
4.852E+08
10
22026.466
22026.466
0
1.0000000
1.0000000
-10
0.0000454
0.0000454
-20
2.061E-09
5.622E-09
-30
9.358E-14
-0.0000307
-40
4.248E-18
-3.1657319
x 0. 1 . exp_series(−x)
Oenbar versagt unser Algorithmus für Abhilfe: Berechne für
x 0.
für
x ∈ V , λ ∈ R.
kx + yk ≤ kxk + kyk
für alle
x, y ∈ V .
Beispiele:
•
Auf
Rn
für
1≤p n Sinn macht, und auch, dass wir mit k·k2 die Euklidische Norm n m im R und R bezeichnen). Beweis: Siehe [Golub-van-Loan]. Problem: Wenn man
x
über die Normalgleichung bestimmt, ist die Fortp anzung des
κ2 (At A) = κ2 (A)2 gegeben (wegen der Multiplikation At b muss 3 man insgesamt sogar mit κ2 (A) rechnen). Für viele Anwendungsprobleme ist aber schon 2 κ2 (A) groÿ und κ2 (A) oder gar κ2 (A)3 inakzeptabel! relativen Fehlers durch
4.4 Lösung mittels Orthogonaltransformationen 4.4.1 Orthogonale Matrizen
Q ∈ Rn×n (bzw. die zugehörige lineare Abbildung) heiÿt orthogonal, n X Euklidische Skalarprodukt hx, yi = xi yi erhält, d.h.
De nition: Eine Matrix wenn sie das
i=1
hQx, Qyi = hx, yi ∀x, y ∈ Rn .
Interpretation:
Q
ist längen- und winkelerhaltend.
34
Bemerkung: Äquivalent zur obigen Bedingung ist die Forderung −1 anders geschrieben Q = Qt .
QQt = Qt Q = I
oder
Beispiele: Drehungen, Spiegelungen und Kombinationen davon.
4.4.2 QR-Zerlegung m×n De nition: Eine QR-Zerlegung einer Matrix A ∈ R hat die Form A = QR, wobei m×m m×n Q∈R eine orthogonale Matrix und R ∈ R eine obere Dreiecksmatrix ist (d.h.
Rij = 0
für
i > j ).
Satz: Zu jedem
A ∈ Rm×n
gibt es eine QR-Zerlegung.
Beweis: Man kann die QR-Zerlegung auf verschiedene Weise konstruieren, etwa durch Spiegelungen, Drehungen oder auch die Gram-Schmidt-Orthogonalisierung. Bemerkung:
•
Der Scilab/Matlab-Befehl zur Berechnung einer QR-Zerlegung heiÿt qr.
•
In manchen Fällen (z.B. wenn
A ∈ Rn×n invertierbar ist und wenn Rii > 0 gefordert
wird), ist die QR-Zerlegung sogar eindeutig (siehe Übung).
4.4.3 Anwendung der QR-Zerlegung Beobachtung: Weil
Q
orthogonal ist und
m≥n
ist, gilt mit
˜b = Qt b
kAx − bk22 = kQRx − bk22 = kRx − Qt bk22 R1:n,1:n x − ˜b1:n 2 = k k2 −˜bn+1:m = kR1:n,1:n x1:n − ˜b1:n k2 + k˜bn+1:m k2 2
Falls
Rang(A) = n,
so ist auch
Rang(R) = n,
2
und das Minimum erhält man oenbar
durch Lösung des Dreieckssystems
R1:n,1:n x1:n = ˜b1:n .
Notation: Wir haben die Matlab-Notation zur Bezeichnung von Untermatrizen durch Indexbereiche verwendet. Bemerkung: Man beachte, dass man zur Lösung von
R1:n,1:n x1:n = ˜b1:n = (Q1:m,1:n )t b1:m
35
nur Teile der QR-Zerlegung braucht, nämlich
Q1:m,1:n
und
R1:n,1:n .
Es gibt daher den
Scilab-Befehl qr(A,"e") (Matlab: qr(A,0)), der genau diese Matrizen berechnet. Programm: Die Lösung des Kleinste-Quadrate-Problems mit Matrix b ∈ Rm erhält man daher in Scilab auch wie folgt: [ Q, R ]
A ∈ Rm×n
und
= q r (A, " e " )
x = R \
(Q' ∗ b )
Beispiel: Die Anwendung auf das Motorenbeispiel zeigt, dass man dieselbe Lösung erhält wie vorher.
4.4.4 Anwendung auf Gleichungssysteme
m = n entspricht das Kleinste-Fehlerquadrate-Problem der Lösung des Gleichungssystems Ax = b. Die QR-Zerlegung ist daher auch für die Lösung eines t Gleichungssystems Ax = b geeignet (man löst einfach das Dreieckssystem Rx = Q b). Beobachtung: Im Fall
Satz: Im Gegensatz zur LR-Zerlegung ist die Stabilität der QR-Zerlegung in der Spektralnorm garantiert: es gilt Beweis: Da
Q
κ2 (Q) = 1
und
κ2 (R) = κ2 (A).
längenerhaltend ist, folgt
κ2 (Q) = 1.
Wegen
kAxk2 = kQRxk2 = kRxk2 folgt
κ2 (A) = κ2 (R)
(siehe Übungsaufgabe).
Aber: Der Rechenaufwand zur Berechnung einer QR-Zerlegung ist 2-3 mal höher als bei der LR-Zerlegung. Man verwendet die QR-Zerlegung daher meist nur bei schlecht konditionierten Problemen.
Übungen Aufgabe: Sei orthogonaler
A ∈ Rn×n invertierbar. Zeigen Sie, dass die Zerlegung A = QR mit Matrix Q und rechter oberer Dreiecksmatrix R mit positiver Diagonale
eindeutig ist.
36
5 Nichtlineare Gleichungssysteme Problem: Gegeben
F : Rn → Rn
und
y ∈ Rn .
Suche
x ∈ Rn
mit
F (x) = y .
Bemerkung: Ohne wesentliche Einschränkung könnte man sich mittels den Fall
y=0
F ← F − y auf y recht
beschränken. In Theorie und Praxis ist allerdings ein allgemeines
praktisch. Anwendungen: In dieser Allgemeinheit fällt fast alles darunter. Oft kommt das Gleichungssystem n aus einem Optimierungsproblem, weil ein Extremum einer Funktion Φ : R → R das System
F (x) := ∇Φ(x) = 0
lösen muss.
5.1 Der eindimensionale Fall Problem: Sei
f :R→R
und
y ∈ R.
Suche ein
x∈R
mit
f (x) = y .
Beispiele:
√ 2
⇒
x∈R
•
Berechne
•
Wann lag der Ölpreis bei 70 Dollar/Barrel?
numerisch
Finde
mit
x2 = 2. ⇒
Finde Zeit mit Ölpreis(Zeit)
= 70.
5.1.1 Halbierungsverfahren Satz: (Zwischenwertsatz) Sei
f (b)
oder
f (b) ≤ y ≤ f (a).
f : R → R stetig und y, a, b ∈ R mit a ≤ b und f (a) ≤ y ≤ x ∈ [a, b] mit f (x) = y .
Dann gibt es
Beweis: Analysis ([HM-Skript], Satz 7.11). Frage: Kann man diesen Satz zur Lösung von
f (x) = y
verwenden?
a+b und berechne f (c). Je nach dem Vorzeichen 2 kann man eines der Teilintervalle [a, c] oder [b, c] so auswählen, dass darin ein
Idee: (Halbierungsverfahren) Nimm
c=
y−f (c) x mit f (x) = c liegt. Dann kann man sich aber auf dieses kleinere Intervall beschränken, und mehrfache Intervallhalbierung bestimmt eine Lösung x beliebig genau. von
Programm: function while
s o l u t i o n = i n t e r v a l _ s o l v e ( f , y , eps , a , b ) a b s ( b−a )>e p s
c = ( a+b ) / 2 ; disp ( c ) ;
37
if
( f ( a )−y )
∗ ( f ( c )−y )e p s
( f , Df , y , x0 , e p s )
disp (x) x=x+r e s / Df ( x ) ; r e s=y− f ( x ) ; end endfunction
√ Anwendungsbeispiel: Berechnung von function
2:
y = square (x)
y=x ∗ x ; endfunction function
y = twice (x ) ;
y=2∗ x ; endfunction newton
( square ,
twice ,
2,
1,
1 e − 12)
Beobachtungen:
•
√ Für positive Startwerte konvergiert das Verfahren gegen
√
− 2. •
Für den Startwert
0
2,
für negative gegen
ist es nicht de niert.
Die Konvergenz ist sehr schnell, wenn der Startwert nahe bei einer der Lösungen
√ ± 2
ist.
Anwendungsbeispiel: Berechnung der Nullstelle function
y = f (x)
y=t a n h ( x ) endfunction
39
x=0
von
f (x) = tanh(x):
function
y = Df ( x )
y=1− t a n h ( x ) ^ 2 endfunction newton
(f ,
Df ,
0,
1 e − 12)
1,
Beobachtung: Hier konvergiert das Verfahren nur, wenn der Startwert nahe bei ist (|x|
≤ 1
reicht). Für groÿe Startwerte (z.B.
x = 2)
0
werden die Iterierten mit
wechselndem Vorzeichen sehr schnell gröÿer.
5.2 Der mehrdimensionale Fall Problem: Zu
F ∈ C 1 (Rn , Rn )
und
Erinnerung: Die Ableitung von stetige matrixwertige Funktion,
y ∈ Rn
wollen wir das Problem
F ∈ C 1 (Rn , Rn ) bezeichnen 0 n n×n also DF ∈ C (R , R ).
F (x) = y
wir mit
DF .
lösen. Es ist eine
5.2.1 Beispiel Aufgabe: Schneide den Kreis
x21 + x22 = 1 mit der Kurve
x 2 = e x1 . Zusammengefasst führt das zum Gleichungssystem
1 F (x) = 0
wobei
2 x1 x1 + x22 F : 7→ . x2 x2 − ex1
Bemerkung: Eine Gra k zeigt, dass es zwei Lösungen gibt: die eine ist der Punkt
die andere liegt in der Nähe von
0 , 1
−1 . 1/e
Idee: Nach De nition der Ableitung gilt wieder die lineare Approximation
F (x) ≈ F (x0 ) + DF (x0 )(x − x0 ) . Daher sollte sich das Newton-Verfahren einfach übertragen lassen. Algorithmus: (Verallgemeinertes Newton-Verfahren für n-dimensionale Probleme) Gegeben 1 n n 0 n n×n seien F ∈ C (R , R ), die Ableitung DF ∈ C (R , R ) und y ∈ Rn , sowie ein n Startwert x0 ∈ R . Dann erhalten wir die Folge der Newton-Iterierten durch
xk+1 = xk + (DF (xk ))−1 (y − F (xk )) .
40
Bemerkung: Die Ableitung
DF (xk )
ist eine Matrix (die sogenannte Jacobi-Matrix) und
muss invertiert werden. Das Verfahren bricht zusammen, wenn es auf einen Punkt mit singulärer Jacobi-Matrix läuft. Programm: (Newton-Verfahren) function
x = newton
( f , Df , y , x0 , e p s )
x=x0 ; r e s=y− f ( x ) ; while
norm ( r e s )>e p s
x=x+(Df ( x ) )
\
res ;
disp (x) ; r e s=y− f ( x ) ; end endfunction
Bemerkung: Es waren nur zwei Änderungen des 1D-Programms nötig: erstens die Verwendung von norm anstelle von abs, zweitens die Verwendung des Backslash-Operators. Das entstandene Programm funktioniert für beliebige Dimensionen. Beispiel: Den zweiten Schnittpunkt des Einheitskreises mit dem Graph
x2 = ex1
berechnet
der Code function
y = f (x)
y =[ x (1)^2+ x ( 2 ) ^ 2 ;
x (2) − e x p ( x ( 1 ) ) ] ;
endfunction f u n c t i o n D = Df ( x ) D=[ 2 ∗ x ( 1 ) , 2 ∗ x ( 2 ) ;
−e x p ( x ( 1 ) )
,1];
endfunction newton
(f ,
Df ,
Man erhält die Lösung
−1;0] , −0.91656258 . 0.39989127
[1;0] ,
x≈
[
1 e − 12)
5.3 Theorie des Newton-Verfahrens 5.3.1 Kondition des Problems
F ∈ C 1 (Rn , Rn ) und y = F (x). Wenn DF (x) invertierbar n n ist, so existieren Umgebungen U ⊂ R von x und V ⊂ R von y sowie eine Abbildung F −1 : V → U , so dass F −1 ◦ F |U = IU und F ◦ F −1 |V = IV . Auÿerdem gilt D(F −1 )(y) = (DF (x))−1 .
Satz: (Umkehrfunktion) Sei
Beweis: Analysis ([HM-Skript, Satz 9.9 und Satz 22.13]).
41
∆x der Lösung x der nichtlinearen Gleichung F (x) = y y um ∆y ?
Frage: Wie groÿ ist die Änderung bei Störungen der rechten Seite Antwort: Wenn
DF (x)
regulär ist, so existiert lokal
F −1
und es gilt nach dem Kapitel
Grundlagen
˙ kD(F −1 )(y)kk∆yk = k(DF (x))−1 kk∆yk . k∆xk ≤ Wieder bezeichnet n×n auf R .
k·k eine beliebige Vektornorm im Rn sowie die zugehörige Operatornorm
Bezeichnung: Wir nennen die Gröÿe
κabs (F (x) = y) = kD(F −1 )(y)k = k(DF (x))−1 k . die Kondition der Lösung
x der Gleichung F (x) = y . Sie hängt von der gewählten Norm
ab. Bemerkungen:
•
Für nichtlineare Probleme kann es zu einem y mehrere Lösungen geben (z.B. hat √ √ x2 = 2 die Lösungsmenge {+ 2, − 2}). Die Kondition einer bestimmten Lösung
x
quanti ziert das lokale Verhalten einer Umkehrfunktion, die auf eine Umgebung
von
•
x
abbildet!
Man beachte, dass wir für die Fehlerfortp anzung bei nichtlinearen Problemen absolute Fehler betrachten, und dementsprechend auch die absolute Kondition verwenden. Der Grund ist einfach: Im allgemeinen gilt nicht
F (0) = 0,
so dass der
Begri des relativen Fehler nicht überall Sinn macht.
Beispiel: Wir betrachten wieder das Problem des Schnitts von Einheitskreis und Graph der Exponentialfunktion. Die Ableitung hatte hier die Form
Die Konditionen der
2x1 2x2 DF (x) = −ex1 1 0 −0.91656258 . . . Schnittpunkte und 1 0.39989127 . . .
berechnen sich dann
bezüglich der Maximumsnorm als
−1 1 3 −1 0 −1 0 2 2 k∞ = • kDF ( ) k∞ = k k∞ = k 1 0 1 −1 1 2 2 −0.91656258 . . . −1 • kDF ( ) k∞ = 1.4755949 . . . 0.39989127 . . . Für beide Schnittpunkte ist der Wert moderat, die Schnittpunkte sind also gut konditioniert.
42
Beispiel: Lineare Probleme sind ein Spezialfall von nichtlinearen. Wir wollen den Schnittpunkt der Geraden
y = kx + m
Die Lösung ist
k
mit der x-Achse bestimmen. Dazu suchen wir die Lösung von
0 x y F (ξ) = , F : 7→ . m y y − kx m 0 1 −k , die Ableitung DF = . natürlich ξ = 0 −k 1
Für kleines
ist diese Matrix fast singulär und die Norm der Inversen wird beliebig groÿ. Der
Schnittpunkt ist also schlecht konditioniert. Geometrisch bedeutet dies, dass die Lage des Schnittpunkts beim schleifenden Schnitt sehr schlecht bestimmt ist, was ja der Erfahrung entspricht. Beispiel: Wir betrachten das Problem
F (x) = y ,
mit
Es hat die eindeutig bestimmte Lösung von
x
F (x) = x2
x = 0,
und
y = 0.
allerdings ist
F 0 (x) = 0,
die Kondition
also nicht de niert (oder unendlich groÿ). Solche Probleme nennt man schlecht
gestellt, und das Verhalten der Lösungsmenge bei Variation von
y
kann beliebig bösartig
sein. Im obigen Fall beobachten wir folgendes:
•
√ ∆y = ε > 0 gibt es zwei Lösungen x = ± ε. Es gilt |∆x| = √1ε |∆y|, kleines ε eine beliebig groÿe Verstärkung des absoluten Fehlers bedeutet.
•
Für
Für
∆y = ε < 0
was für
gibt es gar keine reelle Lösung.
Bemerkung: Wenn man in der Praxis auf ein solch schlecht gestelltes Problem stöÿt, sollte man seine Herkunft zurückverfolgen. Oft wird man feststellen, dass es durch einen Fehler zustande kam, und der Anwender eigentlich die Lösung eines ganz anderen Problems wünscht.
5.3.2 Fixpunktiterationen De nition: Wir betrachten allgemeine Iterationen der Form
x0 ∈ Rn , Wenn die Folge der Grenzwert
(xk )k∈N
x∗
gegen
xk+1 = Φ(xk )
x∗ ∈ Rn
für
k ≥ 1. Φ : Rn → Rn stetig ist, so ist gilt Φ(x∗ ) = x∗ . Wir nennen die
konvergiert und
oenbar ein Fixpunkt von
Φ,
d.h. es
Iteration daher auch eine Fixpunktiteration. Beobachtung: Die Newton-Iteration
xk+1 = xk + (DF (xk ))−1 (y − F (xk )) kann man als spezielle Fixpunkt-Iteration mit Iterationsfunktion
Φ(xk ) := xk + (DF (xk ))−1 (y − F (xk )) auassen.
43
5.3.3 Kontraktionen und Fixpunktsatz De nition: Sei
U ⊂ Rn
und
Φ:U →U
θ0
gewählt werden, damit
F (x) = x2
auf dem Intervall
[−α, α]
eine Kontraktion ist? Antwort: Wir berechnen α < 12 gewählt wird.
F 0 (x) = 2x
und sehen, dass
|F 0 (x)| ≤ θ < 1
nur gilt, wenn
Bemerkung: Wir sehen hier auch, dass die Forderung nach einer Kontraktion nur hinreichend, aber nicht notwendig für die Existenz eines eindeutig bestimmten Fixpunkts ist. Denn oensichtlich hat
F |[−α,α]
ja für alle
α 0, so dass das Newton-Verfahren n in der Kugel Bδ (x∗ ) := {x ∈ R : kx − x∗ k < δ} mit Rate θ konvergiert. Beweisskizze: Der Einfachheit halber zeigen wir es nur für den 1D-Fall und ohne die Konstanten zu quanti zieren (und bezeichnen daher wieder F mit f ). 0 Weil f zweimal stetig dierenzierbar ist, und f (x∗ ) 6= 0, ist Φ in einer Umgebung von
x∗
stetig dierenzierbar mit Ableitung
Φ0 (x) = 1 − f 00 (x)f 0 (x)−2 (y − f (x)) − f 0 (x)−1 f 0 (x) = −f 00 (x)f 0 (x)−2 (y − f (x)) . y = f (x∗ ) gilt Φ0 (x∗ ) = 0. Und weil Φ0 stetig ist, gibt es daher zu einer beliebigen 0 Kontraktionszahl θ ein δ > 0, so dass |Φ (x)| ≤ θ für x ∈ Bδ (x∗ ). Daher ist Φ eine Kontraktion mit Rate θ auf U = Bδ (x∗ ). Mit unserer Version des Fixpunktsatzes folgt Wegen
dann auch die Konvergenz der Fixpunktiteration. De nition: Eine Folge es ein
C>0
(xk ) im Rn
heiÿt quadratisch konvergent gegen ein
x∗ ∈ Rn , wenn
gibt mit
kxk+1 − x∗ k ≤ Ckxk − x∗ k2 . Satz: Unter den obigen Voraussetzungen ist das Newton-Verfahren quadratisch konvergent.
Beweisskizze: In den oben konstruierten Umgebungen von
|y − f (x)| ≤ C1 |x − x∗ | , Die Kontraktionsrate in der Kugel um durch
Ckx − x∗ k
mit
x∗ gelten auch folgende Abschätzungen:
|f 0 (x)−1 | ≤ C2 ,
|f 00 (x)| ≤ C3 .
x∗ mit Radius kx−x∗ k lässt sich daher abschätzen
C = C1 C2 C3 .
Bemerkung: Diese quadratische Konvergenz in der Nähe der Lösung tatsächlich in der Praxis, vorausgesetzt: 1.
F
und
DF
2. die Lösung 3.
F
und
DF
sind ausreichend glatt (F
x∗
∈ C 2 , DF ∈ C 1 ),
ist gut konditioniert (DF (x∗ ) ist regulär)
werden ausreichend genau berechnet und
4. es wurden keine Fehler in der Implementation gemacht.
46
x∗
beobachtet man
5.4 Das globale Newton-Verfahren Problem: In der Praxis ist es oft schwierig, Startwerte zu nden, die so nahe bei einer Lösung sind, dass das unmodi zierte Newton-Verfahren konvergiert. Wir haben am Beispiel des
tanh
gesehen, dass dies sogar schon in sehr einfachen Situationen auftreten
konnte. Beobachtung: Wenn
F
dierenzierbar ist und
DF
regulär ist, so ist die Newton-Richtung
v = DF −1 (xk )(y − F (xk )) auf jeden Fall eine Richtung, in der die Funktion sieht man, weil ja wegen der Dierenzierbarkeit
s 7→ ky − F (xk + s)k kleiner wird. Dies in xk die Approximation
. F (x) = F (xk ) + DF (xk )(x − xk ) gilt. Man kann aber nicht garantieren, dass man den Schritt in voller Länge durchführen kann, ohne den Gültigkeitsbereich dieser Approximation zu verlassen! Idee: Dämpfe den Newton-Schritt mit einem Faktor
vk = DF −1 (xk )(y − F (xk )) , wobei
ωk
ωk ∈]0, 1]:
xk+1 = xk + ωk vk
geeignet gewählt wird.
Bezeichnung: Das entstehende Verfahren heiÿt Newton-Verfahren mit Dämpfungsstrategie oder globales Newton-Verfahren. Bemerkung: Es gibt verschiedene Strategien, den Dämpfungsfaktor zu wählen. Eine oft angewandte Strategie funktioniert im wesentlichen wie folgt: 1. Ein gegebener Dämpfungsfaktor
ωk
wird so lange halbiert, bis er zulässig ist.
Dabei bedeutet zulässig, dass die Verringerung der Fehlernorm im wesentlichen dem erwarteten linearen Abfall entspricht:
ky − F (xk + ωk vk )k ≤ (1 − θωk )ky − F (xk )k . Hierbei ist
θ ∈ (0, 1)
beliebig, eine in der Praxis bewährte Wahl ist
2. Am Anfang wählt man z.B.
ω0 = 1
θ=
1 . 2
(oder aber etwas kleineres, wenn das Problem
schwierig ist). Später führt man obige Prozedur ausgehend vom Startwert
min(2ωk , 1)
ωk+1 =
aus durch. Durch diese anfängliche Verdopplung wird die Dämpfung
allmählich wieder ausgeschaltet, sobald das Verfahren mit gröÿeren Schritten gut funktioniert. Hierdurch man erhält man insbesondere nahe der Lösung die quadratische Konvergenz des ungedämpften Verfahrens zurück. Programm: (Globales Newton-Verfahren)
47
function
x = global_newton
( f , Df , y , x0 , e p s )
x = x0 ; r e s = y− f ( x ) ; omega = 1 ; while
norm ( r e s )>e p s
v = Df ( x )
\
do
res ;
omega = min ( 2 ∗ omega , 1 ) ; while
norm ( y− f ( x+omega ∗ v ) ) > (1 − t h e t a ∗ omega ) ∗ norm ( r e s )
do
omega=omega / 2 ; end x = x+omega ∗ v ; d i s p ( omega ) ; disp (x) ; r e s=y− f ( x ) ; end endfunction function
y = f (x)
y = tanh ( x ) ; endfunction function
y = Df ( x )
y = 1− t a n h ( x ) ^ 2 ; endfunction theta =
0.5;
g l o b a l _ n e w t o n ( f , Df , 0 , 1 0 , 1 e − 10)
Beobachtung: Diese Newton-Variante konvergiert bei der Berechnung der Nullstelle des
tanh
für alle Startwerte
x0 ∈ R
gegen die Lösung
x = 0.
5.5 Iterative Lösung linearer Gleichungssysteme 5.5.1 Einführung Beobachtung: Wenn man das Newton-Verfahren auf ein lineares Gleichungssystem
F (x) := Ax = y , mit invertierbarem
A
A ∈ Rn×n ,
x, y ∈ Rn
anwendet, so ergibt sich der Newton-Schritt
xn+1 = xn + A−1 (y − Axn ) = A−1 y . Das Verfahren liefert also bereits nach einem Schritt die exakte Lösung.
48
Beobachtung: Wenn man allerdings das Gleichungssystem im Newton-Schritt nur approximativ −1 löst, d.h. anstelle von A nur eine Approximation B (den sogenannten Vorkonditionierer) anwendet, so liefert die Vorschrift
xn+1 = xn + B(y − Axn ) ein iteratives Verfahren zur Lösung des linearen Gleichungssystems
Ax = y .
Bemerkung: Eine solche iterative Lösung kann günstig sein, weil insbesondere für hochdimensionale −1 Systeme (d.h. groÿes n) die Berechnung von A oder die Lösung des Gleichungssystems 3 in O(n ) Operationen zu aufwendig ist. Beispiel: Wir betrachten die Matrix
I =: B .
A = I +S
mit
kSk ≤
1 und setzen einfach 2
A−1 ≈
Damit erhalten wir die Iteration
xn+1 = Φ(xn ) := xn + B(y − Axn ) = y − Sxn . Die Abbildung
Φ
ist aber oenbar wegen
kΦ(ξ1 ) − Φ(ξ2 )k = kS(ξ1 − ξ2 )k ≤ kSkkξ1 − ξ2 k eine Kontraktion mit Kontraktionsfaktor Bezeichnung:
1 . 2
k·k bezeichnet hier wieder sowohl eine (beliebige) Vektornorm als auch die
zugehörige Matrixnorm.
A wie oben und ε = 2−k eine vorgegebene Genauigkeit. Dann approximiert die obige Iteration mit Startwert x0 = 0 die Lösung x von Ax = y in höchstens k Schritten bis auf einen Fehler der Gröÿe kxk − xk ≤ εkxk. Folgerung: Sei
Bemerkungen: Der Aufwand für jeden einzelnen Iterationsschritt n×n stark von S ∈ R ab:
•
Für allgemeines 2 also n AO.
•
Für dünnbesetzte Matrizen
hängt
S ist es im wesentlichen der Aufwand für eine Matrix-Vektor-Multiplikation, S
(d.h. nur
Null verschieden) ist der Aufwand
•
xn+1 = y − Sxn
N
N n2
Matrixeinträge von
S
sind von
AO.
7→ y − Sx schnell ∈ Rn . Obwohl S 2n AO, um y − Sx =
Manchmal gibt es eine besondere Technik, um die Abbildung x t durchzuführen. Ein Beispiel wäre S = vw mit Vektoren v, w im allgemeinen voll besetzt ist, braucht man hier nur etwa y − v(wt x) zu berechnen.
49
5.5.2 Konvergenztheorie Satz: Wir betrachten die Iteration
x0 ∈ R n ,
xn+1 = Φ(xn ) := xn + B(y − Axn ) = By + (I − BA)xn
Angenommen, es gibt eine Norm
k·k
(Vektor- und zugehörige Matrixnorm) mit
kI − BAk = ρ < 1 , so konvergiert die Iteration mit der Konvergenzrate Beweis: Oenbar ist
Φ
ρ
bezüglich dieser Norm.
dann eine Kontraktion mit der Rate
Interpretation: Gute Konvergenz erhält man oenbar für bedeutet, dass
B
nahe bei der Inversen von
A
ρ
bezüglich
k·k.
ρ = kI − BAk 1.
Dies
liegt!
5.5.3 Die Beziehung zur Kondition Satz: Es gelte
kI − BAk ≤ ρ < 1 . Dann gilt
κ(BA) ≤ wobei
κ(A) = kAkkA−1 k
1+ρ 1−ρ
die Kondition der Matrix
ρ≥
A
bezeichnet. Umgekehrt gilt
κ(BA) − 1 . κ(BA) + 1
Interpretation: Gute Konvergenz der linearen Iteration ist nur dann möglich, wenn die
BA moderat ist. Dies motiviert auch die Bezeichnung Vorkonditionierer Wahl von B führt nämlich zu κ(BA) κ(A).
Kondition der Matrix für
B:
eine gute
Beweis: Wir setzen
S = BA
und gehen aus von
kI − Sk = ρ < 1.
k(I − S)xk = kx − Sxk ≤ ρkxk ,
Dies bedeutet
x ∈ Rn .
Die Dreiecksungleichung liefert nun einerseits
kxk − kSxk ≤ ρkxk
⇔
kSxk ≥ (1 − ρ)kxk ,
kSxk − kxk ≤ ρkxk
⇔
kSxk ≤ (1 + ρ)kxk ,
andererseits
Nun folgt
κ := κ(S) =
maxkxk=1 kSxk 1+ρ ≤ minkxk=1 kSxk 1−ρ
50
und Au ösen dieser Ungleichung nach
Bemerkung: Wenn
κ 1,
so gilt
ρ
liefert auch
1− κ−1 = κ+1 1+
1 κ 1 κ
ρ≥
≈1−
κ−1 . κ+1
2 . κ
Bemerkung: Die Kunst bei der iterativen Lösung linearer Gleichungssysteme liegt vor allem in der guten Wahl des Vorkonditionierers
B.
1.
B=
1 I für symmetrisch positiv de nites kAk2
2.
B=
1 At für beliebiges kAt Ak2
Für einfache Optionen wie
A
oder
A
konvergiert die Iteration zwar (siehe Übung), aber meist sehr langsam.
5.5.4 Anwendung Beispiel: Sei
A ∈ Rn×n
gegeben als
2
−1 An =
Es gilt
kAn k2 ≈ 4
−1 ..
.
..
.
..
.
. .. . −1 −1 2
(siehe Übung).
Ziel: Obwohl man ein LGS mit der obigen Systemmatrix sehr e zient direkt lösen kann 1 (wie?), wollen wir hier die iterative Lösung betrachten. Wir wählen daher B = I und 4 untersuchen, wie schnell die folgende Iteration konvergiert:
1 xk+1 = xk + (y − Axk ) . 4 Programm: (Iterative Lösung von function n =
y = apply_A
An x = (1, . . . , 1)t )
(x)
size (x , 1 ) ;
y = 2∗ x
−
[x(2:n );0]
−
[ 0 ; x (1: n−1)];
endfunction function
[ x , count ]
= iter_solve
(y)
x = zeros ( size (y , 1 ) , 1 ) ; r e s = y−apply_A ( x ) ; count = 0 ; while
( norm ( r e s )>1 e − 8)
do
x = x +0.25∗ r e s ;
51
r e s = y−apply_A ( x ) ; c o u n t = c o u n t +1; end endfunction Bemerkung: Man beachte, dass die Matrix
A
nicht einmal abgespeichert werden muss!
Zur Durchführung des Verfahrens reicht es vollkommen aus, dass man die Abbildung
x 7→ Ax
berechnen kann.
Beobachtung: Die Zahl der notwendigen Iterationen vervierfacht sich bei Verdopplung der Gröÿe des Systems. Der Grund liegt in der Vervierfachung der Kondition der Matrix
An bei Verdopplung von n, was zu einer entsprechenden Verschlechterung der Konvergenzrate führt.
5.5.5 Optimale Dämpfung Bemerkungen:
•
Auch für die iterative Lösung linearer Gleichungssysteme sind Dämpfungsstrategien sinnvoll und wichtig, insbesondere weil sie gröÿere Freiheit in der Wahl des Vorkonditionierers
B •
erlauben.
Im Gegensatz zur Dämpfung bei nichtlinearen Problemen ist es hier sogar möglich: 1. einen optimalen Dämpfungsfaktor zu berechnen, 2. die Suchrichtung
v = B(y − Ax)
so zu verändern, dass der Schritt in einem
Unterraum, der auch alle alten Suchrichtungen beinhaltet, optimal ist.
•
Beispiele solcher Verfahren sind das CG-Verfahren (conjugate gradient) für symmetrisch positiv de nite Systeme oder das GMRES-Verfahren (generalized minimum residual) für eine beliebige Systemmatrix.
•
In Scilab/Matlab/Octave gibt es dazu die Befehle pcg (preconditioned conjugate gradient) und gmres.
Bemerkung: Ähnliche Matrizen wie die Matrix
An
des vorigen Abschnitts (welche man
aber nicht mehr so einfach direkt lösen kann) entstehen etwa bei der Diskretisierung sogenannter partieller Dierentialgleichungen. Mit diesen kann man etwa Strömungen oder die Deformation von Objekten unter Krafteinwirkung beschreiben. Die hohe Zahl der Unbekannten rührt von der Diskretisierung eines mehrdimensionalen Kontinuums mit Hilfe eines feinen Gitters her. Die schnelle iterative Lösung linearer Systeme mit Hilfe von CG- bzw. GMRES-Verfahren in Kombination mit einem guten Vorkonditionierer ist für eine erfolgreiche numerische Simulation solcher Probleme essentiell.
52
Übungen Aufgabe: Bei der Diskretisierung des Randwertproblems
−u00 (x) = f (x) ,
x ∈ (0, 1) ,
u(0) = u(1) = 0
mit sogenannten Finiten Dierenzen auf einem Gitter der Gitterweite die
h=
n × n-Tridiagonalmatrix 1 A= 2 h
2
−1
−1 ..
.
..
.
1 entsteht n+1
..
.
. .. . −1 −1 2 (k)
k = 1, . . . , n die Eigenvektoren v (k) durch vi = sin(ikπh) und die −2 Eigenwerte als λk = 2h (1−cos(kπh)) gegeben sind. Wie groÿ ist also die Spektralkondition Zeigen Sie, dass für dieser Matrix? Aufgabe: Das einfachste iterative Verfahren zur Lösung eines linearen Gleichungssystems
Ax = y
mit einer symmetrisch positiv de niten Matrix
x(k+1) = x(k) + B(y − Ax(k) ) ,
A
ist die Iteration
B=
ω kAk2
mit einem Dämpfungsfaktor 0 < ω < 2. Zur Konvergenzanalyse betrachten wir den (k) Fehler e = x(k) − x, wobei x die exakte Lösung ist. Zeigen Sie, dass diese Iteration in der Euklidischen Norm
k·k2
konvergiert.
53
6 Interpolation 6.1 Einführung Problem: Eine Funktion
f : [a, b] → R
ist nur an
n
Punkten
t1 , . . . , tn ∈ [a, b]
bekannt
(z.B. durch Messungen), d.h.
f (ti ) = fi , Nichtsdestoweniger wissen wir, dass
f
i = 1, . . . , n .
eigentlich eine stetige (glatte) Gröÿe beschreibt
und wollen daher
•
die Funktion an einem Wert
t ∈ [a, b] \ {t1 , . . . , tn }
•
die Funktion graphisch darstellen oder aber
•
Ableitungen der Funktion berechnen.
approximieren,
n-dimensionalen Vektorraum V bestehend aus Funktionen die n Bedingungen ϕ(ti ) = fi genau ein Element aus V
Idee: Wir wählen einen geeigneten
ϕ,
die so beschaen sind, dass
auswählen. Beispiele: In dieser Vorlesung betrachten wir nur folgende Situationen:
• f
ist glatt und der Raum
V
besteht aus Polynomen
• f
ist glatt und der Raum
V
besteht aus stückweisen Polynomen (Splines)
Bemerkung: Man beachte, dass die Werte von
f
an der Stelle
ti
hier als genau bekannt
angenommen und entsprechend exakt interpoliert werden. Ebenfalls häu g in der Praxis anzutreen ist das Problem, dass die Messungen der
fi = f (ti )
mit Fehlern behaftet
sind. In diesem Fall ist es meist nicht angebracht, die Gleichungen
ϕ(ti ) = fi
exakt zu
n Messungen ein k -dimensionaler Funktionenraum V mit dim(V ) = k n gewählt, und dann durch Lösung eines linearen Ausgleichsproblems ein Element aus V ausgewählt (siehe etwa das Automotoren-Beispiel aus dem Kapitel erfüllen. Stattdessen wird zu
über Ausgleichsprobleme).
54
6.2 Lagrange-Interpolation 6.2.1 Konstruktion Voraussetzung: Gegeben sei ein Intervall mit
ti 6= tj
für
i 6= j
(d.h. die
ti
I = [a, b]
und Stützstellen
Aufgabe: (Lagrangesche Interpolationsaufgabe) Zu einer Funktion n wir nun ein Polynom p ∈ P mit
p(ti ) = f (ti ) , Pn
{t0 , . . . , tn } ∈ [a, b]
sind paarweise verschieden).
f ∈ C 0 ([a, b])
suchen
i = 0, . . . , n .
bezeichnet hier den Vektorraum aller Polynome mit Grad kleiner oder gleich
Beachte: Ab jetzt ist
n
der Polynomgrad und
De nition: Zu den Stützstellen
Satz: Es gilt
Li,n (tj ) = δij ,
ti
n+1
n.
die Zahl der Stützstellen.
i = 0, . . . , n de nieren Q j6=i (t − tj ) Li,n (t) = Q j6=i (ti − tj ) für
i, j = 0, . . . , n. Hieraus P n bilden.
wir die Lagrange-Polynome
folgt insbesondere, dass die
n+1
Lagrange-Polynome eine Basis von Beweis: Wegen
Li,n (tj ) = δij
sieht man sofort, dass man keines der Lagrange-Polynome n aus den anderen durch Linearkombination erhalten kann. Da dim(P ) = n + 1 folgt die
Behauptung. Satz: Die Lagrangesche Interpolationsaufgabe besitzt eine eindeutige Lösung, d.h. zu 0 n jedem f ∈ C ([a, b]) gibt es genau ein Polynom p ∈ P , welches die Interpolationsbedingung erfüllt. Beweis: Existenz: Wir setzen einfach
p(t) =
n X
f (ti )Li,n (t) .
i=0
p˜, welches p˜(ti ) = f (ti ) erfüllt, ebenfalls diese Darstellung p˜ = p sein.
Eindeutigkeit: Weil jedes Polynom in der Lagrange-Basis hat, muss
Beobachtung: Die Interpolationsabbildung
I : C 0 ([a, b]) → P n , die jeder stetigen Funktion
f
f 7→ p
das interpolierende Polynom zu den paarweise verschiedenen
Stützstellen t0 , . . . , tn zuordnet, ist eine lineare Abbildung. (Es ist eine Projektion von C 0 ([a, b]) auf P n ⊂ C 0 ([a, b]).)
55
6.2.2 Gute Basiswahl Frage: Bezüglich welcher Basis von
Pn
sollen wir das Interpolationspolynom darstellen? n (Diese Frage ist unabhängig davon, ob P an sich geeignet ist.) Beobachtung: Die Koe zienten des Interpolationspolynoms in der Lagrange-Basis sind trivialerweise gut konditioniert, weil die lineare Abbildung
I˜ : C 0 ([a, b]) → Rn+1 ,
f 7→ ~f := (f (ti ))i=0,...,n ∈ Rn+1
oenbar der Abschätzung
k~f k∞ ≤ kf kC 0 ([a,b]) genügt, wobei
kf kC 0 ([a,b]) := sup |f (x)| . x∈[a,b]
Aber: Die Kondition der Koe zienten eines Interpolationspolynoms in der üblichen 2 Monombasis {1, x, x , . . .} ist oft schlecht konditioniert. Diese Koe zienten a0 , . . . , an ergeben sich nämlich als Lösung des Vandermonde-Systems
a0 + a1 t0 + . . . + an tn0 = f0 a0 + a1 t1 + . . . + an tn1 = f1 . . .
. . .
. . .
a0 + a1 tn + . . . + an tnn = fn oder in Matrixform
1 t0 · · ·
.. .. . . 1 tn · · ·
a0 f0 . . . . . = . . . . . tnn an fn tn0
Hier ist die Systemmatrix meist sehr schlecht konditioniert. Beispiel: Wir berechnen die Kondition des Vandermonde-Systems für äquidistante Stützstellen mit folgendem Scilab-Code: f u n c t i o n A = vand ( n ) x =[0:1/n : 1 ] ' ; A = z e r o s ( n +1 , n + 1 ) ; for
i =0: n
A ( 1 : n +1 , i +1) = x
.^
i ;
end endfunction for
i =1:10
d i s p ( i ) , d i s p ( c o n d ( vand ( i ) ) ) end
56
Beobachtung: Wir erhalten folgende Werte:
n cond2 (Vn )
1
2
3
4
5
10
2.6180
15.100
98.868
686.43
4924.4
1.1558 · 108
Folgerung: Man stellt Interpolationspolynome daher meist nicht bezüglich der Monombasis dar, sondern arbeitet mit der Lagrange-Basis oder der sogenannten Newton-Basis. De nition: Die Newton-Basis zu den Stützstellen
ωi (t) :=
i−1 Y (t − tj ) ∈ P i ,
t0 , . . . , tn
besteht aus den Polynomen
i = 0, . . . , n .
j=0
Bemerkungen:
•
Die Newton-Basis wird bei Hinzufügen neuer Interpolationspunkte erweitert, ohne dass sich die bisherigen Basiselemente verändern.
•
Die Bestimmung der Koe zienten in der Newton-Basis verwendet die sogenannten dividierten Dierenzen (siehe [Rannacher, Satz 3.1.2]).
•
Die dividierten Dierenzen können auch zur Auswertung des Interpolationspolynoms an einer bestimmten Stelle
x
verwendet werden (sog. Schema von Neville).
6.2.3 Interpolationsfehler
f ∈ C n+1 ([a, b]) und sei p ∈ P n das Interpolationspolynom zu f bezüglich t0 , . . . , tn ∈ [a, b]. Dann gibt es zu jedem x ∈ [a, b] ein ξx ∈ [a, b] mit
Satz: Sei Stellen
f (x) − p(x) =
Hierbei ist
der
f (n+1) (ξx ) ω(x) . (n + 1)!
n Y ω(x) := (x − ti ). i=0
Beweis: Für x ∈ {t0 , . . . , tn } ist die Aussage oenbar richtig. Für f (x)−p(x) wir Kx := und betrachten die Funktion ω(x)
x 6∈ {t0 , . . . , tn }
setzen
ψ(t) = f (t) − p(t) − Kx ω(t) . ψ besitzt mindestens die n+2 Nullstellen {x, t0 , . . . , tn }. Nach dem Satz von Rolle hat ψ 0 00 dann mindestens n + 1 (dazwischen liegende) Nullstellen, ψ mindestens n Nullstellen, (n+1) usw. bis ψ hat mindestens eine Nullstelle ξ = ξx ∈ [a, b]. Für diese gilt 0 = ψ (n+1) (ξx ) = f (n+1) (ξx ) − Kx (n + 1)!
57
woraus die Behauptung folgt. Beispiel: Wir betrachten die Funktion
f (x) =
1 1+x
[0, 1] und interpolieren sie an äquidistanten Stützstellen t0 = 0, t1 = h, . . . , k k! tn = nh mit h = n1 . Dann gilt f (k) (x) = (−1) . Wie man sich leicht überlegt, kann man (1+x)k |ω(x)| für x ∈ [0, 1] (sehr grob!) abschätzen als im Intervall
1 1 2 n h · · ··· ≤ |ω(x)| ≤ h(2h) · · · (nh) = 2 2n n n n
1 n/2 2 n
n
2− 2 = n
Somit gilt für den Interpolationsfehler über das ganze Intervall n
sup |f (x) − p(x)| ≤ x∈[0,1]
2− 2 . n
n und x fest vorgegeben hat, ist die Abschätzung noch einfacher. Z.B. ist der |f (k) (ξ)| 1 wegen supξ∈[0,1] ≤ 1 abschätzbar als Interpolationsfehler für n = 4 und x = 3 k! Wenn man
1 1 1 1 1 1 3 1 · | − | · | − | · | − | · | − 1| 3 3 4 3 2 3 4 3 1·1·1·5·2 ≈ 1.3 · 10−3 = 3 · 12 · 6 · 12 · 3
|f (x) − p(x)| ≤
Aber: Normalerweise fallen die Ableitungen nicht schnell genug ab, um eine Approximation
n → ∞ durch Polynominterpolation mit äquidistanten Stützstellen zu gewährleisten! 1 Ein bekanntes Gegenbeispiel ist auf dem gröÿeren Intervall [−5, 5], was wir im 1+x2 für
nächsten Abschnitt sehen werden.
6.2.4 Kondition der Polynominterpolation Frage: Ist die Interpolationsabbildung
I : C 0 ([a, b]) → P n
gut konditioniert, d.h. ist die
Interpolationsaufgabe gut gestellt? Antwort: Leider nein für groÿe
n, wenn man an einer Approximation in der Maximumsnorm
interessiert ist. Beobachtung: Für die Polynominterpolation
Cn :=
sup 06=f ∈C 0 ([a,b])
:=
I : C 0 ([a, b]) → P n
gilt
kpkC 0 ([a,b]) → ∞ (n → ∞) . kf kC 0 ([a,b]) p = I(f ) für groÿes n erhebliche I linear ist, bedeutet es gleichzeitig,
Dies bedeutet, dass die Lösung der Interpolationsaufgabe Überschwinger aufweisen kann. Weil die Interpolation
58
dass eine kleine Variation von
f
zu einer groÿen Änderung des interpolierenden Polynoms
führen kann. Beispiel: Wir interpolieren die Funktion
f (x) = auf dem Intervall
[a, b] = [−5, 5]
zu den äquidistanten Stützstellen
tk = a + kh , mit Polynomen
n-ten
1 1 + x2
k = 0, . . . , n ,
h=
b−a n
Grades.
Programm: (Polynominterpolation) Beobachtung:
2.0
1.5
1.0
0.5
0.0
−0.5 −6
Für groÿes
n
−4
−2
0
2
4
6
hat das interpolierende Polynom immer gröÿere Oszillationen an den
Intervallrändern. Die Stabilität ist oenbar nicht mehr gewährleistet. Abhilfe:
59
•
Wenn man auf die Äquidistanz der Stützstellen verzichtet, kann man ein wesentlich angenehmeres
•
log(n)-Wachstum
der Stabilitätskonstanten erreichen.
Verwende andere Funktionenräume zur Interpolation! Diesen Weg werden wir in der Folge gehen, wenn wir Räume betrachten, die nur stückweise polynomial (mit niedrigerem Grad) sind.
6.3 Stückweise polynomiale Approximationen (Splines)
∆ = {t0 , . . . , tn } des Intervalls I = ti ∈ [a, b] mit a = t0 < t1 < . . . < tn = b. Wir nennen die h := maxi=1,...,n |ti − ti−1 | die Feinheit der Unterteilung (oder auch Gitterweite).
De nition: Eine Unterteilung (auch Gitter genannt)
[a, b]
besteht aus Punkten
Gröÿe
Ansatz: Sei
∆ = {t0 , . . . , tn }
eine Unterteilung von
I = [a, b].
Wir suchen zu
f ∈ C 0 (I)
interpolierende Funktionen in folgendem Raum
k,r S∆ := {s ∈ C r (I) : s|[ti−1 ,ti ] ∈ P k } . (Also:
r-mal stetig dierenzierbar und auf jedem Teilintervall ein Polynom k -ten Grades.) k,r s ∈ S∆ soll die Bedingungen
Die Interpolationsbedingung ist wie gehabt:
s(ti ) = f (ti ) ,
i = 0, . . . , n
erfüllen.
6.3.1 Lineare Splines Ansatz: Der einfachste Fall ist
k =1
und
r =0
(d.h. die Funktionen sind stückweise
linear und global stetig):
1,0 S∆ := {s ∈ C 0 (I) : s|[ti−1 ,ti ] ∈ P 1 , i = 1, . . . , n} . Diesen Raum nennt man stückweise lineare Funktionen oder auch den Raum linearer Splines (Spline bedeutet etwa Stab auf Englisch). Beobachtung: Der Graph eines linearen Spline ist ein Polygonzug, der an den Stellen Knicke haben darf. Satz:
1,0 S∆
hat eine Basis bestehend aus den Hutfunktionen
ϕi (x) =
x−ti−1 ti −ti−1 ti+1 −x
t −t i+1 i 0
x ∈ [ti−1 , ti ] x ∈ [ti , ti+1 ] , sonst
60
i = 1, . . . , n − 1
ti
zusammen mit
( ϕ0 (x) = Diese Basis erfüllt
t1 −x t1 −t0
(
x ∈ [t0 , t1 ] , sonst
0
ϕi (tj ) = δij
ϕn (x) =
x−tn−1 tn −tn−1
x ∈ [tn−1 , tn ]
0
sonst
.
(sieht man einfach durch Einsetzen).
Beweis: Die Dimension des Raums aller (auch der unstetigen) stückweise linearen Funktionen 1,0 auf n Intervallen ist 2n. S∆ ist nun ein Teilraum davon, der durch n−1 Stetigkeitsbedingungen
i = 1, . . . , n − 1
lim s(t) = lim s(t) , t↑ti
t↓ti
eingeschränkt wird, was zu einer Dimension Zahl der oben de nierten Funktionen
n+1
führt. Dies entspricht aber genau der
ϕi . Diese sind wegen ϕi (tj ) = δij
linear unabhängig
und bilden daher eine Basis. Anwendung: Wegen ϕi (tj ) = 1,0 mit (I∆ f )(ti ) = f (ti ) durch
δij
ist zu
1,0 1,0 f ∈ C 0 ([a, b]) eine Interpolierende s = I∆ f ∈ S∆
1,0 I∆ f (x)
=
n X
f (ti )ϕi (x)
i=0 de niert.
f ∈ C 2 (I) und die Unterteilung ∆ von I habe die maximale h := maxi=1,...,n |ti − ti−1 |. Dann gilt für alle x ∈ I die Abschätzung
Satz: Es sei
1,0 |f (x) − I∆ f (x)| ≤
Gitterweite
h2 sup |f 00 (ξ)| . 8 ξ∈[a,b]
Beweis: Dies folgt sofort aus der Fehlerabschätzung, die wir für die Lagrange-Interpolation kennen. Auf jedem Teilintervall
|f (x) −
[ti−1 , ti ]
1,0 I∆ f (x)|
gilt für ein
ξ ∈ [ti−1 , ti ]
|f 00 (ξ)| = |(x − ti−1 )(x − ti )| . 2
|ω(x)| = |(x − ti−1 )(x − ti )| berechnen. Das ti−1 +ti wird aber oenbar beim Scheitelpunkt der Parabel x = im Intervallmittelpunkt 2 1 h2 2 angenommen und hat den Wert (ti − ti−1 ) ≤ . 4 4 Wir müssen also nur noch das Maximum von
6.3.2 Kubische Splines Beobachtung: Die linearen Splines waren nur stetig, an den Stützstellen können Knicke auftreten. Abgesehen davon, dass die Approximationsordnung relativ niedrig ist, ist das Fehlen an Dierenzierbarkeit für verschiedene Anwendungen problematisch. Insbesondere
61
für Graphikanwendungen ist nicht einmal die stetige Dierenzierbarkeit völlig zufriedenstellend, weil dem Auge auch ein Sprung in der zweiten Ableitung noch auallen kann.
k=3
Abhilfe: Wir betrachten den Ansatzraum für
und
r = 2,
d.h.
3,2 S∆ := {s ∈ C 2 (I) : s|[ti−1 ,ti ] ∈ P 3 , i = 1, . . . , n} . Diese Funktionen nennt man kubische Splines. Frage: Ist der Funktionenraum
•
3,2 S∆
zur Interpolation brauchbar, d.h. insbesondere
Ist das Problem
s(ti ) = f (ti ) , für beliebiges
•
0
f ∈ C (I)
Ist eine Lösung
Lemma: Es sei
lösbar?
3,2 s ∈ S∆
f ∈ C 2 ([a, b])
i = 0, . . . , n
eindeutig?
und
3,2 s ∈ S∆
mit
s(ti ) = f (ti ) ,
i = 0, . . . , n
Dann gelten die Gleichungen
Z
b
b
(f 00 (x) − s00 (x))s00 (x) dx = [(f 0 (x) − s0 (x))s00 (x)]a ,
a woraus folgt
b
Z
00
00
b
Z
2
00
|f (x) − s (x)| dx =
b
Z
2
a
a
b
|s00 (x)|2 dx − 2 [(f 0 (x) − s0 (x))s00 (x)]a .
|f (x)| dx − a
Beweis: Zweimalige partielle Integration auf den Teilintervallen liefert
b
Z
(f 00 (x) − s00 (x))s00 (x) dx
a
=
n X
0
0
[(f (x) − s (x))s
00
t (x)]tii−1
Z −
0
0
= [(f (x) − s (x))s
00
b (x)]a
−
(f 0 (x) − s0 (x))s000 (x) dx
a
i=1 n X
b
[(f (x) − s(x))s
000
t (x)]tii−1
b
Z
(f (x) − s(x))s(4) (x) dx
+ a
i=1 b
= [(f 0 (x) − s0 (x))s00 (x)]a wegen der Interpolationsbedingungen, und weil
s(4) ≡ 0
auf jedem Teilintervall ist. Nun
folgt die zweite Aussage aus der ersten wegen
Z
b 00
00
2
Z
b
|f (x) − s (x)| dx = a
00
2
Z
b
|f (x)| dx −
2
Z
|s (x)| dx − 2
a
Folgerung: Der die Funktion
00
a
f ∈ C 2 ([a, b])
b
(f 00 (x) − s00 (x))s00 (x) dx
a
interpolierende Spline
Randbedingungen
62
3,2 , s ∈ S∆
der eine der
1.
s0 (a) = f 0 (a)
2.
s00 (a) = s00 (b) = 0
und
s0 (b) = f 0 (b)
(vollständige Splineinterpolation)
(natürliche Splineinterpolation)
erfüllt, minimiert das Funktional (Funktional: Abbildung von Funktionen in die reellen Zahlen)
b
Z
|u00 (x)|2 dx
Q(u) := a 2
u ∈ C (I), die 1. für alle ti ∈ ∆ mit f übereinstimmen und 2. im 0 0 0 0 Fall der vollständigen Splineinterpolation auch u (a) = f (a) und u (b) = f (b) erfüllen. (Beweis: s interpoliert auch u.) unter allen Funktionen
Interpretation: Die Gröÿe des Graphen von
|u00 (x)| kann man für kleine Ableitung |u0 (x)| als die Krümmung
u interpretieren. Ein kubischer Spline nimmt also (unter Einhalten der
Interpolationsbedingungen) eine Art minimaler Krümmung an. Folgerung: Der eine Funktion
f ∈ C 2 ([a, b])
interpolierende Spline
3,2 s ∈ S∆ ,
welcher
zusätzlich eine der Randbedingungen 1.
s0 (a) = f 0 (a)
2.
s00 (a) = s00 (b) = 0
und
s0 (b) = f 0 (b)
(vollständige Splineinterpolation)
(natürliche Splineinterpolation)
erfüllt, ist eindeutig bestimmt.
Rb s1 und s2 gäbe, so müssten beide Q(u) = a |u00 (x)|2 dx minimieren, woraus insbesondere Q(s1 ) = Q(s2 ) folgt. Nach dem Lemma ergibt sich dann auch Z b Q(s1 − s2 ) = |s001 (x) − s002 (x)|2 dx = 0
Beweis: Wenn es zwei solche Splines
a
w = s1 − s2 ist damit linear auf jedem Teilintervall [ti−1 , ti ] w(ti ) = 0 für i = 0, . . . , n. Hieraus folgt aber w ≡ 0.
folgt. Die Funktion erfüllt obendrein
Folgerung: Für jedes
und
3,2 s ∈ S∆
zu
ist ein Teilraum der stückweise kubischen Polynome, welcher
4n
f ∈ C 2 ([a, b])
gibt es genau eine Splineinterpolierende
vollständigen (alternativ: natürlichen) Randbedingungen. Beweis: Der Raum
3,2 S∆
3n−3 lineare Übergangsbedingungen 3,2 eine Dimension von S∆
Freiheitsgerade (Dimensionen) hat. Eingeschränkt wird er durch (Stetigkeit von Funktion, erster und zweiter Ableitung), was von
n+3
liefert. Diese werden dann durch
n+1
Interpolationsbedingungen und zwei
Randbedingungen spezi ziert. Dieser letzte Prozess entspricht dem Lösen eines linearen Systems der Dimension
n+3, für welches die Lösungen (wie wir gezeigt haben) eindeutig
sind. Ein Standardergebnis der linearen Algebra sagt, dass dies äquivalent zur Existenz einer Lösung für alle Vorgaben ist.
63
6.3.3 Berechnung eines interpolierenden kubischen Splines 3,2 Beobachtung: Die zweite Ableitung eines kubischen Splines s ∈ S∆ ist ein stückweise 1,0 00 00 00 linearer Spline s ∈ S∆ . s ist daher durch die Werte Mi = s (ti ) eindeutig bestimmt. Man nennt diese Gröÿen die Momente.
3,2 1,0 00 s ∈ S∆ seien die Momente m = s ∈ S∆ und die Werte fi = s(ti ) i−1 = t−thi−1 Teilintervall [ti−1 , ti ] gilt mit ˆ t = tt−t i −ti−1 i
Satz: Für einen Spline gegeben. Auf dem
mi (t) = m ˆ i (ˆt) = (1 − ˆt)Mi−1 + ˆtMi und
h2ˆt(1 − ˆt) si (t) = ˆsi (ˆt) = (1 − ˆt)fi−1 + ˆtfi − i (1 + ˆt)Mi + (1 + (1 − ˆt))Mi−1 6 Beweisskizze: Bei gegebenen Momenten
m = s00
kann der Spline
s auf jedem Teilintervall
durch zweifache Integration erhalten werden: es gilt nämlich
Z
t
Z
τ
s(t) = ti−1 und die Werte von
s00 (r)dr dτ + At + B
ti−1
A und B können so gewählt werden, dass s die Interpolationsbedingungen
an den Intervallenden erfüllt. Für die genaue Rechnung wird auf die Übung verwiesen. Satz: Es sei
I = [a, b] und ∆ = {t0 , . . . , tn } wie oben. Dann sind die Momente (Mi )i=0,...,n 3,2 f ∈ C 1 ([a, b]) vollständig interpolierenden kubischen Splines s ∈ S∆
des eine Funktion
die Lösung des Gleichungssystems
hi hi + hi+1 hi+1 fi+1 − fi fi − fi−1 − , Mi−1 + Mi + Mi+1 = 6 3 6 hi+1 hi h1 h1 f1 − f0 M0 + M1 = − f 0 (a) , 3 6 h1 hn hn fn − fn−1 Mn + Mn−1 = f 0 (b) − , 3 6 hn
i = 1, . . . , n − 1 i=0 i = n.
Beweis: Die Gleichungen ergeben sich durch die Forderung nach stetiger Dierenzierbarkeit und den vollständigen Randbedingungen. Auf
s0i (t) =
[ti−1 , ti ]
gilt
1 0ˆ ˆs (t) = hi fi − fi−1 hi d ˆ d − (t(1 − ˆt)(1 + ˆt))Mi + (ˆt(1 − ˆt)(1 + (1 − ˆt)))Mi−1 hi 6 dˆt dˆt
64
Insbesondere haben wir
1 0 fi − fi−1 hi ˆs (0) = − (Mi + 2Mi−1 ) hi hi 6 1 fi − fi−1 hi s0i (ti ) = ˆs0 (1) = + (2Mi + Mi−1 ) hi hi 6
s0i (ti−1 ) =
s0i (ti ) = s0i+1 (ti ) übersetzt sich dann in die Gleichungen für i = 1, . . . , n−1, 0 0 0 0 die Randbedingung s1 (a) = f (a) liefert die Gleichung für i = 0 und sn (b) = f (b) die Gleichung für i = n. Die Forderung
Bemerkung: Das entstehende Gleichungssystem hat eine stark diagonaldominante Systemmatrix
A,
d.h.
P
j6=i
|Aij | < Aii .
Hieraus folgt insbesondere die Regularität (siehe Übung).
6.3.4 Fehlerabschätzung der kubischen Splines
I , ∆ wie bisher, und k, l ∈ N mit 0 ≤ l ≤ 2 ≤ k ≤ 4 und f ∈ C k (I). Dann 3,2 f vollständig (!) interpolierende kubische Spline s ∈ S∆ die Abschätzung
Satz: Seien erfüllt der
max | x∈I
Spezialfall: Für
k=4
und
dk f dl s dl f k−l (x) − (x)| ≤ Ch max | (x)| . x∈I dxk dxl dxl l=0
ergibt sich, dass für
f ∈ C 4 (I)
max |f (x) − s(x)| ≤ Ch4 max | x∈I
x∈I
gilt
d4 f (x)| dx4
Bemerkung: Für den interpolierenden Spline mit natürlichen Randbedingungen gelten die obigen optimalen Abschätzungen nur, wenn auch f 00 (a) = f 00 (b) = 0 erfüllt.
f
die natürlichen Randbedingungen
6.3.5 Anwendung 1 auf dem Intervall 1+x2 diesmal mit Hilfe von vollständig interpolierenden kubischen Splines.
Beispiel: Wir interpolieren wieder die Funktion
f (x) =
Programm: function y = 1
y = f ./
(x)
(1 + x
.∗
x) ;
endfunction f u n c t i o n M = compute_moments
(t ,
f ,
n = length ( t ) ; A = zeros (n , n) ;
65
Dfa ,
Dfb )
[−5, 5],
b = zeros (n , 1 ) ; //
setup
of
the
interior
rows
i =2: n−1
for
A( i , i ) = ( t ( i +1)− t ( i
− 1) ) / 3 ; A( i , i − 1) = ( t ( i )− t ( i − 1) ) / 6 ; A( i , i +1) = ( t ( i +1)− t ( i ) ) / 6 ; b ( i ) = ( f ( i +1)− f ( i ) ) / ( t ( i +1)− t ( i ) ) −
( f ( i )− f ( i
− 1) ) / ( t ( i )− t ( i − 1) )
; end //
setup
of
the
boundary
rows
h = ( t ( 2 )− t ( 1 ) ) ; A( 1 , 1 ) = h / 3 ; A( 1 , 2 ) = h / 6 ; b ( 1 ) = ( f ( 2 )− f ( 1 ) ) /h−Dfa ; h = ( t ( n )− t ( n − 1) ) ; A( n , n ) = h / 3 ; A( n , n − 1) = h / 6 ; b ( n ) = Dfb −( f ( n )− f ( n − 1) ) / h ; //
finally
we
solve
the
system
M = A\ b ; endfunction function //
y = spline_eval
evaluates
the
(t ,
spline
f , M,
x)
specified
by
t , f ,M a t
x
n = length ( t ) ; //
find
for
the
subinterval
in
which
local
cubic
x
lies
i =2: n
if
( t ( i )>=x ) break ;
end end //
evaluation
of
the
h = t ( i )− t ( i
− 1) ; t a u = ( x−t ( i − 1) ) / h ; y = (1 − t a u ) ∗ f ( i − 1) + t a u ∗ f ( i ) . . . − h ∗ h ∗ t a u ∗ (1 − t a u ) / 6 ∗ (M( i ) ∗ (1+ t a u ) +(2− t a u ) ∗M( i − 1) ) ;
endfunction function //
y = spline_multiple_eval
evaluates
a
spline
for
each
( t , f ,M, x ) vector
n = length (x) ; y = zeros (n , 1 ) ; for
i =1: n
y( i ) = spline_eval
( t , f ,M, x ( i ) ) ;
end endfunction
66
component
of
x
function //
we
//
the
testit use
(t , f )
numerical
boundary
differentiation
values
of
the
for
complete
x= −5;
Dfa =( f ( x+1e − 8)− f ( x ) ) /1 e − 8;
x =5;
Dfb=( f ( x+1e − 8)− f ( x ) ) /1 e − 8;
obtaining spline
m = compute_moments ( t , f ( t ) , Dfa , Dfb ) ;
−5:0.1:5] ';
x=[
y = s p l i n e _ m u l t i p l e _ e v a l ( t , f ( t ) ,m, x ) ; plot (x , [ y , f (x) ] ) ; endfunction c l f () ; testit ([
−5:1:5] , f )
Beobachtung: Schon bei 11 Stützstellen (h
= 1) sieht man kaum einen Unterschied mehr
zwischen der Funktion und dem Spline.
Übungen
I = [a, b], h := b − a und s ∈ P 3 sei ein Polynom dritten Grades, von dem 00 00 Werte s(a) = sa , s(b) = sb , s (a) = Ma , s (b) = Mb kennen. Rechnen Sie nach,
Aufgabe: Sei wir die dass
s
sich darstellen lässt als
s(t) = (1 − ˆt)sa + ˆtsb − wobei
ˆt =
h2 ˆ t(1 − ˆt) (1 + ˆt)Mb + (1 + (1 − ˆt))Ma , 6
t−a . h
Aufgabe: Eine Matrix
A ∈ Rn×n heiÿt stark diagonaldominant, X |Aij | < |Aii | , i = 1, . . . , n .
wenn
j6=i Zeigen Sie, dass
A
dann regulär ist, und dass das sogenannte Jacobi-Verfahren
xk+1 = xk + D−1 (b − Axk )
=Diagonale
Ax = b in P j6=i |Aij | ρ ≤ max i=1,...,n |Aii |
zur iterativen Lösung des Gleichungssystems Konvergenzrate
(D
konvergiert.
67
von
A)
der Maximumsnorm mit einer
7 Numerische Quadratur 7.1 Einführung Das Ziel dieses Kapitels ist die Berechnung des Riemann-Integrals
b
Z
f (x) dx .
I(f ) = a Bezeichnung: Manchmal werden wir statt
I
auch
I[a,b]
schreiben, insbesondere wenn wir
unterschiedliche Integrationsbereiche betrachten. Beobachtung: Die Integration 1.
I
ist linear:
2.
I
ist positiv:
I
ist ein Funktional auf
C 0 ([a, b])
mit
I(αf + βg) = αI(f ) + βI(g) I(f ) ≥ 0
für
f ≥ 0.
3. Es sei
ˆt 7→ ˆt = a + ˆt(b − a) .
Φ : [0, 1] → [a, b] ,
Dann gilt nach dem Transformationssatz
I[a,b] (f ) = (b − a)I[0,1] (f ◦ Φ) . Problem: Normalerweise gibt es zu ausdrückbare Stammfunktion
F,
f ∈ C 0 ([a, b])
keine in elementaren Funktionen
mit Hilfe derer man
I(f ) = F (b) − F (a)
ausrechnen
könnte. Abhilfe: Man berechnet das Integral näherungsweise. Dies nennt man numerische Integration oder numerische Quadratur.
7.2 Quadraturformeln De nition: Eine Quadraturformel
Q
zur Approximation von
Q(f ) =
n X
Rb a
f (t) dt
hat die Form
λi f (ti )
i=0 mit den paarweise verschiedenen Knoten
{t0 , . . . , tn } ⊂ [a, b] und den Gewichten λ0 , . . . , λn ∈
R. Bemerkung: Eine solche Quadraturformel Kriterien: Eine gute Wahl der
λi , ti
Q
ist oenbar linear in
sollte folgendes erfüllen:
68
f.
1. 2.
Pn
i=0
λi = b − a. Das stellt sicher, dass wenigstens f (t) ≡ 1 korrekt integriert wird.
λi > 0
für
i = 0, . . . , n.
Das sichert
λi < 0 ist, Q(f ) < 0. Wie?)
(Umgekehrt: Wenn ein mit
f ≥0
und
f ≥ 0 ⇒ Q(f ) ≥ 0. kann man sehr leicht ein
f ∈ C 0 ([a, b])
nden
Bemerkung: Mit Hilfe des oben de nierten linearen Isomorphismus Φ : [0, 1] → [a, b] Pn ˆ ˆ kann man ausgehend von einer Quadraturformel Q[0,1] (f ) = i=0 λi f (ti ) für I[0,1] eine Quadraturformel
Q[a,b]
für
I[a,b]
Q[a,b] (f ) =
erzeugen als
n X
λi f (ti ) ,
ˆi . ti = Φ(ˆti ) , λi = (b − a)λ
i=0
7.3 Konstruktion durch Polynominterpolation Idee: Interpoliere den Integranden durch ein Polynom
n-ten Grades und integriere dieses
interpolierende Polynom. De nition: Zu den paarweise verschiedenen Stützstellen {t0 , . . . , tn } ⊂ [a, b] seien Li,n Pn n wieder die Lagrange-Polynome, und pf = i=0 f (ti )Li,n ∈ P sei das interpolierende 0 Polynom zu einer Funktion f ∈ C ([a, b]) mit pf (ti ) = f (ti ). Mit
Z Q(f ) =
b
pf (t) dt = a
n X
Z f (ti )
b
Li,n (t) dt a
i=0
haben wir oenbar eine Quadraturformel mit Knoten ti und Gewichten
λi =
Rb a
Li,n (t) dt
gefunden. Satz: Das so konstruierte gilt
Q ist exakt für alle Polynome n-ten Grades, d.h. für alle p ∈ P n
Q(p) = I(p).
P n . Nach De nition ist Q exakt für jedes der Li,n n und daher wegen der Linearität von I und Q auch für alle p ∈ P . Pn Bemerkung: Es gilt oenbar auch die Umkehrung: Wenn Q(f ) = i=0 λi f (ti ) für alle n f ∈ P exakt ist, so folgt Z b λi = Q(Li,n ) = I(Li,n ) = Li,n (t) dt .
Beweis:
{Li,n }i=0,...,n
ist eine Basis von
a
Φ : [0, 1] → [a, b] der oben de nierte lineare Isomorphismus. {ˆt0 , . . . , ˆtn } ⊂ [0, 1] und {t0 , . . . , tn } ⊂ [a, b] seien Knoten, welche Φ(ˆti ) = ti für i = 0, . . . , n erfüllen. ˆ i,n (ˆt) = Li,n (t) für t = Φ(ˆt). Für die Dann gilt für die zugehörigen Lagrange-Polynome L Satz: Es sei wieder
69
über Polynominterpolation zu den obigen Knoten erzeugten Quadraturformeln Pn und Q[a,b] = i=0 λi f (ti ) gilt daher
Q[0,1] =
Pn ˆ ˆ i=0 λi f (ti )
Z
b
Z
1
Li,n (t) dt = (b − a)
λi =
ˆi . ˆ i,n (ˆt) dt = (b − a)λ L
0
a
Φ linear ist, ist Li,n ◦ Φ wieder ein Polynom n-ten Grades, welches (Li,n ◦ Φ)(ˆtj ) = Li,n (tj ) = δij erfüllt. Wegen der eindeutigen Lösbarkeit der Interpolationsaufgabe ˆ i,n gelten. Dann folgt aber aus dem Transformationssatz muss Li,n ◦ Φ = L
Beweis: Da
b
Z
1
Z
ˆ i,n (ˆt) dˆt . L
Li,n (t) dt = (b − a)
λi = a
0
7.3.1 Newton-Cotes-Formeln De nition: Die Quadraturformeln, die durch Polynominterpolation mit Verwendung b−a äquidistanter Stützstellen ti = a+i , i = 0, . . . , n entstehen, heiÿen Newton-Cotes-Formeln. n
Bezeichnung: Für Beispiel: Für
u ∈ C 0 ([a, b])
n = 1, 2, 3, 4
3 4
6 1 8 7 90
x∈[a,b]
erhalten wir mit
ˆ i = λi /h n λ 1 1 1 2 2 1 4 1 2
kuk∞ := max |u(x)|.
setzen wir
h=b−a
Name Trapezregel
6 6 3 3 1 8 8 8 32 12 32 90 90 90
7 90
Simpson-Regel 3 Newtonsche -Regel 8 Milne-Regel
Fehler |I(f ) − h3 kf (2) k∞ 12 5 h kf (4) k∞ 2880 5 h kf (4) k∞ 6480 7 h kf (6) k∞ 1935360
Q(f )|
Bemerkungen:
•
Für die Trapezregel
T (f ) = h2 (f (a)+ f (b)) folgt aus der Interpolationsabschätzung kf − pf k∞ ≤
h2 00 kf k∞ 8
für das an den Endpunkten interpolierende lineare Polynom
pf ∈ P 1
sofort die
schon recht gute Fehlerabschätzung
Z |I(f ) − T (f )| = |
b
Z (pf (x) − f (x)) dx| ≤
a
|pf (x) − f (x)| dx ≤ a
70
b
h3 00 kf k∞ . 8
•
Analog erhält man für die Newton-Cotes-Formel der Ordnung
n die Fehlerabschätzung
|I(f ) − QN C,n (f )| ≤ C(n)hn+2 kf (n+1) k∞ . Man beachte, dass man durch die Integration einen Faktor
h
mehr erhält als bei
der Interpolationsfehlerabschätzung.
•
Durch eine genauere Analyse kann man noch die Fehlerkonstante verbessern. Unten zeigen wir diese Technik exemplarisch für die Trapezregel.
•
Interessant ist zudem die erhöhte Ordnung der Newton-Cotes-Regeln für gerades
n.
Unten beweisen wir diese Ordnungsverbesserung exemplarisch für die Simpson-Regel.
•
Ab
n = 7 treten negative Gewichte auf. Spätestens dann werden diese Newton-Cotes-Formeln
problematisch, weil Auslöschung schon bei der Integration einfacher konstanter Funktionen auftreten kann.
Satz: Für
f ∈ C 2 ([a, b]) genügt die Trapezregel T (f ) =
b−a (f (a)+f (b)) der Abschätzung 2
h3 00 |T (f ) − I(f )| ≤ kf k∞ . 12 Beweis: Es sei
pf ∈ P 1
die lineare Interpolierende von
Nach unserem Approximationssatz gilt für
|f (x) − pf (x)| ≤
f
zu den Stützstellen
a
und
b.
x ∈ [a, b]:
kf 00 k∞ (x − a)(b − x) . 2
Somit folgt aber auch
Z
b
Z
b
Z
b
|f (x) − pf (x)| dx pf (x) dx| ≤ f (x) dx − a a a Z Z 1 00 kf 00 k∞ b kf 00 k∞ 3 kf k∞ ≤ (x − a)(b − x) dx = h x ˆ (1 − x ˆ ) dx = h3 . 2 2 12 a 0
|T (f ) − I(f )| = |
Satz: Die Simpsonregel
S(f ) =
b−a b+a (f (a) + 4f ( ) + f (b)) 6 2
integriert sogar alle Polynome dritten (!) Grades exakt. Hieraus folgt die Fehlerabschätzung
|S(f ) − I(f )| ≤ Ch5 kf (4) k∞ .
71
Beweis: Ohne Einschränkung sei a = −b (eine einfache Verschiebung ändert nichts an I(f ), Q(f ) oder f (4) ). Wir wissen, dass die Simpsonregel S alle Polynome zweiten Grades 3 exakt integriert. Der Clou ist nun, dass sie auch das Polynom x wegen Q(f ) = I(f ) = 0 3 2 3 exakt integriert (Symmetrie!). Nun ist aber P = P ⊕ Span(x ), somit werden alle Polynome dritten Grades exakt integriert. Die Fehlerabschätzung ergibt sich dann wie folgt:
pf
sei die Interpolierende dritten Grades zu äquidistanten Stützstellen. Für diese
gilt
kf − pf k∞ ≤ Ch4 kf (4) k∞ , woraus wie gehabt folgt
b
Z
Z (f − pf ) dx| ≤
|I(f ) − S(f )| = |
b
|f (x) − pf (x)| dx ≤ Ch5 kf (4) k∞ .
a
a
7.4 Zusammengesetzte Formeln Problem: Die Polynominterpolation ist für hohen Polynomgrad (wenigstens für äquidistante Knoten) instabil. Dies überträgt sich auf die Quadraturformeln.
a = x0 < x1 < . . . < xn = b eine Zerlegung des Intervalls [a, b] in Teilintervalle [xi−1 , xi ]. Das Gesamtintegral I[a,b] ergibt sich dann als Summe der I[xi−1 ,xi ] , die lokal mit Idee: Sei
einer Newton-Cotes-Formel niedrigen Grades approximiert werden können. Bezeichnung: Eine solche Quadraturformel nennt man eine zusammengesetzte Quadraturformel.
Beispiele: Wir betrachten äquidistante Unterteilungen
•
mit
h=
b−a . n
Zusammengesetzte Trapezregel:
Th (f ) =
n X
Trapez
Q[xi−1 ,xi ] (f ) =
i=1
•
xi = a + ih
n X h i=1
n−1 f (a) + f (b) X (f (xi−1 ) + f (xi )) = h + f (xi ) 2 2 i=1
Zusammengesetzte Simpsonregel:
Sh (f ) =
n X
QSimpson [xi−1 ,xi ] (f ) =
i=1
n X h i=1
=
6
(f (xi−1 ) + 4f (
xi−1 + xi ) + f (xi )) 2
n−1 n X h f (a) + f (b) X xi−1 + xi + f (xi ) + 2 f( ) . 3 2 2 i=1 i=1
Satz: Die zusammengesetzte Trapezregel für eine äquidistante Unterteilung b−a mit h = genügt der Fehlerabschätzung n
|Th (f ) − I(f )| ≤ (b − a)
72
h2 00 kf k∞ . 12
xi = a + ih
Beweis: Es gibt
n=
h3 b−a Teilintervalle und auf jedem kann der Fehler durch kf 00 k∞ h 12
abgeschätzt werden. Bemerkung: Ebenso ergibt sich für die Simpson-Regel die Abschätzung
|Sh (f ) − I(f )| ≤ (b − a)
Bemerkung: Man beachte den Verlust einer
h4 kf (4) k∞ . 2880
h-Potenz gegenüber der lokalen Abschätzung.
7.5 Ausblicke Wir konnten einige interessante und wichtige Bereiche nicht behandeln. Dazu gehören unter anderem folgende Punkte:
•
Wenn man auf die Äquidistanz der Knoten verzichtet, kann man durch geschickte
s 2s − 1 exakt integrieren. Dies ist die sogenannte
Wahl derselben durch Polynominterpolation Quadraturformeln erzeugen, die bei Knoten noch Polynome vom Grad
Gauÿ-Quadratur. Die einfachste solche Formel ist die Mittelpunktsregel
QMP (f ) = f •
a + b (b − a) . 2
Wenn der Integrand glatt ist und die Unterteilungen regelmäÿig sind, kann man etwa für die zusammengesetzte Trapezregel ein Fehlerverhalten wie
Qh (f ) = I(f ) + C1 h2 + C2 h4 + . . . zeigen. Dies kann man ausnutzen, indem man sukzessive Approximationen Q(h), Q( h2 ), . . . durch ein Polynom in der Variablen x = h2 interpoliert, welches man dann an x = 0 auswertet, um eine Approximation höherer Ordnung zu I(f ) zu erhalten. Diese Technik heiÿt Extrapolation und kann in verschiedenen Situationen angewendet werden, z.B. auch bei der numerischen Dierentiation (siehe Übung 10/1). Die Anwendung auf die numerische Quadratur heiÿt Romberg-Quadratur. Auch Übung 11/1 kann man als eine Stufe einer solchen Extrapolation ansehen.
•
Für viele Integranden ist eine nicht äquidistante Wahl des Gitters von Vorteil. Die automatische Gitteranpassung an einen bestimmten Integranden führt dabei zu adaptiven Quadraturverfahren. Allerdings geht dadurch normalerweise die Möglichkeit zur Extrapolation verloren.
73
Übungen Aufgabe: Man kann die Ableitung einer Funktion
f ∈ C 2 (R)
an der Stelle
x
durch den
Dierenzenquotienten
Dh f (x) =
f (x + h) − f (x) h
approximieren. Wir untersuchen nun, wie gut die Approximation in Abhängigkeit von
h
ist. 1. Wie groÿ ist der Fehler dieser Approximation in führender Ordnung in
h bei exakter
Arithmetik? (Hinweis: Taylor-Entwicklung.) 2. Durch Auslöschung entsteht bei der Auswertung des Zählers in Gleitkommaarithmetik ein Fehler der Gröÿe |f (x)|ε. Wie sollte man demnach f 0 (x) möglichst gut ist?
h wählen, damit die Approximation
von
3. Testen Sie Ihre Vorhersage numerisch am Beispiel
f (x) = ex
und
x = 1.
Aufgabe: Eine genauere Untersuchung zeigt, dass für die Trapezsumme (zusammengesetzte 4 Trapezregel) angewandt auf f ∈ C ([a, b]) mit einer äquidistanten Unterteilung xi = a + ih für i = 0, . . . , n und h = n1 mit n ∈ N folgende Entwicklung gilt
T (h) = I(f ) + Ch2 + O(h4 ) , wobei
I(f ) =
Rb
a um Konstanten
f (x) dx und C ∈ R eine Konstante ist. Verwenden Sie diese Entwicklung, α, β ∈ R zu bestimmen, für welche gilt h αT (h) + βT ( ) = I(f ) + O(h4 ) . 2
Schreiben Sie die entstehende Quadraturformel in der Form
Q(h) = . . . f (x0 ) + . . . f (
x0 + x1 x1 + x2 ) + . . . f (x1 ) + . . . f ( ) + . . . + . . . f (xn ) . 2 2
Was ergibt sich?
74
8 Numerische Lösung gewöhnlicher Dierentialgleichungen 8.1 Grundlagen
f ∈ C 0 (R×Rn , Rn ). Zu einem gegebenen Startwert u0 ∈ Rn 1 n wir eine Funktion u ∈ C ([0, T ], R ) mit
De nition: Es sei
T >0
suchen
du (t) = f (t, u(t)) , dt u(0) = u0 .
und einem
t ∈]0, T [ ,
Dies nennen wir eine gewöhnliche Dierentialgleichung. Bemerkungen:
•
den Fall
•
n=1
Manchmal bezeichnet man nur den Fall
Wenn
f
n>1
als Dierentialgleichung und nennt
ein System gewöhnlicher Dierentialgleichungen.
nicht von
u
abhängt, so kann man die Lösung der Dierentialgleichung
durch eine einfache Integration erhalten:
Z
t
f (τ ) dτ .
u(t) = 0
•
Wenn
f
nicht von
t
abhängt, so nennt man die zugehörige Dierentialgleichung
autonom.
•
In vielen Anwendungen hat die Variable
t
die Bedeutung der Zeit. Die Ableitung nach der Zeit bezeichnen vor allem Physiker oft mit einem Punkt (u(t) ˙ := du (t)). dt
8.2 Beispiele 8.2.1 Pendelgleichung
75
Erinnerung: In der Einführung hatten wir bereits Dierentialgleichungen für das Stabpendel hergeleitet. Diese waren
˙ = ω(t) θ(t) g ω(t) ˙ = − sin θ(t) L θ den Auslenkungswinkel, ω die Winkeländerungsgeschwindigkeit, g die Gravitationskonstante,
wobei und
L
die Länge des Stabs bezeichnet.
Transformation: Wir setzen
u(t) = θ(t), ω(t)
t
und erhalten
u(t) ˙ = f (u(t))
mit
f (u) =
u2 g − L sin u1
.
8.2.2 Räuber-Beute-Modell Modell: Auf einer Insel gebe es Populationen von Raubtieren
B = B(t),
R = R(t)
und Beutetieren
die sich mit der Zeit ändern können. Ein einfaches Modell für die zeitliche
Veränderung der Populationen ist
˙ B(t) = µB B(t) − αB(t)R(t) ˙ R(t) = −µR R(t) + βB(t)R(t) Hierbei sind
µB , µR , α, β
positive reelle Konstanten.
Skalierung: Um die Eigenschaften des Modells zu studieren, wählt man oft die Einheiten so, dass die Gleichung möglichst einfach wird. Im obigen Fall lässt sich durch Wahl der Einheiten von
B, R
und
t
folgende Form erreichen:
˙ B(t) = B(t) − B(t)R(t) ˙ R(t) = −µR(t) + B(t)R(t)
Bemerkungen:
•
Man sieht, dass das Modell in seinem Verhalten im wesentlichen von einem einzigen Parameter
•
µ∈R
abhängt.
In Anwendungen ist es oft praktischer, die üblichen Einheiten zu verwenden. Hier verzichtet man daher meist auf solche Vereinfachungen.
Transformation: Wir setzen
t u(t) = B(t), R(t)
und erhalten
u(t) ˙ = f (u(t))
mit
f (u) =
76
u1 − u1 u2 −µu2 + u1 u2
.
8.3 Geometrische Interpretation Problem: Der Einfachheit halber betrachten wir die autonome Dierentialgleichung
u(t) ˙ = f (u(t)) , t ∈ [0, T ] , mit
u(0) = u0
f ∈ C 0 (Rn , Rn ).
Bemerkung: Sowohl die Pendelgleichung als auch das Räuber-Beute-Modell sind autonom. Nicht autonome Varianten würden sich ergeben, wenn etwa der Aufhängungspunkt des Stabpendels sich bewegen würde, oder wenn Geburts- oder Sterberaten im Räuber-Beute-Modell von der Zeit abhingen. Interpretation: Die Lösung 1. den Anfangspunkt 2. für alle
t
u : [0, T ] → Rn
u(0) = u0
die Ableitung
du (t) dt
beschreibt eine Kurve, die
hat, und die
= f (u(t))
hat.
n Bezeichnung: Die Abbildung f : R → Rn nennt man ein Vektorfeld. Kurven u : (t) = f (u(t)) erfüllen, nennt [a, b] → Rn , die in jedem ihrer Punkte die Gleichung du dt man Integralkurven des Vektorfelds. g = 1) sieht folgendermaÿen aus L (horizontale Achse: Winkel, vertikale Achse: Winkelgeschwindigkeit):
Beispiel: Das Vektorfeld für die Pendelgleichung (
4
2
0
-2
-4
-4 -2 0 2 vf
4
Im Falle des Räuber-Beute-Modells (µ
= 1)
vertikale Achse: Räuber):
77
ergibt sich (horizontale Achse: Beute,
3
2.5
2
1.5
1
0.5
0 0 0.5 1 1.5 2
vf
2.5 3
8.4 Numerische Lösung Problem: Sehr selten kann man für eine gewöhnliche Dierentialgleichung eine geschlossene Lösung der Form
u(t) = . . . Ausdruck
in
t
...
angeben. Genau wie bei der Integration ist man meist auf die numerische Lösung angewiesen.
8.4.1 Das explizite Euler-Verfahren Problem: Wir suchen eine Approximation der Lösung
u : [0, T ] → Rn
der gewöhnlichen
Dierentialgleichung
u(t) ˙ = f (t, u(t)) , t ∈ [0, T ] , mit
u(0) = u0
f ∈ C 0 (R × Rn , Rn ).
Algorithmus: (Explizites Euler-Verfahren) Unterteile [0, T ] in N Intervalle der Länge T n . Approximiere u : [0, T ] → R dann an den Stellen tk = kh durch folgende Werte h= N n yk ∈ R :
y0 = u0 ,
yk = yk−1 + hf (tk−1 , yk−1 ) ,
k = 1, . . . , N .
Bemerkungen:
•
Der Aufwand des Verfahrens ist im wesentlichen gegeben durch der Funktion
f.
78
N
Auswertungen
•
Wenn man ganz
[0, T ]
u : [0, T ] → Rn
nicht nur an den Stellen tk
de nierte Funktion approximieren will,
= kh sondern durch eine auf kann man die Punkte (tk , yk )
einfach linear interpolieren.
•
Angewandt auf ein Integrationsproblem führt das explizite Euler-Verfahren auf die Quadraturformel
Z
T
f (t) dt ∼ h(f (t0 ) + f (t1 ) + . . . + f (tN −1 )) =: QEE (h) . 0 Diese Formel hat im wesentlichen den gleichen Aufwand wie Trapez- oder Mittelpunktregel, besitzt aber im Gegensatz zu diesen Verfahren nur die Fehlerordnung O(h2 ). Es macht daher für die einfache Quadratur wenig Sinn.
O(h)
statt
8.4.2 Verfahren höherer Ordnung Problem: Die Fehlerordnung
O(h)
ist für viele Anwendungen nicht genau genug.
Abhilfe: Man konstruktiert Verfahren höherer Ordnung. Hierfür gibt es verschiedene Ansätze, von denen wir nur wenige Verfahren vorstellen. Algorithmus: Eine Erweiterung der Trapezregel ist das Verfahren von Heun:
y0 = u0 ,
h yk = yk−1 + (f (tk−1 , yk−1 ) +f (tk , yk−1 + hfk−1 )) . {z } 2 | fk−1
Dies hat pro Schritt den doppelten Aufwand wie das explizite Euler-Verfahren, aber 2 auch die Fehlerordnung O(h ). Algorithmus: (Runge-Kutta-Formel vierter Ordnung) Dies ist folgende Erweiterung der Simpson-Regel:
y 0 = u0 ,
h yk = yk−1 + (k1 + 2k2 + 2k3 + k4 ) , 6
wobei
k1 = f (tk−1 , yk−1 )
h h k2 = f (tk−1 + , yk−1 + k1 ) 2 2
h h k3 = f (tk−1 + , yk−1 + k2 ) 2 2
k4 = f (tk , yk−1 + hk3 )
Dieses Verfahren hat pro Schritt im wesentlichen den vierfachen Aufwand wie das explizite 4 Euler-Verfahren, dafür aber auch die Fehlerordnung O(h ).
79
8.4.3 Implementation Programm: (ode.sci) // / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / // /
ordinary
differential
equations
// / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / function
y = EE_step
( f , t , x , dt )
y = x+d t ∗ f ( t , x ) ; endfunction function
y_ = RK4_step
(f ,
t_n ,
y_n ,
h_n )
h2 = h_n / 2 ; k1 = f ( t_n , y_n ) ; k2 = f ( t_n+h2 ,
y_n+h2 ∗ k1 ) ;
k3 = f ( t_n+h2 ,
y_n+h2 ∗ k2 ) ; y_n+h_n ∗ k3 ) ;
k4 = f ( t_n+h_n ,
y_ = y_n + h_n / 6 ∗ ( k1 +2 ∗ ( k2+k3 )+k4 ) ; endfunction global
stepper ;
s t e p p e r = EE_step ; function
[x, t ]
global dim =
=
trajectory
( f , x0 , t 0 , t 1 , N)
stepper ; s i z e ( x0 , 1 ) ;
t = l i n s p a c e ( t 0 , t 1 , N) ; x = zeros x(: ,1) for
( dim , N) ;
= x0 ;
i =2:N
d t = t ( i )− t ( i
− 1) ;
x ( : , i )=s t e p p e r ( f , t ( i ) , x ( : , i
− 1) , d t ) ;
end endfunction function [x, t ]
plot_trajectory =
( f , x0 , t 0 , t 1 , N)
t r a j e c t o r y ( f , x0 , t 0 , t 1 , N) ;
plot (x (1 ,:) ,x (2 ,:) ) ; endfunction function [x, t ]
plot_paths =
( f , x0 , t 0 , t 1 , N)
t r a j e c t o r y ( f , x0 , t 0 , t 1 , N) ;
plot ( t , x) ; endfunction //
P r e d a t o r −p r e y
model
80
global mu =
mu ;
1.0;
function global y =
y = lotka_volterra
(t , x)
mu ;
[ x(1) x(2)
∗(1.0 − x ( 2 ) ) ; ∗ ( x ( 1 )−mu) ] ;
...
endfunction s t e p p e r = EE_step ; // p l o t _ p a t h s //
( lotka_volterra , [ 1 . 0 ; 1 . 0 ] , 0 , 2 0 , 1 0 0 )
Konvergenzuntersuchung
function
f = identity
anhand
e ^x
(t , x)
f = x; endfunction function for
errors
=
errors
(N)
i =2:N
x =
trajectory ( identity ,[1] ,0 ,1 ,2^ i ) ;
e r r o r s ( i ) = x ( 2 ^ i )−e x p ( 1 ) ; end endfunction // s t e p p e r=EE_step ; // e r r o r s ( 8 ) // s t e p p e r=RK4_step ; // e r r o r s ( 8 )
8.4.4 Anwendung Beobachtung: Bei Anwendung auf das Räuber-Beute-Modell mit
1 • 1
µ=1
sehen wir:
ist ein stationärer Punkt, das heiÿt er bleibt für alle Zeiten konstant. Wir
haben ein perfektes biologisches Gleichgewicht.
•
u0 nicht dieses Gleichgewicht ist, so beobachten wir Oszillationen, die umso stärker werden, je weiter u0 vom Gleichgewicht entfernt ist. Solche Oszillationen
Wenn der Anfangswert
beobachtet man auch in der Natur!
•
Das explizite Euler-Verfahren ist viel schlechter als das Runge-Kutta-Verfahren.
Beobachtung: Die Untersuchung des Konvergenzverhaltens von explizitem Euler-Verfahren und Runge-Kutta-Verfahren angewandt auf das Problem
u(t) ˙ = u(t) , t ∈ [0, 1] ,
81
u(0) = 1
mit exakter Lösung
u(t) = et
liefert folgende Tabelle für den Fehler
N
EE
RK-4
4
- 0.3479115
- 0.0002121
8
- 0.1717821
- 0.0000084
16
- 0.0854031
- 0.0000004
32
- 0.0425855
- 2.388E-08
64
- 0.0212644
- 1.419E-09
128
- 0.0106252
- 8.651E-11
256
- 0.0053109
- 5.341E-12
yN − e1 (h =
1 ): N
Wir sehen hier auch zahlenmäÿig die viel bessere Konvergenz des Runge-Kutta-Verfahrens vierter Ordnung.
8.4.5 Ausblicke
•
Die numerische Quadratur ist ein Spezialfall der numerischen Lösung von gewöhnlichen Dierentialgleichungen. Vieles dort im Ausblick erwähnte macht daher auch hier Sinn, insbesondere Extrapolation und adaptive Verfahren.
•
Es gibt gut gestellte Probleme (sogenannte steife Dierentialgleichungen), bei denen die hier vorgestellten expliziten Zeitschrittverfahren versagen (genauer: viel zu kleine Zeitschritte benötigen). Die Abhilfe sind hier implizite Zeitschrittverfahren, bei denen in jedem Schritt ein nichtlineares Gleichungssystem gelöst werden muss.
•
Eine weitaus komplexere Klasse von Problemen als die gewöhnlichen Dierentialgleichungen sind partielle Dierentialgleichungen, bei denen Ableitungen in mehrere Koordinatenrichtungen auftreten. Deren numerische Behandlung übersteigt den Rahmen dieser Vorlesung leider erheblich.
82
9 Übersetzungstabelle Englisch hat sich als Wissenschaftssprache etabliert. Somit ist es nützlich, auch die englischen Fachausdrücke zu kennen. Die folgende Tabelle übersetzt daher den Index dieses Skripts. äquidistante Unterteilungen
equidistant subdivisions
Stützstelle
support point
absolut konvergente Reihe
absolutely convergent series
absolute Kondition
absolute condition
adaptives Quadraturverfahren
adaptive quadrature scheme
Anfangsbedingung
initial condition
Approximation
approximation
arithmetische Operation
arithmetic operation
Assoziativitätsgesetz
associativity law
Auslöschung
cancelation
autonom
autonomous
Basis
base
blockweise
blockwise
Cache
cache
Cache-Optimierung
cache-optimisation
Cauchy-Folge
Cauchy sequence
CG-Verfahren
CG method
Cholesky-Zerlegung
Cholesky decompostion
dünnbesetzte Matrizen
sparse matrices
De nitheit
de niteness
Determinante
determinant
Dierentialgleichung
dierential equation
dierenzierbar
dierentiable
direkte Problem
direct problem
Diskretisierung
discretisation
Distributivgesetz
distributivity law
dividierten Dierenzen
divided dierences
doppelt-genaue IEEE-Gleitkommazahl
double precision oating point number
Drehungen
rotations
Dreiecksgestalt
triangular shape
83
Dreiecksungleichung
triangle inequality
euklidische Norm
Euclidean norm
Euklidische Skalarprodukt
Euclidean scalar product
exakt
exact
expliziten Zeitschrittverfahren
explicit time-stepping scheme
explizites Euler-Verfahren
explicit Euler method
Exponent
exponent
Extrapolation
extrapolation
Fehlerfortp anzung
error propagation
Feinheit
neness
Fixpunkt
xed point
Fixpunkt-Iteration
xed point iteration
Frobenius-Norm
Frobenius norm
Funktional
functional
Gauÿ-Quadratur
Gauss quadrature
Gauÿsches Eliminationsverfahren
Gauss elimination
gewöhnliche Dierentialgleichung
ordinary dierential equation
Gewicht
(quadrature) weight
Gitter
mesh
Gitterweite
meshwidth
Gleitkommazahlen
oating point number
globales Newton-Verfahren
global Newton method
GMRES-Verfahren
GMRES method
gut konditioniert
well-conditioned
Halbierungsverfahren
bisection method
Homogenität
homogeneity
Hutfunktion
hat function
implizite Zeitschrittverfahren
implicit time-stepping scheme
injektiv
injective
instabil
unstable
Integralkurven
integral curve
interpoliert
interpolated
Intervallschachtelung
nested intervals
Iteration
iteration
Iterationsvorschrift
iteration function
iteratives Verfahren
iterative method
Jacobi-Matrix
Jacobian
84
Knoten
node(s)
komponentenweise
component-wise
komponentenweisen Fehlerabschätzungen
component-wise error estimates
Kondition
condition
Kontinuum
continuum
Kontraktion
contraction
Kontraktionsrate
contraction rate
Konvergenzrate
convergence rate
konvex
convex
Krümmung
curvature
kubischer Spline
cubic spline
lösbar
solvable
Laufzeitfehler
runtime error
linear
linear
linear konvergent
linearly convergent
lineare Funktion
linear function
lineares Ausgleichsproblem
linear least-squares problem
linearen Operators
linear operator
linearer Splines
linear spline
lineares Gleichungssystem
linear system of equations
LR-Zerlegung
LU decomposition
Mantissenlänge
mantissa length
Matrix-Matrix-Multiplikation
matrix-matrix multiplication
matrixwertige
matrix-valued
Maximumsnorm
maximum norm
Modulo-Arithmetik
modulo arithmetic
Monombasis
monomial basis
natürliche Splineinterpolation
natural spline interpolation
Newton-Basis
Newton basis
Newton-Cotes-Formeln
Newton-Cotes formulas
Newton-Verfahren
Newton's method
Newton-Verfahren mit Dämpfungsstrategie
Newton method with
damping strategy Normalgleichung
normal equation
Normen
norms
numerisch
numerical
numerische Integration
numerical integration
numerischen Dierentiation
numerical dierentiation
obere Dreiecksgestalt
upper-triangular shape
Operatornorm
operator norm
Optimierungsproblem
optimisation problem
85
orthogonal
orthogonal
orthogonalen Zerlegungen
orthogonal decomposition
paarweise verschieden
paarwise dierent
partielle Dierentialgleichungen
partial dierential equations
partiellen Ableitungen
partial derivative
Permutationsmatrix
permutation matrix
Polygonzug
polygon
positiv
positive
Projektion
projection
QR-Zerlegung
QR decomposition
Quadratur
quadrature
quadratische Konvergenz
quadratic convergence
Quadraturformel
quadrature rule
Rückwärtssubstitution
backward substitution
regulär
regular
relative Kondition
relative condition
relativer Fehler
relative error
Restmatrix
rest matrix
Romberg-Quadratur
Romberg quadrature
Rundung
rounding
Schema von Neville
Neville's scheme
schlecht gestellt
ill-posed
schlecht konditioniert
ill-conditioned
Signum
signum
Simpson-Regel
Simpson's rule
Spaltenpivotsuche
column pivoting
Spektralnorm
spectral norm
Spektralradius
spectral radius
Spiegelung
re ection
Spline
spline
stückweise lineare Funktionen
piecewise linear function
Stützstellen
support point
stabil
stable
Stammfunktion
antiderivative
stark (strikt) diagonaldominant
strongly (strictly) diagonally dominant
stationärer Punkt
stationary point
steife Dierentialgleichung
sti dierential equation
stetig
continuous
stetig dierenzierbar
continuously dierentiable
Summennorm
sum norm
86
Supremumsnorm
supremum norm
surjektiv
surjective
symmetrisch
symmetric
symmetrisch positiv de nit
symmetric positive de nite
System gewöhnlicher Dierentialgleichungen
system of ordinary dierential equations
Trapezregel
trapecoidal rule
Unterteilung
subdivision
Überlauf
over ow
Vektorfeld
vector eld
Verfahren von Heun
method of Heun
vollständige Splineinterpolation
complete spline interpolation
Vorkonditionierer
preconditioner
Vorwärtssubstitution
forward substitution
wissenschaftliche Darstellung
scienti c notation
Zier
digit
zusammengesetzte Quadraturformel
composite quadrature rule
87
