Einführung in die numerische Mathematik für Chemiker

Mathematik für Chemiker

Bernd Hartke

Henrik R. Larsson

Theoretische Chemie

Institut für Physikalische Chemie Max-Eyth-Straße 2

Erdgeschoss, Raum 29 Tel.: 0431/880-2753

hartke@pctc.uni-kiel.de larsson@pctc.uni-kiel.de http://ravel.pctc.uni-kiel.de

„Sprechstunde“: Jederzeit nach Vereinbarung

Inhaltsverzeichnis

0 Vorbemerkungen 4

0.1 Inhaltsplan . . . 6

1 Minimaleinführung in octave 9 1.1 Installation . . . 10

1.2 Literatur . . . 10

1.3 Weitere Bibliotheken . . . 11

1.4 Ein erster Einstieg . . . 11

1.5 Programmierdisziplin! . . . 23

2 Zahlendarstellung 24 2.1 Natürliche Zahlen . . . 24

2.2 Ganze Zahlen . . . 24

2.3 Zahlenbereich und Über-/Unterschreitung . . . 25

2.4 Gleitkommazahlen . . . 26

3 Fehler und Kondition 30 3.1 Fehlerarten . . . 30

3.2 Einschub: Norm . . . 31

3.3 Konditionierung mathematischer Probleme . . . 32

4 Numerische Differentiation 36 4.1 Vorwärtsdifferenzen . . . 36

4.2 Zentrale Differenzen . . . 36

4.3 Höhere Ordnungen: Richardsons Extrapolation . . . 37

4.4 Wahl der Schrittweite . . . 38

4.5 Anmerkungen . . . 38

5 Integration 39 5.1 Vorbemerkungen: . . . 39

5.2 Trapezregel (Euler-Verfahren) . . . 40

5.3 Simpsonregel und andere Verfahren . . . 41

5.4 Adaptive Quadratur . . . 42

5.5 Einschub: Orthogonale Polynome . . . 43

5.6 Gauß-Integration . . . 44

5.7 Berechnungen der Wichtungen und Punkte . . . 46

5.8 Schwierige Integranden . . . 47

5.9 Mehrdimensionale Integration . . . 48

5.10 Monte-Carlo Verfahren . . . 49

6 Suche nach Nullstellen und Extremwerten 52 6.1 Nullstellensuche (root finding) . . . 52

6.2 Verfahren im Überblick . . . 60

6.3 Extremwertsuche = Optimierung . . . 61

6.4 Multidimensionale Minimierung . . . 66

6.5 Mehrdimensionale nichtlineare Nullstellensuche . . . 74

7 Gewöhnliche Differentialgleichungen 76

7.1 Einführung . . . 76

7.2 Explizites Euler-Cauchysches Polygonzugverfahren . . . 78

7.3 Runge-Kutta-Verfahren . . . 79

7.4 Adaptive Schrittweitenkontrolle . . . 82

7.5 Weitere Verfahren . . . 84

7.6 Steife Anfangswertprobleme . . . 86

7.7 Zusammenfassung und Praxistipps . . . 89

7.8 Partielle Differentialgleichungen . . . 91

8 Lineare Gleichungssysteme 95 8.1 Grundlegendes . . . 95

8.2 Konditionierung . . . 98

8.3 Gaußsches Eliminationsverfahren . . . 100

8.4 LU-Zerlegung . . . 103

8.5 Cholesky-Zerlegung . . . 105

8.6 Nachiteration . . . 107

8.7 Kahan-Summation . . . 108

8.8 Direkte Verfahren . . . 109

9 Matrixeigenwertproblem 112 9.1 Grundlagen . . . 112

9.2 Orthogonale Matrizen in zwei Dimensionen . . . 114

9.3 Jacobi-Verfahren . . . 116

9.4 Andere Verfahren . . . 120

9.5 Eigenwertabschätzung: Gerschgorin-Kreise . . . 122

9.6 Iterative Verfahren . . . 124

9.7 Singulärwertzerlegung/singular value decomposition (SVD) . . . 127

9.8 Bilbiotheksroutinen . . . 128

10 Interpolation 130 10.1 Taylor-Interpolation . . . 132

10.2 Lagrange-Interpolation . . . 133

10.3 Interpolation mit rationalen Funktionen . . . 137

10.4 Stückweise/Spline-Interpolation . . . 138

10.5 Trigonometrische Interpolation . . . 141

10.6 Mehrdimensionale Interpolation . . . 145

11 Ausgleichungsrechnung 147 11.1 Hintergrund: Statistik . . . 148

11.2 Allgemeine lineare Regression . . . 151

11.3 Nichtlineare Regression . . . 156

3

0 Vorbemerkungen

Ziele der Veranstaltung:

• wie bringe ich meinen Computer dazu, nicht nur vorgefertigte Programme anderer Leute abzuspulen, sondern meine eigenen Aufgaben auf meine eigene Weise zu lösen;

• erstes Grundverständnis von numerischer i.Ggs. zu analytischer Mathematik;

• Erwerb nötigen Hintergrundwissens für den Umgang mit beliebigen vorgefertigten Program- men aus dem Bereich angewandte Mathematik und numerische Simulation in Physik, Che- mie,...

• Erkennung von Fehlerquellen bei numerischen Rechnungen

KEINE Ziele der Veranstaltung:

• detaillierte Einführung in spezielle tools;

• insbesondere: keine Einführung in Orca oder andere TheoChem-Pakete

• physikalisch- oder theoretisch-chemischer Stoff

• Komplette Abdeckung aller high-end numerischen Methoden

Anmerkung: Teile der Vorlesung basieren auf der Vorlesung von Prof. E. Pehlke, „Numerische Methoden für Physiker“.

procedere:

• 11 verschiedene thematische Kapitel; zu jedem Kapitel gibt es:

(üblicherweise abwechselnd; aber: Ansagen beachten!)

– Zwei Vorlesungsstunden: Einführung in die theoretischen Grundlagen;

Besprechung des zugehörigen Aufgabenblatts mit mehreren Aufgaben – Eine Übungsstunde + eine Vorlesungsstunde, oder zwei Übungsstunden;

Besprechung der Aufgaben

• alle Skriptseiten und alle Aufgabenblätter von meiner homepage als PDF-Dateien herunterlad- bar:http://ravel.pctc.uni-kiel.de→Teaching

→Einführung in die numerische Mathematik für Chemiker (Achtung: Inhalt kann sich während des Semesters noch ändern)

• die Benotung erfolgt über die selbstständige Lösung der Übungen.

– abzugeben sind:

∗ Skripte/Funktionen (^∗.m-Dateien)

∗ Erläuterungen/Antworten auf die Fragen (als pdf-Datei)

– Abgabe der Aufgaben am Besten per e-mail:larsson@pctc.uni-kiel.de – BeiFragen, bitte melden!

5

0.1 Inhaltsplan

Voraussetzungen: MNF-chem0102/0202 „Mathematik für Chemiker 1&2“

Inhalt: (optional)

-1- Octave/Matlab: Einführung in die grundlegenden Konzepte -2- Zahlendarstellung: Ganze Zahlen, Gleitkommazahlen (IEEE 754) -3- Fehlerquellen und -fortpflanzung, Konditionierung

-4- Differentiation

-5- Integration: Trapezregel, Simpson, Gauß, adaptiv,mehrdimensional, Monte-Carlo

-6- Suche nach Nullstellen und Extremwerten: Intervallschachtelung, Bisektion, regula falsi, Newton-Raphson, quasiNewton, mehrdimensional, conjugate gradient,DIIS

-7- Differentialgleichungen Euler, Runge-Kutta, Steifheit,adaptiv, Bulirsch-Stoer, Adams-Moulton -8- Lösung linearer Gleichungssysteme: Zerlegungen (LU, Cholesky), Kondition, Matrixinversion,

direkte Berechnung (conjugate gradient, MINRES, GMRES), Präkonditionierung

-9- Matrixdiagonalisierung: Jacobi-Rotation, Singulärwertzerlegung,direkte Berechnung (Vektori- teration, Lanczos-Verfahren), Abschätzung (Gerschgorin-Kreise)

-10- Interpolation: Polynomial, Spline, Trigonometrisch, Fouriertransformation (FFT), mehrdimen- sionale Interpolation.

-11- Regression:χ₂-Statistik, lineare und nichtlineare Regression (Levenberg-Marquardt),

Wochenplan:

Woche Thema

27.10.2015 Vorbesprechung

