Einführung in die Numerische Mathematik für Chemiker

Einführung in die Numerische
Mathematik für Chemiker
1. Auflage, Stand vom 8. November 2016
Bernd Hartke
Theoretische Chemie
Institut für Physikalische Chemie
Max-Eyth-Straße 2
Erdgeschoß, Raum 29
Tel.: 0431/880-2753
[email protected]
http://ravel.pctc.uni-kiel.de
“Sprechstunde” : jederzeit nach Vereinbarung
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Einleitung
Ziele der Veranstaltung
• erstes Grundverständnis von numerischer im Gegensatz zu analytischer
Mathematik
• Erstellung und Verwendung eigener Programme statt fertiger “apps”
• Verwendung von professionellen numerischen Bibliotheken
• Hintergrundwissen für vorgefertigte Programme der angewandten Mathematik und numerischen Simulation in Physik, Chemie,. . .
KEINE Ziele der Veranstaltung
• detaillierte Einführung in spezielle tools, insbesondere:
– TheoChem-Pakete (Gaussian, Turbomole, Molpro, Orca, . . . )
– algebraisch-numerische Mathe-Programme (Mathematica, Maple, . . . )
– Programmiersprachen (Fortran, C, C++, Java, Python, . . . )
• physikalisch- oder theoretisch-chemischer Stoff
2
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Veranstaltungsformat
• mehrere thematische Kapitel; zu jedem Kapitel gibt es:
– eine Vorlesungsstunde: Einführung in die theoretischen Grundlagen;
Besprechung des zugehörigen Aufgabenblatts
– eine betreute Übungsstunde
• alle Skriptseiten und alle Aufgabenblätter von meiner homepage als
pdf-Dateien herunterladbar:
http://ravel.pctc.uni-kiel.de → teaching
→ Numerical Mathematics for Chemists
(Achtung: Inhalt wird sich während des Semesters fortlaufend ändern)
• die Benotung erfolgt über die Bearbeitung der Übungsaufgaben:
– abzugeben sind:
∗ mindestens: lauffähiger Programmtext, in Octave/Matlab oder
Fortran
∗ keine ausführbaren Programme
∗ ggf. (als pdf-Datei): output, Antworten auf Fragen, weitere Erläuterungen (z.B. bei Schwierigkeiten), etc.
– Funktionierende Programme sind anzustreben.
In Ausnahmefällen reicht es, wenn für mich aus dem Programmtext
klar erkennbar ist, daß die Prinzipien verstanden wurden.
– Abgabe der Aufgaben am besten per e-mail:
[email protected]
[email protected]
3
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Inhaltsplan
1. Einführung in grundlegende Programmierkonzepte
(Octave/Matlab, Fortran)
2. Zahlendarstellung, Fehler, Konditionszahl
3. Differentiation
4. Integration: Trapezregel, Simpson, Gauß;
Anwendung: überall
5. Zufallszahlen und Monte-Carlo-(MC)-Integration
Anwendung: MC-Simulationen von Gasen, Flüssigkeiten,. . .
6. Suche nach Nullstellen und Extremwerten: Newton et al.
Anwendung: Lösung analytisch unlösbarer Gleichungen; zahlreiche Optimierungsprobleme
7. gewöhnliche Differentialgleichungen (Euler, Runge-Kutta)
Anwendung: z.B. einzelne Gleichungen der Kinetik
8. DGL-Systeme und partielle DGLs
Anwendung: komplizierte Kinetiken, klassische Mechanik, Quantenmechanik, Wärmeleitung, Elektrodynamik, Akustik, Strömungsdynamik,. . .
9. Lösung linearer Gleichungssysteme (inkl. Matrixinversion);
Anwendung: überall, z.B. Regression ( Fit“) und Interpolation
”
10. Matrixdiagonalisierung: Eigenwertproblem;
Anwendung: Quantenmechanik, stationäre Lösungen der Schrödingergleichung
11. lineare Regression, χ2 -Statistik; nichtlineare Regression (Levenberg-Marquardt);
Spline-Interpolation;
Anwendung: z.B. Anpassung von Modellfunktionen an Meßdaten
Vorwissen-Umfrage
4
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Literaturempfehlungen:
• W. H. Press, B. P. Flannery, S. A. Teukolsky und W. T. Vetterling:
Numerical Recipes“, Cambridge University Press, Cambridge, 1990:
”
sehr gut geschriebene Einführung in die Grundalgorithmen der Numerischen Mathematik, mit Beispielprogrammen; weit verbreitetes Standardwerk.
• Roland W. Freund und Ronald W. Hoppe: “Stoer/Bulirsch: Numerische
Mathematik 1 und 2”, Springer, Berlin, 10. Auflage, 2007:
bewährtes Standardwerk deutscher Mathe-Fakultäten zur Numerischen
Mathematik; klare Darstellung, aber etwas “mathematischerer Stil”.
• M. Hermann: “Numerische Mathematik”, Oldenburg, 2006:
gut zugängliches Buch, leider ohne Differentialgleichungen und Optimierung.
• Peter Deuflhard und Andreas Hohmann: “Numerische Mathematik 1,2,3:
Eine algorithmisch orientierte Einführung”, deGruyter, Berlin, 4. Auflage,
2008:
Band 1 bietet eine allgemeine Einführung in diverse Themen, Band 2/3
beschäftigen sich mit Differentialgleichungen; klarer Text, leider keine
Behandlung von Optimierung.
• A. Neumaier: Introduction to Numerical Analysis“, Cambridge Univer”
sity Press, Cambridge, 2001:
gutes Buch eines Genies in der Numerischen Mathematik.
• H. R. Schwarz und J. Waldvogel (Hrsg.): Numerische Mathematik“,
”
Teubner-Verlag, 1997:
umfassender Klassiker, aber etwas schwerer zugänglich.
• G. H. Golub und C. F. van Loan: Matrix Computations“, Johns Hop”
kins University Press, 1989, 1993:
Klassiker zur linearen Algebra; ziemlich mathematischer“ Stil; für De”
tailinformationen zu linearen Gleichungssystemen und zum Eigenwertproblem.
• P. E. Gill, W. Murray and M. H. Wright: “Practical Optimization”,
Academic Press, 1981:
Spezialbuch zur Optimierung mit vielen Praxistips; Herleitungen leider
nicht immer selbsterklärend.
5
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
1
Minimaleinführung ins Programmieren
1.1
Fortran
siehe separates Fortran-Einführungsskript
1.2
Octave/Matlab
siehe Octave-Einführung in Henrik Larssons Numerikskript WS15/16
1.42
Welche Programmiersprache soll ich nehmen?
Das ist weitestgehend egal, u.a.:
• alle Numerikalgorithmen sind programmiersprachenunabhängig;
• die grundlegenden Sprachstrukturen (Felder, Schleifen, Verzweigungen)
gibt es überall in sehr ähnlicher Form.
Andere Entscheidungshilfen:
• kleines Hilfsprogramm vs. großes software-Projekt
• eigene Entwicklung für Eigenbedarf vs. kollaboratives multi-user-Paket
• Syntax:
– langatmig aber gut lesbar vs. kurz aber kryptisch
?
– frei = versteckte Komplexität
• interpretiert vs. compiliert
• virtual machine vs. OS/hardware-spezifische compiler/interpreter
• Verfügbarkeit ausgefeilter Entwicklungsumgebungen (Editor mit Syntaxhighlighting und Inhaltshilfen, integrierte debugger/profiler, integrierte
Versionsverwaltung, . . . )
• Verfügbarkeit professioneller Bibliotheken (BLAS, LAPACK, MKL, etc.;
daran hängt letztlich die Laufzeiteffizienz und wieviel code man selbst
schreiben muß)
6
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
2
Elementare Grundlagen
2.1
Zahlendarstellung
Grundproblem:
“So gut wie alle” reellen Zahlen haben unendlich viele signifikante Stellen
(egal in welcher Basis). ABER: Der Computerspeicher ist endlich.
Grundsätzliche Folgen für die numerische Praxis:
• begrenzte Anzahl Stellen
• es gibt eine betragsmäßig größte darstellbare Zahl
• miteinander verrechnete Zahlen sollten ähnlich groß sein
Zahlendarstellung im Computer ist binär, aber mit weiteren technischen
Tricks. Alle Details: siehe Skript WS15/16 (Henrik Larsson)
2.1.1
Ganze Zahlen (integers)
# bits
8
16
32
64
8
16
32
64
Vorzeichen octave-Name
unsigned
unsigned
unsigned
unsigned
signed
signed
signed
signed
uint8
uint16
uint32
uint64
int8
int16
int32
int64
min
max
0
255
0
65 535
0
4 294 967 295
0 18 446 744 073 709 551 615
−128
127
−32 768
32 767
−2 147 483 648
2 147 483 647
−9 223 372 036 854 775 808 9 223 372 036 854 775 807
Durch die Zweierkomplementdarstellung (siehe Skript WS15/16 Henrik Larsson) sind vorzeichenbehaftete ganze Zahlen “seltsam angeordnet”:
Binärwert
8-bit-Zweierkomplement
00000000
00000001
···
01111110
01111111
10000000
10000001
···
11111110
11111111
0
1
···
126
127
-128
-127
···
-2
-1
7
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Was bei “Überlauf” des darstellbaren Zahlenbereichs passiert, hängt u.a. ab
von der hardware, der Programmiersprache und ggf. den Compileroptionen:
• ohne geeignete Compileroptionen verhalten sich Sprachen wie C, Fortran
und viele andere so wie dieses C-Programm:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#include<stdio.h>
int main(){
signed char a; //8 bit lang, entspricht int8 in octave
unsigned char b; //entspricht uint8 in octave
a = 127;
b = 255;
printf (”a = %d\n”,a);
printf (”b = %d\n”,b);
a = a + 1;
b = b + 1;
printf (”a + 1 = %d\n”,a);
printf (”b + 1 = %d\n”,b);
a = a + 999;
b = b + 999;
printf (”a + 1000 = %d\n”,a);
printf (”b + 1000 = %d\n”,b);
}
führt nach der Übersetzung in Maschinensprache zur Ausgabe von
a
b
a
b
a
b
=
=
+
+
+
+
127
255
1 = -128
1 = 0
1000 = 103
1000 = 231
Ob und wie das detektiert bzw. verhindert werden kann, hängt sehr
stark von der Programmiersprache und dem Compiler ab. Schwierig ist
auch, daß es nach dem Überlauf i.d.R. für Vorsichtsmaßnahmen zu spät
ist. Machbar und portabel können aber schon einfache Lösungen wie
diese sein:
– Tests der Art if (a+1 > a) ...
– temporäre Verwendung von real-Variablen
(s.u.; Vorsicht: ggf. Genauigkeitsverlust!)
– temporäre Verwendung von integer-Variablen mit mehr bits.
8
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• das Verhalten von octave/matlab ist anders:
– keine Warnung bei Überlauf
– stattdessen wird der größte/kleinste mögliche Wert verwendet
1
2
3
4
5
6
7
8
9
>> a = int8(127) % ganze Zahl mit einer Darstellung von 8 Bits
a = 127
>> a = a + 1 %addiere
a = 127
>> a = int8(−128);
>> a − 1
ans = −128
>> a + 20000
ans = 127
Ob das “besser” ist, ist Geschmackssache. . .
Da ganze Zahlen (direkt oder indirekt) häufig als “Zeiger” auf Speicheradressen
verwendet werden, ist integer-overflow ein Einfallstor für Computerviren.
Selbst Profis vergessen gelegentlich die Endlichkeit der integers:
• 2014 war die Zahl der Klicks des Youtube-Videos “Gangnam Style” auf
einmal negativ, da die Zahl zu groß für die verwendete signed-integerDarstellung geworden war.
• 1996 stürzte eine unbemannte Ariane5-Rakete ab, da eine Gleitkommazahl in eine ganze Zahl konvertiert wurde und dann jenseits des
ganzzahlig darstellbaren Bereichs lag.
2.1.2
Reelle Zahlen (reals)
Binäre Darstellung einer reellen Zahl:
(−1)v · M · 2e
(1)
mit Vorzeichenbit v, m-stelliger Mantisse M und Exponent e.
• Für x 6= 0 wird e so gewählt, dass das erste Bit in der Mantisse ungleich
Null ist (normalisierte Darstellung); dann muß es nicht gespeichert
werden (implizites erstes Bit)
• Ist die darzustellende Zahl zu klein für die normalisierte Form (e < emin ),
so kommt es zum Exponentenunterlauf (underflow) in den subnormalen
Zahlenbereich. Dort sinkt die Darstellungsgenauigkeit stetig, bis alle
Ziffern 0 sind.
9
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• spezielle Zeichen/Formen:
Exponent
normalisierte Form
denormalisierte Form
±0
±inf
NaN
Mantisse
M ∈ [2−1 , 1 − 2−m ]
M ∈ [2−m , 2−1 − 2−m ]
0
0
6= 0
e ∈ [emin + 1, emax − 1]
e = emin + 1
e = emin
e = emax
e = emax + 1
• NaN (not a number ) ergibt sich u.a. bei Division durch Null
• ±inf tritt bei einem Exponentenüberlauf (overflow) auf. Die Zahl ist
zu groß, um dargestellt zu werden.
• Die meisten reellen Zahlen lassen sich nicht exakt darstellen:
1
2
3
4
>>
ans
>>
ans
printf(’%20.17f\n’,0.3) %0.3 kann nicht exakt wiedergegeben werden
= 0.29999999999999999
printf(’%20.17f\n’,0.5) %genaue Darstellung, da 0.5 = 1/2
= 0.50000000000000000
• Dies führt zu entsprechenden Rundungsfehlern:
1
2
3
4
5
6
7
8
>> if 1.0 − 0.8 − 0.2 == 0
disp( ’ 1.0 − 0.8 − 0.2 ist genau Null’);
else
disp( ’ 1.0 − 0.8 − 0.2 ist NICHT genau Null’);
end
1.0 − 0.8 − 0.2 ist NICHT genau Null
>> printf(’%20.17f\n’,1.0 − 0.8 − 0.2);
−0.00000000000000006
Deshalb sollte man bei real-Zahlen nie auf exakte Gleichheit testen:
1
2
3
4
5
6
>> if 1.0 −
disp( ’ 1.0
else
disp( ’ 1.0
end
1.0 − 0.8 −
0.8 − 0.2 < 1.0e−8
− 0.8 − 0.2 ist ziemlich genau Null’);
− 0.8 − 0.2 ist deutlich verschieden von Null’);
0.2 ist ziemlich genau Null
• Aufgrund der halbexponentiellen Darstellung Gl. 1 sind die exakt darstellbaren Zahlen logarithmisch auf der Zahlenskala verteilt. Beispiel für
M bestehend aus drei Bits und e ∈ [−1, 2]:
0
5 7
16 16
1 3 1
4 8 2
5
8
7
8
3
4
7
4
1
5
2
2
3
2
10
7
2
3
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• Maschinengenauigkeit (“Maschinenepsilon”): maximaler relativer Abstand
zweier benachbarter, normalisierter real-Zahlen (m Mantissenstellen,
Unterschied von 1 in der letzten Stelle):
2−m
= 21−m = 2
−1
2
(2)
• Skript berechne_eps.m zum Berechnen von :
1
2
3
4
5
6
7
8
disp( ’Berechne das Maschinenepsilon’);
meps = 1.0;
while 1.0 + 0.5 ∗ meps > 1.0
meps = meps ∗ 0.5;
end
fprintf ( ’Berechnetes Maschinenepsilon:
%20.17e\n’,meps);
fprintf ( ’Analytisch berechnetes Maschinenepsilon: %20.17e\n’,2 ∗ 2ˆ(−53)); %
Mantisse: 52 bits + implizites erste bit = 53
fprintf ( ’ Intern verwendetes Maschinenepsilon: %20.17e\n’,eps); %eps ist eine
Konstante in octave
• Ausgabe:
1
2
3
4
5
>> berechne eps
Berechne das Maschinenepsilon
Berechnetes Maschinenepsilon:
2.22044604925031308e−16
Analytisch berechnetes Maschinenepsilon: 2.22044604925031308e−16
Intern verwendetes Maschinenepsilon: 2.22044604925031308e−16
Gebräuchliche Darstellungen in den meisten Sprachen:
• single precision:
– 32 bit (4 byte):
23(+1) bits Mantisse, 8 bits Exp. (e ∈ [−126, 129]), 1 bit Vorzeichen
– Wertebereich der normalisierten Zahlen: 1.17 · 10−38 − 3.4 · 1038
– Wertebereich der denormalisierten Zahlen: 1.4 · 10−45 − 1.17 · 10−38
– Genauigkeit: = 2−24 = 6 · 10−8
• double precision:
– 64 bit (8 byte):
52(+1) bits Mantisse, 11 bits Exp. (e ∈ [−1022, 1025]), 1 bit Vorz.
– Wertebereich der normalisierten Zahlen: 2.22 · 10−308 − 1.797 · 10308
– Wertebereich der denormalisierten Zahlen: 3.9 · 10−324 − 2.22 · 10−308
– Genauigkeit: = 2−53 = 1.1 · 10−16
11
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Auch bei reals sind schon gravierende Praxisfehler aufgetreten:
• 1992 wurde bei der Schleswig-Holsteiner Landtagswahl aufgrund eines
Rundungsfehlers fälschlicherweise angenommen, daß die Partei Die Grünen die 5 %-Hürde überschritten hat.
• Im ersten Golfkrieg hat eine Patriot-Rakete 28 US-Soldaten getötet, da
die Zeit in Inkrementen von 1/10 gezählt wurde. 1/10 hat allerdings keine
endliche binäre Zahlendarstellung, sodaß Rundungsfehler auftraten.
2.2
Fehlerarten
Es gibt verschiedene Arten von Fehlern:
Modellfehler
• Fehler aufgrund des Unterschiedes zwischen dem realen System und dem
mathematischen Modell
• Beispiele:
– Pendel mit Luftreibung, Corioliskraft, etc. versus Newtonsche Gleichung im Inertialsystem
– Molekül in Lösung versus Molekulardynamiksimulation in der Gasphase mit klassischer Beschreibung der Atomkerndynamik und approximativer quantenmechanischer Beschreibung der Elektronen
Datenfehler
• Abweichung der realen Parameterwerte von den eingegebenen Daten
• Beispiele:
– Signifikante Stellen von Naturkonstanten
– Messungenauigkeiten im Experiment
Verfahrensfehler
• Intrinsischer Fehler des numerischen Verfahrens
• Beispiele:
– Endliche Taylorreihenentwicklung
– Diskretisierungsfehler
12
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• In diesem Zusammenhang heißt ein Verfahren konvergent, wenn der
Verfahrensfehler gegen null gehen kann
• Beispiel: Verfeinerung der Diskretisierung, z.B. bei einer Integration
• Gegenbeispiel: Divergente Reihenentwicklung, z.B. bei quantenmechanischer Vielteilchenstörungstheorie
Rundungsfehler
• Fehler aufgrund der Darstellung von Gleitkommazahlen auf dem Computer
• Beispiele
– Das Distributivgesetz gilt nicht: a × (b + c) 6= a × b + a × c
– Das Assoziativgesetz gilt nicht: a + (b + c) 6= (a + b) + c
Andere nützliche Unterscheidungen:
f (x) sei das mathematische Problem in Abhängigkeit von der Eingabe x, f˜
der numerische Algorithmus und x̃ fehlerbehaftete Eingabedaten:
Kondition: ||f (x̃) − f (x)||: Wie stark schwankt das Problem bei Störung?
Stabilität: ||f˜(x̃) − f˜(x)||: Wie stark schwankt der numerische Algorithmus
bei Störung?
Konsistenz: ||f˜(x) − f (x)||: Wie gut löst der numerische Algorithmus (mit
exakter Eingabe) das Problem?
Konvergenz: ||f˜(x̃) − f (x)||: Wie gut löst der numerische Algorithmus (mit
gestörter Eingabe) das Problem?
13
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
2.3
Normen
2.3.1
Vektoren
Die Norm k · k in einem Vektorraum X erfüllt die Bedingungen
1. Definitheit: k~xk ≥ 0, ~x ∈ X; insbes. k~xk = 0 ⇒ ~x = ~0
2. Homogenität: kα~xk = |α|k~xk, α ∈ C
3. Dreiecksungleichung, Subadditivität: k~x + ~y k ≤ k~xk + k~y k
Dies ist für zahlreiche unterschiedliche Normvorschriften gegeben, für X = Rn
(gewöhnliche Vektoren) z.B.:
1-Norm k~xk1 =
Pn
i
|xi |
Euklidische Norm k~xk2 =
pP
p-Norm k~xkp = p i |xi |p
pP
|xi |2
i
Maximums-Norm k~xk∞ = maxi |xi |
Diese Vielfalt stört i.d.R. nicht, weil sich alle Normen gegeneinander abschätzen lassen:
C1 k~xkα ≤ k~xkβ ≤ C2 k~xkα
(3)
Die Konstanten C1 , C2 lassen sich in den meisten Fällen genau angeben.
2.3.2
Matrizen, Abbildungen
Matrixnormen können über Vektornormen definiert werden oder können auch
andere Charakteristika von Matrizen involvieren, die bei Vektoren nicht
auftreten (z.B. Eigenwerte, Singulärwerte).
Sinnvollerweise sollten Matrixnormen k · k mit Vektornormen k · kV verträglich
(kompatibel) sein; das ist dann der Fall, wenn
kA~xkV ≤ kAk k~xkV
(4)
Normen und Eigenwerte hängen jedoch auch zusammen: Ist eine Matrixnorm
mit irgendeiner Vektornorm verträglich, dann gilt für jeden Eigenwert λ einer
quadratischen Matrix A
|λ| ≤ kAk
(5)
weil mit dem zugehörigen Eigenvektor ~x gezeigt werden kann:
|λ| k~xk = kλ~xk = kA~xk ≤ kAk k~xk
14
(6)
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
woraus die Behauptung nach Division durch k~xk folgt. Nach Gl. 5 ist der
Spektralradius (Betrag des betragsgrößten Eigenwerts) einer Matrix niemals
größer als ihre Norm.
Beispiele für Matrixnormen:
• basierend auf der Operatornorm:
kA~xk
kAk = max
= max kA~xk
k~xk
k~xk=1
~x6=~0
1-Norm kAk1 = maxj∈[1,n]
(7)
Pn
|Aij | (Spaltensummennorm)
√
Euklidische (Spektral-)Norm kAk2 = λmax ; λmax ist der größte Eigenwert von
A† A (größter Singulärwert von A, siehe SVD)
P
inf-Norm kAk∞ = maxi∈[1,n] nj=1 |Aij | (Zeilensummennorm)
i=1
• basierend auf einer vektorartigen Anneinanderreihung aller Matrixelemente:
Frobeniusnorm kAkF =
qP P
m
n
i
j
|Aij |2
Auch Matrixnormen lassen sich nach Gl. 3 gegeneinander abschätzen.
15
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
2.4
Kondition
Aufgrund der finiten Zahlendarstellung auf dem Computer sollte man nicht
von einer exakten Eingabe x ausgehen, die zu einem exakten Resultat f (x),
führt, sondern von einer Eingabemenge E, welche zu einer Resultatsmenge
R = f (E) führt:
E
x
f (x)
R
Die Kondition charakterisiert das Verhältnis zwischen E und R.
• Sei ~y = F (~x) unter der Annahme von exakten Zahlen (ohne Rundungsfehler), sowie
~ eine Zahl mit Rundungsfehler ist.
~y˜ = F (~x˜), wobei ~x˜ = ~x + δx
• Der absolute bzw. relative Fehler sind dann
k~y˜ − ~y k
k~y k
(8)
κabs = inf{α | k~y˜ − ~y k ≤ αk~x˜ − ~xk},
(9)
∆abs = k~y˜ − ~y k ,
∆rel =
• Die absolute Konditionszahl ist
also das kleinstmögliche Verhältnis k~y˜ − ~y k × (k~x˜ − ~xk)−1 .
• Sofern F (~x) differenzierbar ist, gilt
κabs = k∇F (~x)k ≡ kF 0 (~x)k.
(10)
(hier wird die Matrixnorm gebraucht)
• Die relative Konditionszahl ist analog
)
( k~y˜ − ~y k
˜
k
~
x
−
~
x
k
≤α
,
κrel = inf α k~y k
k~xk
das kleinstmögliche Verhältnis
k~
y˜−~
y k/k~
yk
(11)
× (k~x˜−~xk/k~xk)−1 .
• Sofern F (~x) differenzierbar ist, gilt
κrel =
kF 0 (~x)k
k~xk
× k~xk = κabs ×
.
kF (~x)k
kF (~x)k
Im Folgenden: Beschränkung auf κ ≡ κrel
• Ist κ - 1, so ist das Problem gut konditioniert.
– Für κ < 1 gibt es sogar eine Fehlerverringerung
– Für κ = 1 treten nur Rundungsfehler auf
16
(12)
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• Für κ = O(100) ist das Ergebnis typischerweise akzeptabel
• Für κ 104 ist das Ergebnis oft unbrauchbar. Das Problem ist schlecht konditioniert
• Falls das Ergebnis eines numerischen Problems unstetig von den Eingabedaten
abhängt, wird es als unsachgemäß gestellt angesehen.
Die Kondition eines Problems hängt nicht nur vom Problem, sondern auch
von der Eingabemenge E ab.
Beispiele:
• Schnittpunktsberechnung zweier Geraden:
– gut konditioniert:
g
g̃
r̃
r
h
1
h̃
Der Schnittpunkt r̃ variiert hier ähnlich stark wie die Geraden g
und g̃ selbst.
– schlecht konditioniert (schleifender Schnitt):
r
Starke Schnittpunktsvariation
• Addition/Subtraktion:
– Fadd (a, b) = a + b bzw. Fsub (a, b) = a − b
– mit Eingabevektor ~x = (a, b)T ; wähle 1-Norm: k(a, b)T k1 = |a| + |b|
– κabs = kF 0 (a, b)k1 = k(1, 1)T k1 = 1. (Matrixnorm!)
– κrel,add =
|a|+|b|
|a+b| ,
κrel,sub =
|a|+|b|
|a−b|
⇒ Die Addition von Zahlen mit gleichem Vorzeichen ist gut konditioniert:
κrel = 1
Aber: Die Subtraktion ist schlecht konditioniert, wenn
|a + b| |a| + |b|,
also wenn |a| und |b| annähernd gleich sind.
17
(13)
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
Zahlenbeispiel:
Rechnung mit 5 Dezimalstellen:
x = 1.2345 · 10−3
y = 1.7899 · 10−3
u = 1.0000 · 100 + x = 1.0012 · 100
v = 1.0000 · 100 + y = 1.0018 · 100
z1 = v − u = 6.000 · 10−4
z2 = y − x = 5.554 · 10−4
Mathematisch ist z1 = z2 . Numerisch kommt es zum Weghebephänomen
und ein erheblicher Rundungsfehler tritt auf.
Subtraktion von annähernd gleichen Zahlen
sollte vermieden werden!
• weitere Beispiele siehe Henrik Larssons Numerikskript WS15/16.
18
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
3
Differentiation
Die erste Ableitung einer Funktion f (x) nach x an der Stelle x0 ist definiert
als
f
(x)
−
f
(x
)
df
∆f
0
f 0 (x)|x=x0 =
=
lim
(14)
=
lim
∆x→0 ∆x x=x
dx x=x0 x→x0
x − x0
0
Zu dem analytischen Problem, daß es sich hier um einen unbestimmten
Ausdruck vom Typ “ 00 ” handelt, kommt aus numerischer Sicht das Problem
hinzu, daß mit kleiner werdendem ∆x die Beträge der Zahlenpaare x, x0 und
f (x), f (x0 ) immer ähnlicher werden, während sich die Größenordnung dieser
vier Zahlen selbst nicht ändert.
Algorithmenunabhängige Folgerungen:
• Numerische Differentiation verstößt zwangsläufig gegen die Regel, Subtraktion annähernd gleicher Zahlen zu vermeiden!
• Der Limes in Gl. 14 ist numerisch (bei begrenzter Zahlendarstellung)
nicht ausführbar.
Die abgebrochene Taylorreihe behebt diese Probleme nicht, ermöglicht aber
• Approximationen an Ableitungen, ohne Ausführung des Limes
(finite Differenzen),
• und mit einer Abschätzung des Fehlers.
3.1
Vorwärtsdifferenzen
Taylorentwicklung von f (x) um x0 (mit ∆x = h):
h2 00
f (x0 + h) = f (x0 ) + hf (x0 ) + f (x0 ) + O h3 f (3) (x0 )
2
f (x0 + h) − f (x0 )
1
⇒
= f 0 (x0 ) + hf 00 (x0 ) + O(h2 f (3) (x0 ))
h
2
0
(15)
(16)
Also gilt
f (x0 + h) − f (x0 )
(17)
h
mit einem Fehler linear in h (also vgl.weise schlechte Konvergenz: halb so
großes h halbiert den Fehler).
f 0 (x0 ) ≈
19
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
3.2
Zentraldifferenzen
h2 00
h3 000
f (xo + h) = f (x0 ) + hf (x0 ) + f (x0 ) + f (x0 ) + · · ·
2
3!
2
h3
h
f (xo − h) = f (x0 ) − hf 0 (x0 ) + f 00 (x0 ) − f 000 (x0 ) + · · ·
2
3!
2
f (xo + h) − f (xo − h) h 000
⇒ f 0 (x0 ) =
− f (x0 ) + · · ·
2h
3!
0
(18)
(19)
(20)
Fehler nur noch quadratisch in h (halb so großes h ⇒ ein Viertel des Fehlers).
Veranschaulichung:
m2
m1
f
m3
x̄ − h
x̄
x̄ + h
Unterschied zwischen Vorwärts- (gepunktete Linie), Rückwärts(durchgezogene Linie) und Zentraldifferenzen (gestrichelte Linie).
3.3
Höhere Ordnungen
Auf ähnliche Weise sind herleitbar:
• Formeln höherer Ordnungen für die erste Ableitung,
• Formeln unterschiedlicher Ordnungen für höhere Ableitungen.
Siehe z.B. https://en.wikipedia.org/wiki/Finite_difference_coefficient.
Alternativer Weg zu Formeln höherer Ordnung: Richardson-Extrapolation
(eine allgemeinere Idee zur Genauigkeitsverbesserung, wenn man zwei Approximationsformeln mit zwei unterschiedlichen h-Werten hat)
• Gl. 20 kann umgeschrieben werden zu
I : f 0 (x0 ) = D(h) + e2 h2 + e4 h4 + . . .
(21)
D(h) ist der Differenzenquotient und ei die von h unabhängigen Fehleranteile.
20
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
• Idee der Extrapolation: Berechne f 0 (x0 ) für Vielfache von h:
II : f 0 (x0 ) = D(2h) + 4e2 h2 + 16e4 h4 + . . .
(22)
4 × I − II
4D(h) − D(2h)
: f 0 (x0 ) =
− 4e4 h4 + . . .
(23)
3
3
−f (x0 + 2h) + 8f (x + h) − 8f (x − h) + f (x − 2h)
⇔ f 0 (x0 ) =
+ O(h4 )
12h
(24)
Der e2 -term fällt weg (ähnlich wie oben bei Zentraldifferenzen), dadurch
steigt die Fehlerordnung nochmals. Diese Methode kann automatisiert werden:
Beliebige Genauigkeit und Fehlerabschätzung möglich.
3.4
Wahl der Schrittweite
Achtung: h und/oder x0 + h ggf. nicht exakt binär darstellbar. Abhilfe:
1
2
temp = x0 + h;
h = temp − x0;
Das erhaltene h ist nun “exakt” in der Binärdarstellung (aufgrund des Fehlers
in der Subtraktion).
h kann sowohl zu groß als auch zu klein gewählt werden. Beispiel Vorwärtsdifferenzen:
• Der Rundungsfehler ist ungefähr ( = Maschinenepsilon)
f (x0 ) r ∼ h (25)
(Der Rundungsfehler ist generell proportional zu 1/h, auch für andere
finite Differenzenformeln.)
• Fehler aufgrund von Trunkierung der Taylorreihe für Vorwärtsdifferenzen:
t ∼ |hf 00 (x)|
Beachte: r ∼
1
h
(26)
aber t ∼ h !
• Minimierung von (h) = r + t liefert
s
f (x0 )
h∼
f 00 (x0 )
(27)
Die optimale Schrittweite hängt von der Maschinengenauigkeit und dem
00
Inversen der relativen Krümmung ff (x(x00)) ab.
21
Numerik für Chemiker • Prof. Dr. B. Hartke, Universität Kiel, [email protected]
3.5
Schlußbemerkungen
• Höhere Ordnung bewirkt nicht immer höhere Genauigkeit.
• Höhere Ordnung bedeutet nicht unbedingt bessere Effizienz.
• Funktionsberechnungen können teuer sein.
Beispiel: Gradienten/Frequenzen in der ab-initio-Quantenchemie.
Goldene Praxisregel: keine finiten Differenzen verwenden, sonderen analytische Ableitungen. Diese müssen jedoch zusätzlich einprogrammiert
werden. Das wird trotzdem gemacht ⇒ das lohnt sich! Beispiel: Anilin:
– eine SCF-Energieberechnung dauere 5 min
– bei 14∗3=42 Freiheitsgraden und Zentraldifferenzen (2 Punkte pro
Freiheitsgrad) braucht der komplette Gradientenvektor
42∗2∗5 min = 420 min = 7 h
– ⇒ eine konvergierte Geometrieoptimierung braucht 1 Tag!
– mit “analytischen Gradienten” braucht der komplette Gradientenvektor nur ca. 5–10 min und die konvergierte Geometrieoptimierung
ca. 1 h.
• Es existieren weitere Alternativen zur Ableitungsberechnung, z.B. via
Fouriertransformation:
dn
1
√
f
(x)
=
dxn
2π
Z∞
−∞
dn ikx
1
F(k) n e dk = √
dx
2π
Z∞
F(k) (ik)n eikx dk = f (n) (x)
−∞
(28)
22