02.11.2015–06.11.2015 -1- -2- 09.11.2015–13.11.2015 -3- -4- 16.11.2015–20.11.2015 -5- 23.11.2015–27.11.2015 -6- 30.11.2015–04.12.2015

07.12.2015–11.12.2015

14.12.2015–18.12.2015 -7- 07.01.2016–08.01.2016

11.01.2016–15.01.2016 -8- 18.01.2016–22.01.2016 -9- 25.01.2016–29.01.2016

01.02.2016–05.02.2016 -10- 08.02.2016–12.02.2016 -11-

Literaturempfehlungen:

Allgemein

• W. H. Press, B. P. Flannery, S. A. Teukolsky und W. T. Vetterling: „Numerical Recipes 3rd Edition“, Cambridge University Press, Cambridge, 2007:

sehr gut geschriebene Einführung in die Grundalgorithmen der Numerischen Mathematik, mit Beispielprogrammen; weit verbreitetes Standardwerk.

• M. Hermann, „Numerische Mathematik“, Oldenburg, 2006

Gut zugängliches Buch. Leider fehlt die Behandlung von Differentialgleichungen und von Minimierung.

• Roland W. Freund, Ronald W. Hoppe: „Stoer/Bulirsch: Numerische Mathematik 1 und 2“, Springer, Berlin, 10. Auflage, 2007:

bewährtes Standardwerk deutscher Mathe-Fakultäten zur Numerischen Mathematik, bekannt für klare Darstellung. Etwas „mathematischer Stil“.

• Peter Deuflhard, Andreas Hohmann: „Numerische Mathematik 1&2&3: Eine algorithmisch orientierte Einführung“, de Gruyter, Berlin, 4. Auflage, 2008;

Band 1 ist eine allgemeine Einführung mit verschiedenen Themen; klarer Text, leider keine Behandlung von Minimierung. Band 2 und 3 beschäftigen sich ausgiebig mit Differentialglei- chungen.

• A. Neumaier: „Introduction to Numerical Analysis“, Cambridge University 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 zugänglich.

• G. H. Golub und C. F. van Loan: „Matrix Computations“, Johns Hopkins University Press, 2013:

Klassiker zur linearen Algebra; ziemlich „mathematischer“ Stil; für Detailinformationen zu linearen Gleichungssystemen und zum Eigenwertproblem.

• P. E. Gill, W. Murray, M. H. Wright: „Practical Optimization“, Academic Press, 1981:

Spezialbuch zur Optimierung mit vielen Tipps für die Praxis. Die Herleitungen sind leider nicht immer selbsterklärend.

Matlab/Octave-spezifisch

• Cleve B. Moler, „Numerical Computing with Matlab“, SIAM, 2013,http://de.mathworks.

com/moler/chapters.html:

vom Gründer von Matlab; sehr gut für den Einstieg und kostenlos verfügbar mit einger großen Ansammlung an Skripten.

• Alfio Quarteroni, Fausto Saleri, „Wissenschaftliches Rechnen mit MATLAB“, Springer, Berlin, 2005Englische Ausgabe (neuere Auflage): Alfio Quarteroni, Fausto Saleri, „Scientific Computing with MATLAB and Octave“, Springer, Berlin, 2006

7

Warum Numerik wichtig ist

Einige Fehler in der realen Welt, die aufgrund von numerischen Problemen resultierten:¹

• Im ersten Golfkrieg hat eine Patriot-Rakete 28 US-Soldaten getötet, da die Zeit in Inkrementen von¹/10gezählt wurde.¹/10hat allerdings keine endliche binäre Zahlendarstellung, sodass Rundungsfehler auftraten.

• 1996 stürzte eine unbemannte Ariane5-Rakete ab, da eine Gleitkommazahl in eine ganze Zahl konvertiert wurde und die Zahlendarstellung der ganzen Zahl nicht mehr ausreichte. Siehe dazu auch Abschnitt 2.3.

• 2014 war die Zahl der Klicks vom Youtube-Video “Gangnam Style” auf einmal negativ, da die Zahl zu groß für die dahinterliegende Binärdarstellung war (siehe Abschnitt 2.3).²

• 1992 wurde bei der Schleswig-Holsteiner Landtagswahl aufgrund eines Rundungsfehlers fälschlicherweise angenommen, dass die ParteiDie Grünendie 5 %-Hürde überschritten haben.

1http://ta.twi.tudelft.nl/users/vuik/wi211/disasters.html

2https://plus.google.com/wm/4/+YouTube/posts/BUXfdWqu86Q

1 Minimaleinführung in octave

• Octave ist eine freie Software/Programmiersprache zur Lösung von numerischen Problemen

• Octave ist ein Klon der proprietären Software matlab und enthält die meisten Funktionen von matlab. Matlab-Skripte sind normalerweise auch in octave ausführbar und umgekehrt.

• Octave bietet allerdings einige Syntaxerweiterungen an und verhält sich im Vergleich zu matlab in kleinen Details unterschiedlich, siehe auchhttp://wiki.octave.org/FAQ#How_

is_Octave_different_from_Matlab.3F

Auf die Syntaxerweiterungen sollte aber verzichtet werden, um Kompatibilität mit matlab zu gewährleisten.

Vorteile

• Die Programmiersprache ist sehr einfach und sehr gut für numerische Probleme geeignet.

• Octave und Matlab sind sehr gut dokumentiert

⇒ Sehr leichter Einstieg, speziell für Programmieranfänger

• Es gibt eine große Bibliothek an bereits implementierten (numerischen) Funktionen

• Octave kann mit der sehr schnellen (proprietären) Intel MKL-Bibliothek oder anderen Bibliotheken benutzt werden.¹ Damit sind viele numerische Rechnungen aufbauend auf BLAS/LAPACK/FFT-Aufrufe genauso schnell wie in entsprechenden Fortran- oder C-Codes.

• Octave kann Problemlos auf verschiedene Betriebssysteme laufen

• Matlab/Octave istdieEntwicklungsumgebung der Numeriker

• Vereinzelt gibt es sogar professionelle Quantenchemie-Programmpakete: RESCU ist ein matlab- basiertes Dichtefunktionaltheorieprogramm, mit dem Berechnungen von sehr großen Syste- men (>5000Atome) möglich sind²

Nachteile

• Die leichte Syntax verhindert eine einfache Programmierung von großen und komplexen Programmen, wie sie in moderneren Programmiersprachen wie C++ oder python möglich ist.

• Octave ist in Vergleich zu matlab etwas weniger schnell (u. a. keine Just-In-Time-Kompilierung).

Siehe aber der Hinweis auf die Benutzung der MKL-Bibliothek in Octave.

• Je nach Anwendungsfall (viele Schleifen, keine Verwendung von Bibliotheksroutinen möglich) kannoctave/matlab sehr langsam sein

1https://software.intel.com/en-us/articles/using-intel-mkl-in-gnu-octave

2http://arxiv.org/abs/1509.05746

9

1.1 Installation

• Eine Dokumentation für die Installation ist unter

https://www.gnu.org/software/octave/download.html zu finden.

• Octave sollte für alle gängigen Linux/Unix-artigen Betriebssysteme über das entsprechende Software-Repositorium zu installieren sein.

• Für windows reicht es, die entsprechende zip-Datei unter https://ftp.gnu.org/gnu/

octave/windows/herunterzuladen. In dem octave-Ordner befindet sich die ausführbare Datei octave.exe.

• Für MacOS: siehe

http://wiki.octave.org/Octave_for_MacOS_X

1.2 Literatur

Octave ist ein mächtiges Werkzeug mit vielen Funktionen. Im Rahmen dieser Vorlesung werden nur die absolut essentiellen Grundlagen gezeigt. Daher ist für eine weitergehende Arbeit mit octave (oder matlab) ein zusätzliches Literaturstudium unerlässlich.

• Für konkrete Fragen bietet sich eine Suchmaschinensuche mit dem zusätzlichen Stichwort

„matlab“ an. Die sehr ausführliche Dokumentation der matlab-Funktionen kann für die meis- ten Funktionen auch für octave benutzt werden. So führt eine Suche mit dem Stichwort

“matlab diagonalisation” direkt zur Seite dereig-Funktion, welche eine Eigenwertzerlegung („Diagonalisierung“) durchführt. Die meisten Beispiele dort sind auch in octave möglich.

• Die offizielle Dokumentationsseite von Matlab ist unterhttp://de.mathworks.com/help/

matlab/index.htmlzu finden.

• Achtung: Ab und an scheinen die Matlab-Seiten nur nach einer Produktregistrierung zugänglich zu sein. Links, die über eine Suchmaschinensuche erhalten werden, scheinen davon aber nicht betroffen zu sein.

• Weitere Quellen:

– Die Dokumentation in octave nebsthelp- unddoc-Befehl (siehe Abschnitt 1.4.8) – Eine etwas längere Einführung:

http://math.jacobs-university.de/oliver/teaching/iub/resources/octave/

octave-intro/octave-intro.html – Das sehr ausführliche, offizielle Handbuch:

https://www.gnu.org/software/octave/octave.pdf – Das offizielle Wiki:

http://wiki.octave.org/Main_Page

– Video-Tutorial (ohne GUI³, welche erst vor kurzem in Octave implementiert wurde):

https://www.youtube.com/watch?v=TWqYSOiSTEA&list=

PLlcAcGxlil6PsXCQQOw1tpJzxSzZBjln-&index=1 – Video-Tutorial für die GUI:

https://www.youtube.com/watch?v=J85b2LLoizc

3GUI=Grapische Benutzerschnittstelle

1.3 Weitere Bibliotheken

• Unter

http://octave.sourceforge.net/index.html

finden sich weitere Pakete, um die Funktionialität von octave zu erweitern. So finden sich dort u. a. spezielle Routinen für genetische Algorithmen, Quaternionen (verallgemeinerte komplexe Zahlen), Bild- und Signalverarbeitung, Statistik und vieles mehr

• Matlab bietet eine ähnliche Seite unter

https://www.mathworks.com/matlabcentral/fileexchange/

ACHTUNG: Jeder kann dort Code hochladen, sodass die Qualität der Pakete sehr stark divergiert.

• Desweiteren bietet sich eine Internetsuche an. Viele Numeriker veröffentlichen ihre Matlab- Skripte/Pakete auf ihren Webseiten.

1.4 Ein erster Einstieg

Die graphische Benutzerschnittstelle (GUI) von octave:

Die einzelnen Fenster:

1 Das derzeitige Verzeichnis 2 Verzeichnisnavigator

3 Momentan verwendete Variablen mit Attributen (Workspace) 4 Vorher benutzte Befehle

5 Ausschnitt aus dem Befehlsfenster mit Ein/Ausgabe

11

6 Befehlsfenster (gerade ausgewählt) 7 Editor (nicht ausgewählt)

8 Dokumentation (nicht ausgewählt)

1.4.1 Das Befehlsfenster

Im Befehlsfenster kann interaktiv gearbeitet werden. Mit „> >“ beginnende Zeilen geben Befehle an. Drücken der Enter-Taste führt die Befehle aus. Herzstück einer jeden Programmiersprache sind die Variablen, denen per „=“ Zahlen oder kompliziertere mathematische Objekte (z. B. Matrizen) zugeordnet werden. Letztere können Resultate mathematischer Ausdrücke oder Funktionsaufrufe sein. Das Gleicheitszeichen ist also alsZuordnungszeichenzu verstehen. Mathematische Gleicheit wird per „==“ überprüft (mehr dazu später). Es ist dringend zu empfehlen, Variablennamen nach ihrer Funktion im Programm sinnvoll zu wählen, auch wenn sich dadurch die Schreibarbeit erhöht.

Der Variablentyp (Zeichenkette/String, Gleitkommazahl, Matrix etc.) wird automatisch bestimmt und kann sich nach einer anderen Zuordnung wieder ändern.

Wenn hinter einem Befehl ein „;“ steht, wird die Befehlsausführung im Hintergrund durchgeführt und das Ergebnis nicht auf dem Bildschirm ausgegeben:

1 >> variable1 = 41 + 1

2 variable1 = 42

3 >> %Alles hinter einem Prozentzeichen ist ein Kommentar und wird ignoriert

4 >> variable2 = 42; %Ausgabe wird unterdrückt

5 >> variable2 = variable2 + 3; % " variable2 + 3" wird berechnet und wieder in variable2 gespeichert

6 >> variable2

7 variable2 = 45

8 >> variable1 + 1 %Ausdrücke ohne Zuweisung werden automatisch der Variable "ans"

zugewiesen

9 ans = 43

10 >> ans

11 ans = 43

12 >> ans = ans ∗ exp(0.4) %alle gängigen Funktionen stehen zur Verfügung

13 ans = 64.148

14 >> pi %auch einige Konstanten sind vordefiniert

15 ans = 3.1416

16 >> sin(pi)

17 ans = 1.2246e−16

18 >> cos(pi)

19 ans = −1

20 >> a = ' Buchstaben' %Auch Zuweisungen zu anderen Objekten wie Text (" strings ") sind mö glich

21 a = Buchstaben

22 >> Variable1 = −1; %auf Groß−und Kleinschreibung wird geachtet!

23 >> variable1 %dies ist eine andere Variable

24 variable1 = 42

25 >> clear variable1 %variable1 wird gelöscht

26 >> variable1

27 error: ' variable1 ' undefined near line 1 column 1

28 >> Variable1 %ist immer noch vorhanden

29 Variable1 = −1

30 >> clear all %alle Variablen werden gelöscht

Nach Eingabe der gezeigten Befehle füllt sich allmählich das “Workspace”-Fenster 3 mit den definierten Variablen und die Befehlsgeschichte 4 mit den ausgeführten Befehlen.

In Zeile 17 ist eine fundamentale Eigenschaft von Gleitkommazahlen zu sehen:sin(π) ergibt nicht genau Null, sondern in diesem Fall1.2246·10⁻¹⁶, was an Rundungsfehlern aufgrund der beschränkten Darstellung von reellen Zahlen auf dem Computer liegt. Für cos(π)kommen die Rundungsfehler zufälligerweise nicht vor. Mehr dazu im Kapitel 2.

Die Fehlermeldung in Zeile 27 ist etwas missverständlich: Es steht dort “line 1”, da es im interakti- ven Modus keine Zeilen gibt. Die Fehlermeldung ist eigentlich für Skripte gedacht (siehe Abschnitt 1.4.6).

1.4.2 Vektoren und Matrizen

Die Definition von Vektoren und Matrizen ist denkbar einfach:

1 >> v = [1,2,3] %ein Zeilenvektor

2 v =

3

4 1 2 3

5

6 >> a = [1; 2; 3;] %ein Spaltenvektor

7 a =

8

9 1

10 2

11 3

12

13 >> v = 1:3 %automatische Generierung eines Zeilenvektors : Starte bei 1 und ende bei 3

14 v =

15

16 1 2 3

17

18 >> 1:2:8 %Starte bei 1 und Ende in Zweierschritten bei 8

19 ans =

20

21 1 3 5 7

22

23 >> size(v) %gibt die Größe des Vektors v aus: Zahl der Zeilen x Zahl der Spalten

24 ans =

25

26 1 3

27

28 >> size(a)

29 ans =

30

31 3 1

32 >> A = [1 2 3; %eine 2x3−Matrix

33 4 5 6]

34 A =

13

35

36 1 2 3

37 4 5 6

38

39 >> size(A)

40 ans =

41

42 2 3

43

44 >> A (:,1) %gibt die erste Spalte aus

45 ans =

46

47 1

48 4

49

50 >> A (:,3) %gibt die dritte Spalte aus

51 ans =

52

53 3

54 6

55 >> A (2,:) %gibt die zweite Zeile aus

56 ans =

57

58 4 5 6

59

60 >> A (2,2:3) %gibt die zweite bis dritte Spalte in der zweiten Zeile aus

61 ans =

62

63 5 6

64 >> A' %transponiert die Matrix

65 ans =

66

67 1 4

68 2 5

69 3 6

70 >> v ∗ A' %berechnet das Matrix−Vektor−Produkt des Zeilenvektors v mit der transponierten Matrix A

71 ans =

72

73 14 32

74 >> A = A' ∗ A%ein Matrix−Matrix−Produkt

75 A =

76

77 17 22 27

78 22 29 36

79 27 36 45

80 >> A .^2 %elementweise Operationen können mit einem vorangestellten "." ausgeführt werden. Hier werden alle Elemente von A mit 2 potenziert

81 ans =

82

83 289 484 729

84 484 841 1296

85 729 1296 2025

86 >> B = zeros(3,2) % erstellt eine 3x2−Matrix mit Nullen

87 B =

88

89 0 0

90 0 0

91 0 0

92 >> ones(2,3)

93 ans =

94

95 1 1 1

96 1 1 1

97 >> eye(2)

98 ans =

99

100 Diagonal Matrix

101

102 1 0

103 0 1

Komplizierte „Scheibenoperationen“ (slicing) wie in Zeile 60 lassen sich am Besten durch Ausprobie- ren bei einer größeren Matrix nachvollziehen. Intern werden Vektoren als Matrizen mit nur einer Spalte- oder Zeile behandelt, wie es in Zeile 23 zu sehen ist.

In octave gibt es neben dersize-Funktion eine Fülle an weiteren Funktionen. Dazu sei auf die Dokumentation und die Literatur (siehe Abschnitte 1.2 und 1.4.8) verwiesen. Hier sei nur als weiteres Beispiel die Eigenwertzerlegung einer Matrix, bestehend aus Zufallszahlen, gezeigt:

1 >> Arand =rand(3) %generiert eine Matrix aus Zufallszahlen

2 Arand =

3

4 0.406312 0.473776 0.744493

5 0.088445 0.420816 0.399403

6 0.131595 0.297414 0.832140

7

8 >> A = 0.5 ∗ (Arand' + Arand) %symmetrisiert die Matrix

9 A =

10

11 0.40631 0.28111 0.43804

12 0.28111 0.42082 0.34841

13 0.43804 0.34841 0.83214

14

15 >> eigVal = eig(A) %gibt die Eigenwerte der nun symmetrischen Matrix A aus

16 eigVal =

17

18 0.10745

19 0.22655

20 1.32527

21

22 >> [eigVecs , eigVal ] = eig(A) %gibt die Eigenvektoren−und−werte aus; letztere als

15

Diagonalmatrix

23 eigVecs =

24

25 0.85370 0.16906 0.49256

26 −0.41269 0.79650 0.44190

27 −0.31761 −0.58053 0.74974

28

29 eigVal =

30

31 Diagonal Matrix

32

33 0.10745 0 0

34 0 0.22655 0

35 0 0 1.32527

36 >> eigVecs ∗ eigVal ∗ eigVecs ' %generiert wieder die Matrix A

37 ans =

38

39 0.40631 0.28111 0.43804

40 0.28111 0.42082 0.34841

41 0.43804 0.34841 0.83214

Zeilen 15 und 22 zeigen, dass sich dieAusgabeeiner Funktion je nach Zuweisung unterscheiden kann:

Wirdeignur einer Variable zugewiesen, so werden nur die Eigenwerte ausgegeben.[U,D] = eig(A) gibt dahingegen die EigenvektormatrixUnebst der Diagonalmatrix mit den EigenwertenDaus. Die Notation[U,D]ist alsVektor von Variablenzu verstehen. In den Variablen können alle möglichen Ausdrücke (Zahlen, Vektoren, Matrizen, Strings,...) gespeichert werden.

1.4.3 Kontrollstrukturen und Schleifen

Neben den Variablen sind Kontrollstrukturen und Schleifen ein weiteres Herzstück einer jeden Programmiersprache. Mittels Kontrollstrukturen ist das Verhalten eines Programmes je nach Eingabe manipulierbar:

1 >> a = −4;

2 >> if a > 0

3 disp('a ist echt positiv . ') ; %Zeige eine Nachricht auf dem Display an

4 elseif a < 0

5 disp('a ist echt negativ . ') ;

6 elseif a == 0 %solch eine Abfrage ist fehlerbehaftet ! Siehe den Text unten.

7 disp('a ist genau Null. ') ;

8 else

9 error(' a ist unmöglich!') ; %gebe einen Fehler aus

10 end %ende der if−Abfrage

11 a ist negativ .

12 >> v = [0, 9, 8];

13 >> if any(v > 8) %any ist eine weitere nü tzliche Funktion

14 disp('v beinhaltet mind. ein Element größer 8') ;

15 end

16 v beinhaltet mind. ein Element größer 8

Zeile 6 ist fehlerbehaftet, da aufgrund der finiten Genauigkeit von Gleitkommazahlen exakte Nullen selten vorkommen. In Abschnitt 1.4.1 haben wir schon gesehen, dass auf dem Computersin(π)nicht unbedingt Null ist. Besser wäre z.B. die Abfrage abs(a) < 1e-14, welche testet, ob der Absolutwert vonakleiner als10⁻¹⁴ist.

Boolsche Ausdrücke sind durch die Benutzung von&&für UND, sowie||für ODER möglich. Ein Ausrufezeichen!stellt eine Verneinung dar. Da das Gleichheitszeichen schon für die Variablenzu- weisung belegt ist, wird==für die Überprüfung auf Gleichheit benutzt. In octave steht1für WAHR und0für FALSCH:

1 >> a = 9;

2 >> b = −7;

3 >> if a > 0 && b > 0

4 disp('a und b sind beide positiv ') ;

5 elseif a ~= 9

6 disp('a ist nicht 9. ') ;

7 elseif a > 0 || b < 0

8 disp('Entweder ist a positiv oder b ist negativ oder beides . ') ;

9 end

10 Entweder ist a positiv oder b ist negativ oder beides .

Schleifen, welche es erlauben, Code-Zeilen mehrmals zu durchlaufen, sind mit demfor-Konstrukt möglich:

1 >> for i = 4:5

2 disp( i ) ;

3 end

4 4

5 5

6 >> for i = 9:−1:6 %Gehe von 9 bis 6 und dekrementiere jeweils um 1

7 disp( i ) ;

8 end

9 9

10 8

11 7

12 6

13 >> v = [4,9,2,3];

14 >> for i = v %Da 4:5 ein Vektor ist ( siehe den Abschnitt zu Vektoren) , kann genauso gut ein explizit gespeicherter Vektor benutzt werden

15 disp( i ) ;

16 end

17 4

18 9

19 2

20 3

21 >> for i = v

22 disp( i ) ;

23 if i == 9

24 break; %beende die gesamte Schleife

25 end

26 end

27 4

28 9

17

29 >> for i = v

30 if i == 2

31 continue; %beende diesen Schleifendurchlauf

32 end

33 disp( i ) ;

34 end

35 4

36 9

37 3

Im letzten Beispiel wird die Ausgabe füri= 2aufgrund dercontinue-Anweisung übersprungen.

Eine weitere Möglichkeit bietet die while-Schleife an, die den Code innerhalb der Schleife so lange ausführt, bis eine Bedingung nicht mehr wahr ist:

1 >> n = 10;

2 >> f = n; %berechne die Fakultät von n = 10

3 >> while n > 1

4 n = n−1;

5 f = f∗n;

6 end

7 >> f

8 f = 3628800

9 >> factorial (10)

10 ans = 3628800

Wenn auf Schnelligkeit geachtet werden muss, sollten Schleifen in octave wenn möglich vermieden und durch geschickte Anwendung der schon in octave implementierten Funktionen ersetzt werden.⁴ Folgendes Beispiel invertiert die Reihenfolge der Elemente eines Vektors, was allerdings schon in der Funktionfliplrimplementiert ist:

1 >> v = 1:4

2 v =

3

4 1 2 3 4

5

6 >> v2 = v

7 >> for i = 1:length(v)

8 v( i ) = v2(length(v) −i + 1) ;

9 end

10 >> v

11 v =

12

13 4 3 2 1

14 >> fliplr (1:4)

15 ans =

16

17 4 3 2 1

4Die Bibliotheksroutinen von ovtave sind optimiert und in der deutlich schnelleren Sprache C++ geschrieben bzw.

greifen selbst auf optimierte Bibliotheken zu.

1.4.4 Input/Output

Mitformatändert man die Formatierung der Standardausgabe:

1 >> 100 ∗ pi

2 ans = 314.16

3 >> format long

4 >> 100 ∗ pi

5 ans = 314.159265358979

6 >> format long e

7 >> 100 ∗ pi

8 ans = 3.14159265358979e+02

9 >> format %Rückstellung auf Standardwerte

10 >> 100 ∗ pi

11 ans = 314.16

Mehr dazu findet sich unterhelp format(siehe Abschnitt 1.4.8). Mitfprintfsind auch kompli- ziertere Ausgaben möglich:

1 >> fprintf ('%6d ist eine große Zahl\n',5e6) ; % '\ n' stellt einen Zeilenumbruch dar

2 5000000 ist eine große Zahl

3 >> fprintf ('%3.2f ist pi auf zwei Nachkommastellen genau\n',pi)

4 3.14 ist pi auf zwei Nachkommastellen genau

5 >> fprintf ('%6.5f ist pi auf fünf Nachkommastellen genau\n',pi)

6 3.14159 ist pi auf fünf Nachkommastellen genau

7 >> fprintf ('%f\n' ,[2.0,4.0,2.124]) ; % '%f ' nimmt ein Standardformat

8 2.000000

9 4.000000

10 2.124000

Für eine detailliere Dokumentation, siehe

https://de.mathworks.com/help/matlab/ref/fprintf.html

Dateien lassen sich überfopen,fprintfundfcloseverwalten und beschreiben:

1 >> v = [2.0, 1.999, pi, cos(0.1) ]

2 v =

3

4 2.00000 1.99900 3.14159 0.99500

5

6 >> dateiID = fopen(' werte.dat',' w') ; % 'w' heißt, dass die Datei beschrieben wird (' write ')

7 >> fprintf (dateiID ,' #einige Werte\n') ; % schreibe eine Kommentarzeile, üblicherweise mit vorangestellter Raute

8 >> fprintf (dateiID ,' %f\n',v) ; %gebe die Werte aus

9 >> fclose(dateiID) ; %schließe die Datei

Anschließend sollte im derzeitigen Verzeichnis (siehe 1 ) die Datei „werte.dat“ erscheinen. Das Einlesen erfolgt analog mitfscanfund'r'anstelle von'w'infopen. Große Matrizen oder ähnliches lassen sich aber am Besten mittelssaveundloadein- und auslesen. Siehe dazu die Dokumentation.

1.4.5 Graphiken erstellen

Graphiken können mit demplot-Befehl erstellt werden:

19

1 >> graphics_toolkit ('gnuplot') ; %evt. nötig , siehe Text

2 >> x = linspace(−10,10,100) ; %gibt ein ä quidistantes Gitter von −10 bis 10 mit 100 Punkten aus

3 >> y = sin(x) .∗ exp(−abs(x)) ; %elementweise Multiplikation ist '.∗' und nicht '∗'!

4 >> y2 = cos(x) .∗ exp(−abs(x)) ;

5 >> plot(x,y,x,y2,' .. g') ; %Zeichnet y(x) und y2(x) ; letzteres grün und gepunktet

6 >> grid; %plotted ein Gitter

7 >> legend(' y1',' y2') ;

8 >> xlabel(' x') ;

9 >> ylabel(' y') ;

Der Plot kann über die Kommandozeile im Plotfenster oben Links gespeichert werden:

-0.4 -0.2 0 0.2 0.4 0.6 0.8 1

-10 -5 0 5 10

y

x

y1 y2

Mit der neuen Version ovtave 4.0 hat sich die Graphikbibliothek von gnuplot auf Qt geändert. Dies bereitet unter Linux anscheinend noch etwas Probleme, daher der Befehl in Zeile 1. Unter Windows ist dieser Befehl nicht nötig und bereitet evt. Probleme. Weiterhin zu beachten ist die elementweise Multiplikation in Zeile 3. Eine normale Multiplikation würde eine Matrix-Matrix-Multiplikation zwischen den Zeilenvektoren sin(x) und exp(-abs(x)) durchführen, was allerdings zu einem Fehler führt, da zwei Zeilenvektoren nicht miteinander multipliziert werden können.

3D-Graphiken sind analog möglich, z. B. mit den Befehlensurf(3D-Oberflächen),mesh(3D-Gitter)

odercontour(Konturplots). Dafür sei auf die Dokumentation verwiesen.

1.4.6 Skripte schreiben

Eine interaktive Eingabe von vielen Befehlen ist auf Dauer recht mühsam. Vor allem kann man so keine Programme speichern. Die GUI stellt auch einen Editor für das Schreiben von Skripten zur Verfügung. Klicken auf 7 öffnet das Editor-Fenster:

Ein Klick auf den rot markierten Knopf oder das Drücken der TasteF5führt das Skript aus und man sieht das Ergebnis im Kommando-Fenster, 6 . Ist die Datei noch nicht gespeichert, so wird zuerst ein Speicherdialog aufgerufen. Es empfiehlt sich, die Datei mit der Endung „.m“ zu speichern, was die typische Endung von matlab-Skripten ist. Das Skript kann auch ausgeführt werden, indem der Name des Skripts (ohne Endung „.m“) in der Kommandozeile eingegeben wird.

Es ist zu empfehlen, das Editor-Fenster so zu verschieben, dass man sowohl das Editor-Fenster als auch das Kommandozeilen-Fenster sieht. Dies kann man durch Drücken der linken Maustaste auf die „Editor-Leiste“ (grün markiert im obigen Bild) und anschließendes Verschieben des Fensters (bei gedrückter Taste) erreichen.

Achtung: Das Aufrufen von Skripten fügt Variablen zum Workspace ( 3 ) hinzu. Innerhalb von Skripten kann man auch auf Variablen im Workspace zugreifen, was meistens nicht gewünscht ist.

Einclear all;in der ersten Zeile löscht alle vorhandenen Variablen im Workspace.⁵

1.4.7 Funktionen schreiben

Neben Skripten sind Funktionen sehr hilfreich für eine bessere Gliederung eines Programms. In octave mussjedeFunktion in eine separate Datei gespeichert werden. Der Name der Funktion muss mit dem Dateinamen übereinstimmen. So könnte die Datei „meine_funktion.m“ wie folgt aussehen:

1 function z = meine_funktion(x,y)

2 z = x + y;

3 end

In Zeile 1 ist angegeben, dass die Funktion die Varriablen x und y annimmt und ein Ergebnis ausliefert, was intern in der Funktion als Variablezabgespeichert wird. Die Funktion kann nun in der Kommandozeile oder auch in Skripten benutzt werden:

1 >> ergebnis = meine_funktion (2,3)

2 ergebnis = 5

5Wenn die Skripte auf Schnelligkeit ausgerichtet sein sollen, sollte davon abgesehen werden.

21

Auch komplizierte Funktionen sind möglich:

1 function [r1 , r2] = funk2(x,y)

2 a = x + y;

3 r1 = sqrt(a) ;

4 r2 = [r1 ∗2, r1 +2];

5 end

Der Aufruf erfolgt über[a,b] = funk2(ein1,ein2). Der Ausgabevektor ist wieder als Vektor von Variablen zu verstehen.

Es können alle Funktionen und Skripte aufgerufen werden, die im derzeitigen Suchpfad vorhanden sind. Der Suchpfad kann mitpathausgegeben werden und besteht aus einigen Standardverzeichnis- sen und dem aktuellen Verzeichnis ( 1 ). Überaddpathkönnen weitere Verzeichnisse hinzugefügt werden. Es bietet sich z.B. an, eigene numerische Funktionen in einem Verzeichnis zu speichern.

Skripte, welche die Funktionen aufrufen, können dann woanders gespeichert und der Pfad mittels addpathim Skript entsprechend verändert werden.

1.4.8 Die Hilfefunktion

Die mit Abstand wichtigste Funktion in octave isthelp. Für jede eingebaute Funktion in octave gibt helpeinen kurzen Eintrag über die Funktionsweise aus:

1 >> help(' eig ') %Eingabe von "help eig" ist auch möglich

2 ' eig ' is a built−in function from the file libinterp / corefcn /eig. cc

3

4 −−Built−in Function: LAMBDA =eig(A)

5 −−Built−in Function: LAMBDA =eig(A, B)

6 −−Built−in Function: [V, LAMBDA] =eig(A)

7 −−Built−in Function: [V, LAMBDA] =eig(A, B)

8 Compute the eigenvalues (and optionally the eigenvectors ) of a

9 matrix or a pair of matrices

10

11 The algorithm used depends on whether there are one or twoinput

12 matrices , if they are real or complex, and if they are symmetric

13 (Hermitian if complex) or non−symmetric.

14

15 The eigenvalues returned by ' eig ' are not ordered.

16

17 See also : eigs , svd.

18

19 Additional help for built−in functions and operators is

20 available in the online version of the manual. Use the command

21 ' doc <topic>' to search the manual index.

22

23 Help and information about Octave is also available on the WWW

24 at http :// www.octave.org and via the help@octave.org

25 mailing list .

Das Dokumentationsfenster, 8 , ist mitunter auch sehr hilfreich. Eingabe vondoc('eig')führt z.B. direkt zum Eintrag voneigin der Dokumentation. Angegeben ist derselbe Text, wie er auch mittelshelpzu sehen ist, allerdings sind im Fenster noch weitere Funktionen und ein allgemeiner Text zu sehen.

1.4.9 Debuggen

Siehe die Live-Vorführung.

1.5 Programmierdisziplin!

• Warum kein „Freistil-Programmieren“?

– Programme anderen verständlich machen – eigene Programme später wieder verstehen – Schreiben längerer Programme erleichtern

• „Strukturiertes Programmieren“ = Programme übersichtlich und verständlich gestalten; insbe- sondere:

– sinnvolle Variablennamen

– gleiche Variablen in verschiedenen Programmblöcken gleich nennen

– längere, abgeschlossene Unteraufgaben in Funktionen verlegen, auch wenn nur einmal verwendet

– for- und if-Blöcke usw. systematisch und immer gleich einrücken

• Kommentare!Natürlich keine trivialen, aber im Zweifelsfall immer viel mehr als man zunächst denkt.

• Tips zurProgrammentwicklung:

– modular! = Niemals ein Programm einfach von Anfang bis Ende herunterschreiben, sondern jeden sinnvollen Abschnitt sofort testen. Nachträgliche Fehlersuche immer viel schwerer.

– erleichternde Hilfsmittel wie Debugging ausnutzen

• Verzeichnisse übersichtlich halten:

– Extra-file mit Verzeichnis-Inhaltsangabe

– nicht mehr benötigte Programmvarianten sofort löschen

– sinnvolle file-Namen: je länger, desto besser; Varianten eines Programms und Zusam- mengehörigkeit Programm/Input/Output sollte erkennbar sein.

23

2 Zahlendarstellung

2.1 Natürliche Zahlen

• Eine Zahln ∈N⁰ wird als Summe n=

m−1

X

i=0

d_ibⁱ (2.1)

mit einerBasisb∈N^≥2undZiffernd₀, d₁, . . . , dm−1 ∈[0, b−1]dargestellt.

• Die Ziffern lassen sich ausndurch wiederholte Division durch die Basisbgewinnen. Wirdb nicht angezeigt, so ist die Zehnerbasis gemeint.

• Notation:nreffiZBasis.

• Die Ziffern lesen sichvon rechts nach links.

• Beispiel 23:

– 10₂₃ = 0·23⁰ + 1·23¹ = 23

– 23₁₀ = 3·10⁰ + 2·10¹ = 3 + 20 = 23 – 27₈ = 7·8⁰+ 2·8¹ = 7 + 16 = 23

– 113₄ = 3·4⁰+ 1·4¹+ 1·4² = 3 + 4 + 16 = 23

– 10111₂ = 1·2⁰+ 1·2¹+ 1·2²+ 0·2³+ 1·2⁴ = 1 + 2 + 4 + 0 + 16 = 23

• Auf dem Computer wird die Zweier-/Binärdarstellung verwendet.¹

• In diesem Kontext werden die ZiffernBits genannt.

• 8 Bits werden zu einemBytegebündelt.

• Die maximal darstellbare Binärzahl beträgt11. . .1

| {z }

mmal

2 = 2^m−1(−1wegen der Darstellung der Null)

2.2 Ganze Zahlen

Im weiteren Verlauf: Beschränkung auf die Binärdarstellung.

• Einfachste Darstellung: Produkt einer natürlichen Zahl mit einem Vorzeichen-Bitσ ∈ {−1,1}.

• Nachteil: Die Null kommt doppelt vor (±0).²

• Besser:Zweierkomplementdarstellung(als2T abgekürzt)

n=

m−2

X

i=0

d_i·2ⁱ−d_m−1·2^m−1 (2.2)

1Die bei der Systemprogrammierung häufig verwendete Hexadezimaldarstellung wird hier nicht behandelt.

2Für Gleitkommazahlen ist dies sogar erwünscht, nicht aber für ganze Zahlen.

• Beispiel010111_2T = 1·2⁰+ 1·2¹+ 1·2²+ 0·2³+ 1·2⁴−0·2⁵ = 23

• Beispiel101001_2T = 1·2⁰+ 0·2¹+ 0·2²+ 1·2³+ 0·2⁴−1·2⁵ = 1 + 8−32 = −23

• 0000101112T = 23

• 111101001_2T = 1·2⁰ + 0 ·2¹ + 0·2² + 1·2³ + 0 ·2⁴ + 1 ·2⁵ + 1·2⁶ + 1·2⁷−1·2⁸ = 1 + 8 + 32 + 64 + 128−256 =−23

y Die Darstellung negativer Zahlen ändert sich je nach Länge der Zifferndarsellung: Es werden Einsen hinzugefügt. Für positive Zahlen werden Nullen hinzugefügt.

• Zahlenbereich fürmBits:[10. . .0_2T,0111. . .1_2T] = [−2^m−1,2^m−1−1]

• −1ist immer111. . .1_2T

• Vorzeichenbehaftete Darstellungen werden auch als signed bezeichnet. Das Gegenteil ist entsprechendunsigned.

2.3 Zahlenbereich und Über-/Unterschreitung

• Zahlenbereich von gebräuchlichen Zahlendarstellungen:

# bits Vorzeichen octave-Name min max

8 unsigned uint8 0 255

16 unsigned uint16 0 65 535

32 unsigned uint32 0 4 294 967 295

64 unsigned uint64 0 18 446 744 073 709 551 615

8 signed int8 −128 127

16 signed int16 −32 768 32 767

32 signed int32 −2 147 483 648 2 147 483 647

64 signed int64 −9 223 372 036 854 775 808 9 223 372 036 854 775 807

• Octave-Routine:intmax,intmin

• Eine Zahlenbereichsüberschreitung wird in octavenicht als Warnung ausgegeben. Anstelle dessen wird nur der größte/kleinste Wert angenommen:

1 >> a = int8 (127) %gebe eine ganze Zahl mit einer Darstellung von 8 Bits zurück

2 a = 127

3 >> a = a + 1 %addiere

4 a = 127

5 >> a = int8(−128);

6 >> a −1

7 ans = −128

8 >> a + 20000

9 ans = 127

• In anderen Programmiersprachen (C, Fortran,...) ist das Verhalten fataler: Additions/Subtrak- tionsoperationen werden so ausgeführt, als wenn noch weitere Bits vorhanden wären.

• Das führt dazu, dass unsigned integers wieder bei Null beginnen und signed integers bei negativen Zahlen beginnen

25

• Beispiel: Das folgende C-Programm

1 #include<stdio .h>

2 int main(){

3 signed char a; //8 bit lang , entspricht int8 in octave

4 unsigned char b; // entspricht uint8 in octave

5 a = 127;

6 b = 255;

7 printf ("a = %d\n",a) ;

8 printf ("b = %d\n",b) ;

9 a = a + 1;

10 b = b + 1;

11 printf ("a + 1 = %d\n",a) ;

12 printf ("b + 1 = %d\n",b) ;

13 a = a + 999;

14 b = b + 999;

15 printf ("a + 1000 = %d\n",a) ;

16 printf ("b + 1000 = %d\n",b) ;

17 }

führt nach der Übersetzung in Maschinensprache zur Ausgabe von a = 127

b = 255 a + 1 = -128 b + 1 = 0 a + 1000 = 103 b + 1000 = 231

• Dies wird alsoverflowbezeichnet und kann zu eklatanten Fehlern führen (siehe „Warum Numerik wichtig ist“)

y Ein overflow ist nicht immer leicht erkennbar.

• Compileroptionen führen zu einer Überprüfung auf overflows, was das Programm allerdings langsamer macht.

y Man sollte im vornherein abschätzen, ob der Zahlenbereich ausreicht oder nicht. Üblicherweise werden 32- oder 64-Bit-Darstellungen verwendet.

• Octave benutzt standardmäßig die 64-Bit-Darstellung

2.4 Gleitkommazahlen

• Da ein Computer nicht unendlich lange Zahlen ausRabspeichern kann, wird immer mit einer UntermengeFⁱ ⊂Rgerechnet.Fⁱist eine der Mengen der Gleitkommazahlen Welche Zahlen inFⁱenthalten sind, hängt von der jeweiligen Darstellung auf dem Computer ab.

• Eine weithin gebräuchliche Darstellung von Gleitkommazahlen richtet sich nach dem AN- SI/IEEE 754-Standard, der im folgenden behandelt wird.

• Darstellung einer Zahlxin Basisb(im folgenden Binärdarstellung):

x= (−1)^v·M ·b^e (2.3)

• v ist dasVorzeichenbit

• M is dieMantissemit Zifferndi,i∈[0, m−1]: M = d₀

b¹ + d₁

b² +· · ·+ dm−1

b^m (2.4)

• eis derExponent: Eine ganze Zahle∈[e_min, e_max].³

2.4.1 Regeln

• Fürx6= 0wirdeso gewählt, dass died₀ in der Mantisse ungleich Null ist. Man spricht von normalisierterGleitkommadarstellung.

• In Binärdarstellung wirdd₀nicht gespeichert, da es in normalisierter Darstellung immer1ist:

implizites erste Bit

• Ist die darzustellende Zahl zu klein für die normalisierte Form (e < e_min wird benötigt), so kommt es zumExponentenunterlauf (underflow)für Zahlen imsubnormalen Zahlen- bereich

• Die Darstellungsgenauigkeit sinkt im subnormalen Zahlenbereich stetig, bis alle Ziffernd_i = 0 sind.

• spezielle Zeichen/Formen:

Exponent Mantisse

normalisierte Form e∈[e_min+ 1, e_max−1] M ∈[b⁻¹,1−b^−m] denormalisierte Form e=e_min+ 1 M ∈[b^−m, b⁻¹−b^−m]

±0 e=e_min 0

±inf e=e_max 0

NaN e=e_max+ 1 6= 0

• NaN(not a number) ergibt sich z. B., wenn man durch Null teilt

• ±inf tritt bei einemExponentenüberlauf (overflow)auf. Die Zahl ist zu groß, um darge- stellt zu werden

2.4.2 Gebräuchliche Darstellungen

Einfache Genauigkeit (single precision)

• Mantissenlänge: 23 bit (+ implizites erste Bit)

• e∈[−126,129], 8 bit

• Mit Vorzeichenbit 32 bit, 4 byte

3Meist wird einBiasaufeaddiert damiteimmer positiv ist und eine Zweierkomplementdarstellung vermieden werden kann.

27

• Wertebereich der normalisierten Zahlen: 1.17 · 10⁻³⁸ − 3.4 · 10³⁸ (realmin- oder realmax(’single’))

• Wertebereich der denormalisierten Zahlen:1.4·10⁻⁴⁵−1.17·10⁻³⁸

• Genauigkeit:= 2⁻²⁴ = 6·10⁻⁸

• Darstellung im Speicher:

Doppelte Genauigkeit (double precision)

• Standardformat in octave und in den meisten numerischen Programmen

• Mantissenlänge: 52 bit (+ implizites erste Bit)

• e∈[−1022,1025], 11 bit

• Mit Vorzeichenbit 64 bit, 8 byte

• Wertebereich der normalisierten Zahlen:2.22·10⁻³⁰⁸−1.797·10³⁰⁸

• Wertebereich der denormalisierten Zahlen:3.9·10⁻³²⁴−2.22·10⁻³⁰⁸

• Genauigkeit:= 2⁻⁵³ = 1.1·10⁻¹⁶

• Auch diese Genauigkeit reicht für manche Anwendungen nochnichtaus!

2.4.3 Anmerkungen

• Die meisten reellen Zahlen lassen sich nur mit einer bestimmten Genauigkeit darstellen:

1 >> printf(' %20.17f\n' ,0.3) %0.3 kann nicht in endlicher Binä rdarstellung exakt wiedergegeben werden

2 ans = 0.29999999999999999

3 >> printf(' %20.17f\n' ,0.5) %genaue Darstellung, da 0.5 = 1/2

4 ans = 0.50000000000000000

• Dies führt zu entsprechendenRundungsfehlern, die nicht leicht zu entdecken sind:

1 >> if 1.0 −0.8 −0.2 == 0

2 disp(' 1.0 −0.8 −0.2 ist genau Null') ;

3 else

4 disp(' 1.0 −0.8 −0.2 ist NICHT genau Null');

5 end

6 1.0 −0.8 −0.2 ist NICHT genau Null

7 >> printf(' %20.17f\n',1.0 −0.8 − 0.2) ;

8 −0.00000000000000006

• Aufgrund der halbexponentiellen Darstellung (Gl.(2.3)) sind die exakt darstellbaren Zahlen logarithmisch auf der Zahlenskala verteilt.

• Beispiel fürM bestehend aus drei Bits unde∈[−1,2]:

1 4

3 8

1 2

3 4

3

1 2 2 3

5 16

7 16

5 8

7 8

7 4

5 2

7

0 2

• Spezielle Softwarepakete wie Sympy, Mathematica oder Maple sind auf das Rechnen mit beliebig hoher Genauigkeit spezialisiert.⁴

2.4.4 Maschinengenauigkeit

• Maximaler relativer Fehler beim Runden einer Zahlxaufgrund von finiter Genauigkeit:

|x−rd(x)|

|x| ≤ 0.5−2^−m

2⁻¹ = 2^−m ≡ (2.5)

rd(x)ist eine Funktion, welchexrundet, sodass sie auf dem Computer dargestellt werden kann.

• Als „Maschinenepsilon“ wird der maximale relative Abstand zweier benachbarter normalisierter Gleitkommazahlen bezeichnet (Unterscheidung in der letzten darstellbaren Stelle um 1):

2^−m

2⁻¹ = 2^1−m = 2 (2.6)

• Skriptberechne_eps.mzum Berechnen von:

1 disp(' Berechne das Maschinenepsilon') ;

2 meps = 1.0;

3 while 1.0 + 0.5 ∗ meps > 1.0

4 meps = meps ∗ 0.5;

5 end

6 fprintf (' Berechnetes Maschinenepsilon: %20.17e\n',meps);

7 fprintf (' Analytisch berechnetes Maschinenepsilon: %20.17e\n',2 ∗ 2^(−53)); % Mantissenlänge: 52 + implizites erste Bit = 53

8 fprintf (' Intern verwendetes Maschinenepsilon: %20.17e\n',eps) ; %eps ist eine Konstante in octave

• Ausgabe:

1 >> berechne_eps

2 Berechne das Maschinenepsilon

3 Berechnetes Maschinenepsilon: 2.22044604925031308e−16

4 Analytisch berechnetes Maschinenepsilon: 2.22044604925031308e−16

5 Intern verwendetes Maschinenepsilon: 2.22044604925031308e−16

• Die Ausgabe zeigt, dass das Rechnen auf dem Computer zwar mit Ungenauigkeiten versehen, der Fehler aber genau vorher bestimmbar ist. Ein Computer rechnet deterministisch.

4Auch in octave/matlab gibt es Pakete für das Rechnen mit beliebiger Genauigkeit.

29

3 Fehler und Kondition

3.1 Fehlerarten

Es gibt verschiedene Arten von Fehlern:

3.1.1 Modellfehler

• Fehler aufgrund des Unterschiedes zwischen dem realen System und dem mathematischen Modell

• Beispiele:

– Pendel mit Luftreibung, Corioliskraft etc. versus Newtonsche Gleichung im Inertialsystem – Molekül in Lösung versus Molekulardynamiksimulation in der Gasphase mit klassi-

scher Beschreibung der Atomkerndynamik und approximativer quantenmechanischer Beschreibung der Elektronen

3.1.2 Datenfehler

• Abweichung der realen Parameterwerte von den eingegebenen Daten

• Beispiele:

– Signifikante Stellen von Naturkonstanten – Messungenauigkeiten im Experiment

3.1.3 Verfahrensfehler

• Intrinsischer Fehler des numerischen Verfahrens

• Beispiele:

– Endliche Taylorreihenentwicklung – Diskretisierungsfehler

• In diesem Zusammenhang heißt ein Verfahrenkonvergent, wenn der Verfahrensfehler gegen null gehen kann

• Beispiel: Verfeinerung der Diskretisierung, z. B. bei einer Integration

• Gegenbeispiel: Divergente Reihenentwicklung, z. B. bei quantenmechanischer Vielteilchenstö- rungstheorie

3.1.4 Rundungsfehler

• Fehler aufgrund der Darstellung von Gleitkommazahlen auf dem Computer

• Beispiele

– Das Distributivgesetz gilt nicht:a·(b+c)6=a·b+a·c – Das Assoziativitätsgesetz gilt nicht:a+ (b+c)6= (a+b) +c

3.2 Einschub: Norm

• DieNormk · kin einem Vektorraum¹Xerfüllt die Bedingungen 1. kxk ≥0, x∈X

2. kαxk=|α|kxk, α∈C 3. kx+yk ≤ kxk+kyk

• BeispielX=Rⁿ(gewöhnliche Vektoren):

1-Norm kxk₁ =Pn i |x_i|

Euklidische Norm kxk₂ =pP

i|x_i|² p-Norm kxk_p =pP^p

i|x_i|^p

Maximums-Norm kxk∞= maxi|xi|

• Alle Normen aufRⁿlassen sich gegeneinander abschätzen:

C₁kxk_α ≤ kxk_β ≤C₂kxk_α, (3.1)

wobeiC₁ undC₂Konstanten sind.

3.2.1 Matrixnorm/Abbildungsnorm

• Die Abbildungsnorm einer linearen AbbildungAbzw. MatrixAist definiert als kAk ≡max

x6=0

kAxk kxk

, (3.2)

hängt also von der Definition der Norm vonxab.

• Beispiele für einen×n-Matrix:

1-Norm kAk₁ = max_j∈[1,n]Pn

i=1|A_ij|(Spaltensummennorm) Euklidische (Spektral-)Norm kAk2 =√

λ_max;λ_maxist der größte Eigenwert vonA^†A(größ- ter Singulärwert vonA, siehe SVD)

inf-Norm kAk_∞ = max_i∈[1,n]Pn

j=1|A_ij|(Zeilensummennorm)

• Der Aufwand für die numerische Berechnung von Matrixnormen variiert demnach sehr stark.

• Normdefinitionen, indem einem×n-Matrix als(m·n)-langer Vektor angesehen wird, sind auch möglich. Beispiel:Frobeniusnorm

kAk_F = v u u t

m

X

i n

X

j

|A_ij|², (3.3)

welche der euklidischen Norm eines Vektors entspricht.

1Siehe Mathematik für Chemiker II, Kapitel 2.1 und 2.1.2

31

3.3 Konditionierung mathematischer Probleme

• Aufgrund der finiten Darstellung von Zahlen auf dem Computer, sollte man bei der Fehlerana- lyse statt einer exakten Eingabexin eine Funktionf(x)von einer EingabemengeE ausgehen, welche zu einer ResultatsmengeR=f(E)führt:

E x

f(x) R

• Eine Maschinenzahlxwird demnach durchn

˜

x∈R: ^|˜x−x|_|x| ≤o

dargestellt (aus Gl. 2.5 oder, bei Datenfehler, größer, s. o.).

• DieKonditioncharakterisiert das Verhältnis zwischenE undR.

• Beispiel 1: Gutkonditionierte Schnittpunktberechnung zweier Geraden:

g

˜ g

h

˜h

r r˜

Der Schnittpunkt˜rvariiert hier ähnlich stark wie die Geradeng undg˜.

• Beispiel 1: Schlechtkonditionierte Schnittpunktberechnung (schleifender Schnitt):

1

r

Starke Schnittpunktsvariation

3.3.1 Definitionen

• Seiy=F(x)unter der Annahme vonexaktenZahlen (ohne Rundungsfehler)

• Seiy˜ =F(˜x), wobeix˜=x+δxeine geänderte Eingabe ist

• Derabsolute Fehlerist entsprechend

∆_abs =k˜y−yk. (3.4)

• Derrelative Fehlerist

∆_rel = k˜y−yk

kyk . (3.5)

• Dieabsolute Konditionszahlist

κ_abs= inf{α| k˜y−yk ≤αk˜x−xk}, (3.6)

also letztlich daskleinstmöglicheVerhältnisk˜y−yk ·(k˜x−xk)⁻¹.

• SofernF(x)differenzierbar ist, gilt

κ_abs=k∇F(x)k ≡ kF⁰(x)k. (3.7)

(hier wird die Matrixnorm gebraucht)

• Dierelative Konditionszahlist analog κ_rel= inf

α

k˜y−yk

kyk ≤αk˜x−xk kxk

, (3.8)

das kleinstmögliche Verhältnis^k˜y−yk/^kyk·(^k˜x−xk/^kxk)⁻¹.

• SofernF(x)differenzierbar ist, gilt κ_rel= kF⁰(x)k

kF(x)k · kxk=κ_abs· kxk

kF(x)k. (3.9)

• Im Folgenden: Beschränkung aufκ≡κ_rel

• Istκ-1, so ist das Problemgut konditioniert. – Fürκ <1gibt es sogar eine Fehlerverringerung – Fürκ= 1treten nur Rundungsfehler auf

• Fürκ =O(100)ist das Ergebnis typischerweise akzeptabel

• Fürκ 10⁴ ist das Ergebnis oft unbrauchbar. Das Problem istschlecht konditioniert

• Falls das Ergebnis eines numerischen Problemsunstetig von den Eingabedaten abhängt, wird es i. A. alsunsachgemäßgestellt angesehen.

• Die Kondition eines Problems hängt i. A. von der EingabemengeE ab (siehe obiges Beispiel zur Schnittpunktberechnung)

3.3.2 Beispiele

Addition/Subtraktion

• F(a, b) = a+b.

y Eingabevektorx= (a, b)^T.

• Wähle 1-Norm:k(a, b)^Tk₁ =|a|+|b|

• κ_abs=kF⁰(a, b)k₁ =k(1,1)^Tk₁ = 1. (Matrixnorm!)

• κ_rel= ^|a|+|b|_|a+b|

y Für die Addition von Zahlen mit gleichem Vorzeichen ergibt sich alsoκ_rel= 1 y DieSubtraktionist allerdings schlecht konditioniert, wenn

|a+b| |a|+|b|, (3.10)

also wenn|a|und|b|annähernd gleich sind.

33

Zahlenbeispiel:

– Rechnung auf fünf dezimale Stellen:

x= 1.2345·10⁻³ y= 1.7899·10⁻³

u= 1.0000·10⁰+x= 1.0012·10⁰ v = 1.0000·10⁰+y = 1.0018·10⁰ z1 =v−u= 6.000·10⁻⁴

z₂ =y−x= 5.554·10⁻⁴ – Mathematisch istz1 =z2.

– Numerisch kommt es zum Weghebephänomen und ein erheblicher Rundungsfehler tritt auf

• Subtraktion von annähernd gleichen Zahlen sollte vermieden werden!

Quadratische Gleichung

• Die Lösung vonx²−2px+q = 0mitq < p² istF±=p±p

p²−q. Kondition vonF±?

• Wähle wieder die 1-Norm:k(p, q)k₁ =|p|+|q|

• kF⁰(p, q)k₁:

∂F

∂p = 1± p

pp²−q (3.11)

∂F

∂q = −1

2p

p² −q (3.12)

ykF⁰(p, q)k1 = max 1± p

pp²−q, −1 2p

p²−q

!

(3.13)

• Also:

κ= max 1± p

pp²−q, −1 2p

p²−q

!

·|p|+|q|

|F±| (3.14)

• Problem 1: FürF± →0(alsoq→0) geht die Kondition gegen∞.

• Dann:p

p²−q→ |p|undF− →p− |p|, was zu großen Rundungsfehlern führt.

• Problem 2: Fürq→p² ist es analog.

Abhilfe zu Problem 1

• Satz von Vieta:

F₁ =p+ sign(p)·p

p²−q (3.15)

F₂ = q

F₁ (3.16)

(3.17)

• F_1/2 ist je nach VorzeichenF₊oderF−.

• Dadurch wird die fehlerbehaftete Subtraktion vermieden.

3.3.3 Einfache Funktion

• Kondition vonF(x) = ^1−cos(x)_x : κ=

x·F⁰(x) F(x)

=

x· x

1−cos(x) ·xsin(x) + cos(x) x²

=

xsin(x) + cos(x) 1−cos(x)

(3.18)

• Starker Fehler fürx→0.

• Besser: In diesen Fällen Reihenentwicklung:

1−cos(x)

x = 1

x

1−

1−x² 2 + x⁴

24±. . .

= x 2 ·

1− x²

12±. . .

. (3.19)

• Fürx= 10⁻⁶ ist ^x₁₂² <10⁻¹³und^x/2entsprechend auf 12 Stellen genau.

• Die Kondition wäre fürF(x)nach Gl.(3.18)2·10¹²und der Fehler groß:

1 >> x = 10^(−6);

2 >> F1 = (1−cos(x) ) /x

3 F1 = 5.00044450291171e−07

4 >> F2 = x/2

5 F2 = 5.00000000000000e−07

6 >> F3 = 4.9999999999995833333∗10^−7; %"Exakt" berechnet mit Mathematica

7 >> (F1 −F3)/F3 % relativer Fehler für die Funktion

8 ans = 8.89005824242853e−05

9 >> (F2 −F3)/F3 % relativer Fehler für die Reihenentwicklung

10 ans = 8.32209870677419e−14

35