Einf¨uhrung in die Numerische Mathematik f¨ur Chemiker

(1)

1. Auflage, Stand vom 26. April 2021

Henrik R. Larsson, Bernd Hartke Theoretische Chemie

Institut f¨ ur Physikalische Chemie Max-Eyth-Straße 2

Erdgeschoß, Raum 29 Tel.: 0431/880-2753

hartke@pctc.uni-kiel.de

http://ravel.pctc.uni-kiel.de

(2)

Einleitung

Ziele der Veranstaltung

• erstes Grundverst¨andnis von numerischer im Gegensatz zu analytischer Mathematik

• Erstellung und Verwendung eigener Programme statt fertiger “apps”

• Verwendung von professionellen numerischen Bibliotheken

• Hintergrundwissen f¨ur vorgefertigte Programme der angewandten Mathematik und numerischen Simulation in Physik, Chemie,. . .

KEINE Ziele der Veranstaltung

• detaillierte Einf¨uhrung 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

(3)

Veranstaltungsformat

• mehrere thematische Kapitel; zu jedem Kapitel gibt es:

– eine Vorlesungsstunde: Einf¨uhrung in die theoretischen Grundlagen;

Besprechung des zugeh¨origen Aufgabenblatts – eine betreute ¨Ubungsstunde

• alle Skriptseiten und alle Aufgabenbl¨atter von meiner homepage als pdf-Dateien herunterladbar:

http://ravel.pctc.uni-kiel.de → teaching

→ Numerical Mathematics for Chemists

(Achtung: Inhalt kann sich w¨ahrend des Semesters ¨andern)

• die Benotung erfolgt ¨uber die Bearbeitung der ¨Ubungsaufgaben:

– abzugeben sind:

∗ mindestens: lauff¨ahiger Programmtext, vorzugsweise in Fortran

∗ keine ausf¨uhrbaren Programme

∗ als pdf-Datei: eigener Programm-output, Antworten auf Fragen, ggf. weitere Erl¨auterungen (z.B. bei Schwierigkeiten), etc.

– Funktionierende Programme sind anzustreben.

In Ausnahmef¨allen reicht es, wenn f¨ur mich aus dem Programmtext klar erkennbar ist, daß die Prinzipien verstanden wurden.

– Abgabe der Aufgaben am besten per e-mail, an:

spenke@pctc.uni-kiel.de

• bestanden mit ≥ 50% der maximal m¨oglichen Punkte (angegeben bei jeder (Teil)Aufgabe)

(4)

Selber programmieren!

• Lehrinhalte = numerische Standardalgorithmen

⇒ fertige Programme in fast allen Programmiersprachen online auffindbar.

• ABER: nur selber programmieren bringt Lerneffekte!

• Abgabe von online gefundenen Programmen als Eigenleistung ist ein Plagiat

• . . . es sei denn, die Quelle wird angegeben.

Dann gibt es aber keine Programmierpunkte.

Plagiate

Die Programmieraufgaben sind die (einzige) Prüfungsleistung dieses Moduls, also formal zu behandeln wie z.B. Abschlußklausu- ren. Daher gilt dafür die Prüfungsverfahrensordnung (PVO).

• Laut PVO sind Plagiate schlimme T¨auschungsversuche, direkt gleichzusetzen mit einer wiederholten(!) leichteren T¨auschung (z.B. Benutzung unerlaubter Hilfsmittel).

• Täuschung bei einer (Teil-)Prüfung führt zwingend zu einer Wertung von “nicht ausreichend” oder “nicht bestanden” für diese (Teil-)Prüfung.

• T¨auschungen & Plagiate sollen immer dem Pr¨ufungsausschuß gemeldet werden; damit sind und bleiben sie aktenkundig.

• Der T¨ater ist anzuh¨oren, aber:

• In schwerwiegenden Fällen (u.a. bei Plagiaten) kann der PrüfA dann entscheiden, daß die Prüfung/das Modul als endgültig nicht bestanden gilt.

(5)

Wann ist ein Plagiat ein Plagiat?

Einordnung als Plagiat nicht immer eindeutig oder unstrittig.

Definitiv erlaubt (keine Quellenangabe n¨otig) ist z.B.:

• sich online und/oder in Lehrb¨uchern Pseudocode oder sogar Beispielprogramme anzuschauen, aber dann unabh¨angig davon eine eigene Implementation zu schreiben;

• eigenen Programmcode (z.B. aus fr¨uheren Aufgaben) wiederzuverwerten.

Ohne Quellenangabe nicht akzeptabel ist es z.B.:

• ein fremdes Beispielprogramm oder signifikante Teile davon unver¨andert als Eigenleistung abzugeben,

• oder nur unwesentlichen Kleinigkeiten darin zu ¨andern (z.B. Variablennamen).

Um daraus eine Eigenleistung zu machen, sind größere Änderungen nötig. Dann wäre es leichter gewesen, das Programm von Anfang an selbst zu schreiben.

Wenn Sie sich unsicher sind, ob das Verwenden einer bestimmten Quelle als Plagiat gewertet werden k¨onnte, wenden Sie sich im Vorfeld an den betreuenden Assistenten oder Professor. Die meisten unabsichtlichen Plagiate entstehen durch nachl¨assiges/feh- lendes Kennzeichnen der verwendeten Quellen.

Wir k¨onnen nicht alle T¨auschungsversuche immer entdecken, behalten uns aber vor,

• Verdachtsf¨allen nachzugehen

• und n¨otigenfalls die obigen Konsequenzen folgen zu lassen.

(6)

Inhaltsplan

1. Einf¨uhrung in grundlegende Programmierkonzepte: (Aufg.1) (Octave/Matlab, Fortran)

2. Zahlendarstellung, Fehler, Konditionszahl 3. Differentiation (Aufg.2)

4. Integration: Trapezregel, Simpson, Gauß; (Aufg.3) Anwendung: ¨uberall

5. Zufallszahlen und Monte-Carlo-(MC)-Integration (Aufg.3) Anwendung: MC-Simulationen von Gasen, Fl¨ussigkeiten,. . .

6. Suche nach Nullstellen und Extremwerten: Newton et al. (Aufg.4)

Anwendung: Lösung analytisch unlösbarer Gleichungen; zahlreiche Optimierungsprobleme 7. gewöhnliche Differentialgleichungen (Euler, Runge-Kutta) (Aufg.5)

Anwendung: z.B. einzelne Gleichungen der Kinetik 8. DGL-Systeme und partielle DGLs (Aufg.5)

Anwendung: komplizierte Kinetiken, klassische Mechanik, Quantenmechanik, W¨armeleitung, Elektrodynamik, Akustik, Str¨omungsdynamik,. . .

9. L¨osung linearer Gleichungssysteme (inkl. Matrixinversion); (Aufg.6) Anwendung: ¨uberall, z.B. Regression (

”Fit“) und Interpolation

10. Spline-Interpolation; lineare Regression, χ₂-Statistik; nichtlineare Regression (Levenberg-Marquardt); (Aufg.7) Anwendung: z.B. Anpassung von Modellfunktionen an Meßdaten

11. Matrixdiagonalisierung: Eigenwertproblem; (Aufg.8)

Anwendung: Quantenmechanik, stationäre Lösungen der Schrödingergleichung

(7)

Literaturempfehlungen:

• W. H. Press, B. P. Flannery, S. A. Teukolsky und W. T. Vetterling:

”Numerical Recipes“, Cambridge University Press, Cambridge, 1990:

sehr gut geschriebene Einf¨uhrung 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:

bew¨ahrtes Standardwerk deutscher Mathe-Fakult¨aten zur Numerischen Mathematik; klare Darstellung, aber etwas “mathe- matischerer Stil”.

• M. Hermann: “Numerische Mathematik”, Oldenburg, 2006:

gut zug¨angliches Buch, leider ohne Differentialgleichungen und Optimierung.

• Peter Deuflhard und Andreas Hohmann: “Numerische Mathematik 1,2,3: Eine algorithmisch orientierte Einf¨uhrung”, de- Gruyter, Berlin, 4. Auflage, 2008:

Band 1 bietet eine allgemeine Einf¨uhrung in diverse Themen, Band 2/3 besch¨aftigen sich mit Differentialgleichungen; klarer Text, leider keine Behandlung von Optimierung.

• 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¨anglich.

• G. H. Golub und C. F. van Loan:

”Matrix Computations“, Johns Hopkins University Press, 1989, 1993:

Klassiker zur linearen Algebra; ziemlich

”mathematischer“ Stil; f¨ur Detailinformationen 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 selbsterkl¨arend.

(8)

1 Minimaleinf¨ uhrung ins Programmieren

1.1 Fortran

siehe separates Fortran-Einf¨uhrungsskript 1.2 Octave/Matlab

siehe Octave-Einf¨uhrung in Henrik Larssons Numerikskript WS15/16

(9)

1.42 Welche Programmiersprache soll ich nehmen?

Das ist weitestgehend egal, u.a.:

• alle Numerikalgorithmen sind programmiersprachenunabh¨angig;

• die grundlegenden Sprachstrukturen (Felder, Schleifen, Verzweigungen) gibt es ¨uberall in sehr ¨ahnlicher Form.

Andere Entscheidungshilfen:

• kleines Hilfsprogramm vs. großes software-Projekt

• eigene Entwicklung f¨ur Eigenbedarf vs. kollaboratives multi-user-Paket

• Syntax:

– langatmig aber gut lesbar vs. kurz aber kryptisch – frei = versteckte Komplexit¨^? at

• interpretiert vs. compiliert

• virtual machine vs. OS/hardware-spezifische compiler/interpreter

• Verf¨ugbarkeit ausgefeilter Entwicklungsumgebungen (Editor mit Syntax-highlighting und Inhaltshilfen, integrierte debug- ger/profiler, integrierte Versionsverwaltung, . . . )

• Verf¨ugbarkeit professioneller Bibliotheken (BLAS, LAPACK, MKL, etc.; daran h¨angt letztlich die Laufzeiteffizienz und wieviel code man selbst schreiben muß)

(10)

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.

Grunds¨atzliche Folgen f¨ur die numerische Praxis:

• begrenzte Anzahl Stellen,

• es gibt eine betragsm¨aßig gr¨oßte darstellbare Zahl,

• miteinander verrechnete Zahlen sollten ¨ahnlich groß sein.

Zahlendarstellung im Computer ist bin¨ar, aber mit weiteren technischen Tricks.

Alle Details: siehe Skript WS15/16 (Henrik Larsson)

2.1.1 Ganze Zahlen (integers)

# 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

(11)

Durch die Zweierkomplementdarstellung (siehe Skript WS15/16 Henrik Larsson) sind vorzeichenbehaftete ganze Zahlen “seltsam angeordnet”:

Bin¨arwert 8-bit-Zweierkomplement

00000000 0

00000001 1

· · · · · ·

01111110 126

01111111 127

10000000 -128

10000001 -127

· · · · · ·

11111110 -2

11111111 -1

Was bei “ ¨Uberlauf” des darstellbaren Zahlenbereichs passiert, h¨angt 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 #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;

(12)

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

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

17 }

f¨uhrt nach der ¨Ubersetzung in Maschinensprache zur Ausgabe von a = 127

b = 255

a + 1 = -128 b + 1 = 0

a + 1000 = 103 b + 1000 = 231

Ob und wie das detektiert bzw. verhindert werden kann, hängt sehr stark von der Programmiersprache und dem Compiler ab. Schwierig ist auch, daß es nach dem Überlauf i.d.R. für Vorsichtsmaßnahmen zu spät ist. Machbar und portabel können aber schon einfache Lösungen wie diese sein:

– Tests der Art if (a+1 > a) ...

– tempor¨are Verwendung von real-Variablen (s.u.; Vorsicht: ggf. Genauigkeitsverlust!)

– tempor¨are Verwendung von integer-Variablen mit mehr bits.

(13)

• das Verhalten von octave/matlab ist anders:

– keine Warnung bei ¨Uberlauf

– stattdessen wird der gr¨oßte/kleinste m¨ogliche Wert verwendet

1 >> a = int8(127) % ganze Zahl mit einer Darstellung von 8 Bits

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

Ob das “besser” ist, ist Geschmackssache. . .

Da ganze Zahlen (direkt oder indirekt) h¨aufig als “Zeiger” auf Speicheradressen verwendet werden, ist integer-overflow ein Einfallstor f¨ur 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ß f¨ur die verwendete signed-integer-Darstellung geworden war.

• 1996 st¨urzte 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)

Bin¨are Darstellung einer reellen Zahl:

(−1)^v ·M ·2^e (1)

mit Vorzeichenbit v, m-stelliger Mantisse M und Exponent e.

(14)

• F¨ur x 6= 0 wird e so gew¨ahlt, 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 f¨ur die normalisierte Form (e < e_min), so kommt es zumExponentenunterlauf (underflow) in den subnormalen Zahlenbereich. Dort sinkt die Darstellungsgenauigkeit stetig, bis alle Ziffern 0 sind.

• spezielle Zeichen/Formen:

Exponent Mantisse

normalisierte Form e ∈ [emin+ 1, emax−1] M ∈ [2⁻¹,1−2^−m] denormalisierte Form e = e_min+ 1 M ∈ [2^−m,2⁻¹ −2^−m]

±0 e = e_min 0

±inf e = e_max 0

NaN e = e_max+ 1 6= 0

• NaN (not a number) ergibt sich u.a. bei Division durch Null

• ±inf tritt bei einem Exponenten¨uberlauf (overflow) auf. Die Zahl ist zu groß, um dargestellt zu werden.

• Die meisten reellen Zahlen lassen sich nicht exakt darstellen:

1 > c a t p r i n t r e a l s . f 9 0

2 program p r i n t r e a l s

3 i m p l i c i t none

4 w r i t e(∗, ’ ( a , f 2 2 . 2 0 ) ’) ’ 0 . 3 i n s i n g l e−p r e c i s i o n : ’ , 0 . 3

5 w r i t e(∗, ’ ( a , f 2 2 . 2 0 ) ’) ’ 0 . 3 i n double−p r e c i s i o n : ’ , 0 . 3 d0

6 w r i t e(∗, ’ ( a , f 2 2 . 2 0 ) ’) ’ 0 . 5 i s e x a c t : ’ , 0 . 5

7 end program p r i n t r e a l s

8 > . / p r i n t r e a l s

9 0 . 3 i n s i n g l e−p r e c i s i o n : 0 . 3 0 0 0 0 0 0 1 1 9 2 0 9 2 8 9 5 5 0 8

10 0 . 3 i n double−p r e c i s i o n : 0 . 2 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 8 8 9 0

11 0 . 5 i s e x a c t : 0 . 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

• Dies f¨uhrt zu entsprechenden Rundungsfehlern:

(15)

1 > c a t r o u n d e r r . f 9 0

2 program r o u n d e r r

4 i f ( 1 . 0−0 . 8−0 . 2 == 0 . 0 ) t h e n

5 w r i t e(∗,∗) ’ e x a c t l y z e r o ’

6 e l s e

7 w r i t e(∗,∗) ’ not e x a c t l y z e r o , but ’ , 1 . 0−0 . 8−0 . 2

8 end i f

9 end program r o u n d e r r

10 > . / r o u n d e r r

11 not e x a c t l y z e r o , but −1.4901161E−08

Deshalb sollte man bei real-Zahlen nie auf exakte Gleichheit testen:

1 r e a l : : a

2 i f ( a == 2 . 5 ) . . . ! f a l s c h

3 i f (abs( a−2.5) < 1 . d−7) . . . ! b e s s e r

• Aufgrund der halbexponentiellen Darstellung Gl. 1 sind die exakt darstellbaren Zahlen logarithmisch auf der Zahlenskala verteilt. Beispiel f¨ur M bestehend aus drei Bits und e∈ [−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

(16)

• Maschinengenauigkeit (“Maschinenepsilon”): maximaler relativer Abstand zweier benachbarter, normalisierter real-Zahlen (m Mantissenstellen, Unterschied von 1 in der letzten Stelle):

2^−m

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

• Programm zur Berechnung von :

1 > c a t c a l c e p s . f 9 0

2 program c a l c e p s

4 r e a l(k i n d=8) : : e p s

5 e p s =1. d0

6 do

7 e p s =0.5 d0∗e p s

8 i f ( .not. 1 . d0+e p s > 1 . d0 ) e x i t

9 end do

10 w r i t e(∗,∗) ’ double−p r e c i s i o n machine e p s i l o n : ’

11 w r i t e(∗,∗) ’ c a l c u l a t e d : ’, e p s

12 w r i t e(∗,∗) ’ a n a l y t i c : ’ , 2 . d0∗∗(−53)

13 end program c a l c e p s

14 > . / c a l c e p s

15 double−p r e c i s i o n machine e p s i l o n :

16 c a l c u l a t e d : 1 . 1 1 0 2 2 3 0 2 4 6 2 5 1 5 7E−016

17 a n a l y t i c : 1 . 1 1 0 2 2 3 0 2 4 6 2 5 1 5 7E−016

(17)

Gebr¨auchliche 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⁻³⁸ −3.4·10³⁸ – Wertebereich der denormalisierten Zahlen: 1.4·10⁻⁴⁵−1.17·10⁻³⁸ – Genauigkeit: = 2⁻²⁴ = 6·10⁻⁸

• 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⁻³⁰⁸ −1.797·10³⁰⁸ – Wertebereich der denormalisierten Zahlen: 3.9·10⁻³²⁴ −2.22·10⁻³⁰⁸ – Genauigkeit: = 2⁻⁵³ = 1.1·10⁻¹⁶

(18)

Auch bei reals sind schon gravierende Praxisfehler aufgetreten:

• 1992 wurde bei der Schleswig-Holsteiner Landtagswahl aufgrund eines Rundungsfehlers fälschlicherweise angenommen, daß die Partei Die Grünen die 5 %-Hürde überschritten hat.

• Im ersten Golfkrieg hat eine Patriot-Rakete 28 US-Soldaten getötet, da die Zeit in Inkrementen von ¹/¹⁰ gezählt wurde. ¹/¹⁰ hat allerdings keine endliche binä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

– Molek¨ul in L¨osung versus Molekulardynamiksimulation in der Gasphase mit klassischer Beschreibung der Atomkerndy- namik und approximativer quantenmechanischer Beschreibung der Elektronen

Datenfehler

• Abweichung der realen Parameterwerte von den eingegebenen Daten

• Beispiele:

– Signifikante Stellen von Naturkonstanten – Messungenauigkeiten im Experiment

(19)

Verfahrensfehler

• Intrinsischer Fehler des numerischen Verfahrens

• Beispiele:

– Endliche Taylorreihenentwicklung – Diskretisierungsfehler

• 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 Vielteilchenst¨orungstheorie 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 n¨utzliche Unterscheidungen:

f(x) sei das mathematische Problem in Abh¨angigkeit von der Eingabe x, ˜f der numerische Algorithmus und ˜x fehlerbehaftete Eingabedaten:

Kondition: ||f(˜x)−f(x)||: Wie stark schwankt das Problem bei St¨orung?

Stabilit¨at: ||f˜(˜x)−f˜(x)||: Wie stark schwankt der numerische Algorithmus bei St¨orung?

Konsistenz: ||f˜(x)−f(x)||: Wie gut l¨ost der numerische Algorithmus (mit exakter Eingabe) das Problem?

Konvergenz: ||f˜(˜x)−f(x)||: Wie gut l¨ost der numerische Algorithmus (mit gest¨orter Eingabe) das Problem?

(20)

2.3 Normen

2.3.1 Vektoren

Die Norm k · k in einem Vektorraum X erf¨ullt die Bedingungen 1. Definitheit: k~xk ≥ 0, ~x ∈ X; insbes. k~xk = 0 ⇒~x =~0 2. Homogenit¨at: kα~xk = |α|k~xk, α ∈ C

3. Dreiecksungleichung, Subadditivit¨at: k~x+~yk ≤ k~xk+k~yk

Dies ist für zahlreiche unterschiedliche Normvorschriften gegeben, für X = Rⁿ (gewöhnliche Vektoren) z.B.:

1-Norm k~xk1 =Pn i |xi|

Euklidische Norm k~xk₂ =pP

i|x_i|² p-Norm k~xk_p =pP^p

i|x_i|^p

Maximums-Norm k~xk_∞ = max_i|x_i|

Diese Vielfalt st¨ort i.d.R. nicht, weil sich alle Normen gegeneinander absch¨atzen lassen:

C1k~xk_α ≤ k~xk_β ≤ C2k~xk_α (3)

Die Konstanten C₁, C₂ lassen sich in den meisten F¨allen genau angeben.

2.3.2 Matrizen, Abbildungen

Matrixnormen können über Vektornormen definiert werden oder können auch andere Charakteristika von Matrizen involvieren, die bei Vektoren nicht auftreten (z.B. Eigenwerte, Singulärwerte).

Sinnvollerweise sollten Matrixnormen k · k mit Vektornormen k · k_V vertr¨aglich (kompatibel) sein; das ist dann der Fall, wenn

kA~xk_V ≤ kAk k~xk_V (4)

(21)

Normen und Eigenwerte hängen jedoch auch zusammen: Ist eine Matrixnorm mit irgendeiner Vektornorm verträglich, dann gilt für jeden Eigenwert λ einer quadratischen Matrix A

|λ| ≤ kAk (5)

weil mit dem zugeh¨origen Eigenvektor ~x gezeigt werden kann:

|λ| k~xk = kλ~xk= kA~xk ≤ kAk k~xk (6) woraus die Behauptung nach Division durch k~xk folgt. Nach Gl. 5 ist der Spektralradius (Betrag des betragsgr¨oßten Eigenwerts) einer Matrix niemals gr¨oßer als ihre Norm.

Beispiele f¨ur Matrixnormen:

• basierend auf der Operatornorm:

kAk = max

~x6=~0

kA~xk k~xk

= max

k~xk=1kA~xk (7)

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

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

λ_max;λ_max ist der größte Eigenwert von A^†A (größter Singulärwert von A, siehe SVD) inf-Norm kAk∞= maxi∈[1,n]Pn

j=1|A_ij|(Zeilensummennorm)

• basierend auf einer vektorartigen Anneinanderreihung aller Matrixelemente:

Frobeniusnorm kAk_F =q Pm

i

Pn j |A_ij|²

Auch Matrixnormen lassen sich nach Gl. 3 gegeneinander absch¨atzen.

(22)

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), f¨uhrt, sondern von einer Eingabemenge E, welche zu einer Resultatsmenge R = f(E) f¨uhrt:

E x

f(x) R

Die Kondition charakterisiert das Verh¨altnis zwischen E und R.

• Sei ~y=F(~x) unter der Annahme von exakten Zahlen (ohne Rundungsfehler), sowie ˜~y=F(˜~x), wobei ˜~x=~x+δx~ eine Zahl mit Rundungsfehler ist.

• Der absolute bzw. relative Fehler sind dann

∆_abs =k~y˜−~yk , ∆_rel= k~y˜−~yk

k~yk (8)

• Die absolute Konditionszahl ist

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

also das kleinstm¨ogliche Verh¨altnis k~y˜−~yk ×(k~x˜−~xk)⁻¹.

• Sofern F(~x) differenzierbar ist, gilt

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

(hier wird die Matrixnorm gebraucht)

• Die relative Konditionszahl ist analog

κ_rel= inf (

α

k~y˜−~yk

k~yk ≤αk~x˜−~xk k~xk

)

, (11)

das kleinstm¨ogliche Verh¨altnis^k^~^y−~^˜ ^yk/k~yk×(^k^˜^~^x−~^xk/k~xk)⁻¹.

• Sofern F(~x) differenzierbar ist, gilt

κrel = kF⁰(~x)k

kF(~x)k × k~xk=κabs× k~xk

kF(~x)k. (12)

Im Folgenden: Beschr¨ankung auf κ≡κ_rel

(23)

• Ist κ-1, so ist das Problem gut konditioniert.

– F¨ur κ <1 gibt es sogar eine Fehlerverringerung – F¨ur κ= 1 treten nur Rundungsfehler auf

• F¨urκ=O(100) ist das Ergebnis typischerweise akzeptabel

• F¨urκ10⁴ ist das Ergebnis oft unbrauchbar. Das Problem istschlecht konditioniert

• Falls das Ergebnis eines numerischen Problems unstetig von den Eingabedaten abh¨angt, wird es als unsachgem¨aß gestellt angesehen.

Die Kondition eines Problems h¨angt nicht nur vom Problem, sondern auch von der Eingabemenge E ab.

Beispiele:

• Schnittpunktsberechnung zweier Geraden:

– gut konditioniert:

g g˜

h h˜

r r˜

Der Schnittpunkt ˜r variiert hier ¨ahnlich stark wie die Geraden g und ˜g selbst.

– schlecht konditioniert (schleifender Schnitt):

1

r

Starke Schnittpunktsvariation

(24)

• Addition/Subtraktion:

– F_add(a, b) =a+b bzw. F_sub(a, b) = a−b

– mit Eingabevektor ~x = (a, b)^T; w¨ahle 1-Norm: k(a, b)^Tk₁ = |a|+|b|

– κ_abs = kF⁰(a, b)k₁ = k(1,1)^Tk₁ = 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|, (13)

also wenn |a| und |b| ann¨ahernd gleich sind.

Zahlenbeispiel: Rechnung mit 5 Dezimalstellen:

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

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

z₂ = y −x = 5.554·10⁻⁴

Mathematisch ist z₁ = z₂. Numerisch kommt es zum Weghebeph¨anomen, und ein erheblicher Rundungsfehler tritt auf.

Subtraktion von ann¨ahernd gleichen Zahlen sollte vermieden werden!

• weitere Beispiele siehe Henrik Larssons Numerikskript WS15/16.

(25)

2.5 Zeitaufwand f¨ur Rechenoperationen und Programmstrukturen Faustregeln f¨ur Operationen mit reals:

• Addition/Subtraktion ist sehr schnell

• Multiplikation ist zwischen 50% und Faktor 4 langsamer

• Division ist etwa doppelt so aufwendig wie Multiplikation (daher kann x*0.5 schneller sein als x/2.0)

• praktisch alle Funktionen (exp, ln, sqrt, sin, cos, aber auch x^y) sind mindestens 10mal so aufwendig wie eine Multiplikation (Grund: sie werden über kurze, Taylor-artige Reihenentwicklungen berechnet; manchmal möglicher Ausweg: table-lookup) Faustregeln für schnell ausführbare Programme:

• geschachtelte Schleifen ¨uber mehrdimensionale Felder: innerste Schleife immer ¨uber den “schnellsten” Index (kann je nach Programmiersprache der erste oder der letzte sein)

• keine bzw. m¨oglichst wenige if-Verzweigungen in Schleifen.

(26)

3 Differentiation

Die erste Ableitung einer Funktion f(x) nach x an der Stelle x₀ ist definiert als f⁰(x)|_x=x₀ = df

dx x=x0

= lim

x→x0

f(x)−f(x0)

x−x₀ = lim

∆x→0

∆f

∆x x=x0

(14) Zu dem analytischen Problem, daß es sich hier um einen unbestimmten Ausdruck vom Typ “⁰₀” handelt, kommt aus numerischer Sicht das Problem hinzu, daß mit kleiner werdendem ∆x die Beträge der Zahlenpaare x, x₀ und f(x), f(x₀) immer ähnlicher werden, während sich die Größenordnung dieser vier Zahlen selbst nicht ändert.

Algorithmenunabh¨angige Folgerungen:

• Numerische Differentiation verstößt zwangsläufig gegen die Regel, Subtraktion annähernd gleicher Zahlen zu vermeiden!

• Der Limes in Gl. 14 ist numerisch (bei begrenzter Zahlendarstellung) nicht ausf¨uhrbar.

Abgebrochene Taylorreihen (s.u.) beheben das nicht, erm¨oglichen aber

• Approximationen an Ableitungen, ohne Ausf¨uhrung des Limes (finite Differenzen),

• und mit einer Absch¨atzung des Fehlers.

Daher gilt:

• numerische Ableitung einfacher: braucht nur f(x), nicht f⁰(x), ABER:

• wenn irgend m¨oglich, analytisch ableiten bzw. die analytische Ableitung im Programm verwenden;

• numerische Ableitung nur als Notl¨osung;

• analytische Ableitung mit numerischer Ableitung überprüfbar, aber Erfahrung in der Wahl der Schrittweite h = ∆x nötig!

(27)

3.1 Vorw¨artsdifferenzen

Taylorentwicklung von f(x) um x₀ (mit ∆x = h):

f(x₀ +h) =f(x₀) + hf⁰(x₀) + h²

2 f⁰⁰(x₀) +O h³f⁽³⁾(x₀)

(15)

⇒ f(x₀ +h)−f(x₀)

h = f⁰(x₀) + 1

2hf⁰⁰(x₀) +O(h²f⁽³⁾(x₀)) (16) Also gilt mit einem Fehler linear in h (also vgl.weise schlechte Konvergenz: halb so großes h halbiert den Fehler):

f⁰(x0) ≈ f(x₀ +h)−f(x₀)

h (17)

3.2 Zentraldifferenzen

f(x_o +h) = f(x₀) +hf⁰(x₀) + h²

2 f⁰⁰(x₀) + h³

3!f⁰⁰⁰(x₀) +· · · (18) f(x_o−h) = f(x₀)−hf⁰(x₀) + h²

2 f⁰⁰(x₀)− h³

3!f⁰⁰⁰(x₀) +· · · (19)

⇒f⁰(x₀) = f(x_o +h)−f(x_o−h)

2h − h²

3!f⁰⁰⁰(x₀) + · · · (20) Fehler nur noch quadratisch in h (halb so großes h ⇒ ein Viertel des Fehlers). Veranschaulichung:

¯ x

¯

x−h x¯+h m1

m2

m3 f

Unterschied zwischen Vorwärts- (gepunktete Linie), Rückwärts-(durchgezogene Linie) und Zentraldifferenzen (gestrichelte Linie).

(28)

3.3 H¨ohere Ordnungen

Auf ¨ahnliche Weise sind herleitbar:

• Formeln h¨oherer Ordnungen f¨ur die erste Ableitung,

• Formeln unterschiedlicher Ordnungen f¨ur h¨ohere Ableitungen.

Siehe z.B. https://en.wikipedia.org/wiki/Finite_difference_coefficient.

Alternativer Weg zu Formeln h¨oherer 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⁰(x0) = D(h) +e2h² +e4h⁴ +. . . (21) D(h) ist der Differenzenquotient und e_i die von h unabh¨angigen Fehleranteile.

• Idee der Extrapolation: Berechne f⁰(x₀) f¨ur Vielfache von h:

II : f⁰(x₀) = D(2h) + 4e₂h² + 16e₄h⁴ +. . . (22) 4×I −II

3 : f⁰(x₀) = 4D(h)−D(2h)

3 −4e₄h⁴ +. . . (23)

⇔f⁰(x₀) = −f(x0 + 2h) + 8f(x+h)−8f(x−h) +f(x−2h)

12h + O(h⁴) (24)

Der e₂-term fällt weg (ähnlich wie oben bei Zentraldifferenzen), dadurch steigt die Fehlerordnung nochmals. Diese Methode kann automatisiert werden: Beliebige Genauigkeit und Fehlerabschätzung möglich.

(29)

3.4 Wahl der Schrittweite

Achtung: h und/oder x₀ +h ggf. nicht exakt bin¨ar darstellbar. Abhilfe:

1 temp = x0 + h;

2 h = temp − x0;

Das erhaltene h ist nun “exakt” in der Bin¨ardarstellung (aufgrund des Fehlers in der Subtraktion).

h kann sowohl zu groß als auch zu klein gew¨ahlt werden. Beispiel Vorw¨artsdifferenzen:

• Der Rundungsfehler ist ungef¨ahr ( = Maschinenepsilon) _r ∼

f(x₀) h

(25) (Der Rundungsfehler ist generell proportional zu 1/h, auch f¨ur andere finite Differenzenformeln.)

• Fehler aufgrund von Trunkierung der Taylorreihe f¨ur Vorw¨artsdifferenzen:

_t ∼ |hf⁰⁰(x)| (26)

Beachte: _r ∼ _h¹ aber _t ∼h !

• Minimierung von (h) = _r+_t liefert

h ∼ s

f(x0)

f⁰⁰(x₀) (27)

Die optimale Schrittweite h¨angt von der Maschinengenauigkeit und dem Inversen der relativen Kr¨ummung ^f_f⁰⁰₍^(x_x⁰⁾

0) ab.

(30)

h kann sowohl zu groß als auch zu klein gew¨ahlt werden.

Fehler

Zeitschritt Schrittw

h

(31)

3.5 Imagin¨arschrittableitung

Tatsächlich kann man Differenzen (und damit Rundungsfehler bzw. Weghebephänomene) zumindest bei ersten Ableitungen mit einem Trick völlig vermeiden: Für einen rein imaginären Schritt ih lautet die Taylorreihe:

f(x+ih) =f(x) +ihf⁰(x)− h²

2 f⁰⁰(x)− ih³

3! f⁰⁰⁰(x) +· · · (28)

Division durch h und Betrachtung nur des Imagin¨arteils auf beiden Seiten liefert:

f⁰(x) = Im[f(x+ ih)]

h + O(h²) (29)

Vorteile der Ableitungsformel Gl. 29:

• Fehlerordnung h² (statt h¹ wie bei simplen Vorw¨artsdifferenzen), bzw. selbe Fehlerordnung wie bei Zentraldifferenzen, aber mit nur einer Funktionsberechnung;

• keine differenzbasierten Rundungsfehler/Weghebung, damit sehr viel kleinere h-Werte m¨oglich;

• bei Sprachen mit “eingebauten” komplexen Variablen und Funktionen (wie z.B. Fortran) ist die Implementation sehr einfach.

(32)

3.6 Schlußbemerkungen

• H¨ohere Ordnung bewirkt nicht immer h¨ohere Genauigkeit.

• H¨ohere Ordnung bedeutet nicht unbedingt bessere Effizienz.

• Funktionsberechnungen k¨onnen teuer sein.

Beispiel: Gradienten/Frequenzen in der ab-initio-Quantenchemie.

Goldene Praxisregel: keine finiten Differenzen verwenden, sonderen analytische Ableitungen. Diese m¨ussen jedoch zus¨atzlich 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 Gradientenvek- tor

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 Geome- trieoptimierung ca. 1 h.

• Es existieren weitere Alternativen zur Ableitungsberechnung, z.B. via Fouriertransformation:

dⁿ

dxⁿ f(x) = 1

√2π

∞

Z

−∞

F(k) dⁿ

dxⁿ e^ikxdk = 1

√2π

∞

Z

−∞

F(k) (ik)ⁿe^ikxdk = f⁽ⁿ⁾(x) (30)

(33)

4 Integration

4.1 Vorbemerkungen:

Wann brauchen wir numerische Integration?

• nicht bei analytisch integrierbaren Funktionen, sondern

• bei analytisch gegebenen, aber nicht analytisch integrierbaren Funktionen,

• bei nur numerisch gegebenen Funktionen.

Einzige Voraussetzung an zu integrierende Funktion f(x):

f(x) muß berechenbar sein, f¨ur einige diskrete, vorgegebene Werte von x innerhalb des Integrationsintervalls [a, b].

Grundprinzip:

Zb

a

f(x)dx =

N

X

i=1

c_if(x_i) + Fehler(N, . . .) (31)

• unabh¨angig von der Form von f(),

• Lage der St¨utzstellen x_i bestimmt Werte der Koeffizienten c_i,

• N m¨oglichst klein,

• Fehler für gegebenes N möglichst klein und abschätzbar, soll für wachsendes N schnell gegen Null gehen.

Alle Integrationsverfahren dieses Typs machen implizite Annahmen ¨uber den Verlauf von f(x) zwischen den Integrationsst¨utzstellen x_i. Das ist

• einerseits ein Grundproblem bei allen numerischen Algorithmen, die mit Diskretisierungen arbeiten;

• andererseits theoretisch und praktisch nicht ganz so gef¨ahrlich wie es scheint.

(34)

4.2 Trapezregel (Euler-Verfahren)

(bei ¨aquidistanten St¨utzstellen: Newton-Cotes-Verfahren)

a x b

i x

i+1

h

b

Z

a

f(x)dx ≈ hfi|^b_a(b−a) (32)

= h

N

X

i=1

hfi|^x_xⁱ⁺¹

i (33)

= h

N

X

i=1

1

2(f_i+f_i+1) (34)

= h 1

2f₁+f₂+· · ·+f_N−1+1 2f_N

(35)

In jedem Intervall der Breite h wird f(x) durch Gerade (Polynom 1.Ordnung) approximiert. Daraus folgt:

• die Formel ist exakt nur dann, wenn f(x) eine Gerade ist;

• Wenn wir f(x_i+1) durch eine Taylorreihe um f(x_i) approximieren:

f(x_i+1) = f(x_i) + hf⁰(x_i) + h²

2!f⁰⁰(x_i) +· · · (36) brechen wir hier also nach dem linearen Glied ab

⇒ der Fehler ist von der Ordnung O 1/N²

= O(h²) Kontrolle der Konvergenz in der Praxis:

• einmalige Berechnung des Integrals (mit einem geratenen Wert für h) ist gefährlich; nur mit weiteren Informationen zu rechtfertigen (z.B. hinreichend genaue Resultate mit diesem h-Wert bei sehr ähnlichen Funktionen); sonst immer nötig:

• mehrmalige Berechnung des Integrals mit jeweils kleiner werdenden h-Werten, bis ¨Anderung der Integralapproximation kleiner als vorgegebene Schranke .

(36)

4.3 Programmiertips zur Integration

Ein typisches Integrationsprogramm hat drei wesentliche Teile:

1. Funktion, die f(x) f¨ur einen beliebigen input-Wert x berechnet;

2. Subroutine f¨ur die eigentliche Integration (bei festem N):

• Parameter:

– input: Integrationsgrenzen a, b

– input: (feste) Anzahl von St¨utzstellen N – output: Integralwert s

• berechne h aus a, b, N

• setze s = ¹₂f₁ + ¹₂f_N

• Schleife ¨uber alle i = 2, N −1:

s= s+f_i (wobei f_i ein Funktionsaufruf ist)

• setze s = s×h 3. Hauptprogramm

• user interface

• Schleife mit Z¨ahlvariable i:

– berechne Integral s_i (Aufruf der subroutine) – Konvergenztest: ist |s_i −s_i−1| < ϑs_i ?

wenn ja, Ergebnisausgabe und Stop; sonst weiter – (zum Testen s_i ausgeben)

– erh¨ohe N

Fehleranf¨alligste Stelle: Konvergenztest:

• Absolutwert der Differenz beachten!

• relativer Test besser als absoluter |s_i −s_i−1| < ϑ (dann ist ϑ abh¨angig von f(x) und Integrationsin- tervall)

• die Computergenauigkeit ist begrenzt

⇒ man kann ϑ auch zu klein w¨ahlen!

(37)

4.4 Simpsonregel und andere Verfahren

• Simpson/Keplersche Fassregel: f(x) wird im Intervall h nicht durch Gerade, sondern durch Parabel angen¨ahert; dies ergibt:

Zb

a

f(x)dx ≈ h

3(f1 + 4f2 + 2f3 + 4f4 +· · ·+ 2f_N₋₂ + 4f_N−1 +f_N) (37)

https://de.wikipedia.org/wiki/Datei:Composite_Simpsons_rule2.png

• Implementation sehr ¨ahnlich zur Trapezregel. Achtung: Zahl der Punkte muss ungerade sein!

• Diese Formel ist exakt f¨ur Polynome 3. Ordnung; der Fehler ist von h¨oherer Ordnung als bei der Trapezregel:

O 1

N⁴

= O(h⁴) (38)

• Zahlreiche weitere Formeln dieser Art verf¨ugbar: Bis zur 7. Ordnung (https://de.wikipedia.org/wiki/Newton-Cotes-Formeln).

Danach aufgrund von negativen Gewichten nicht mehr sinnvoll (Rundungsfehler bei Differenzbildung ¨ahnlich großer Zahlen!)

• In manchen Situationen wichtig (s.u.): Es gibt auch sog.

”offene“ Formeln, bei denen die Funktionswerte f(a) und/oder f(b) an den Intervallgrenzen nicht berechnet werden m¨ussen, sondern nur Funktionswerte an Punkten streng innerhalb des Intervalls.

(38)

4.5 Adaptive Quadratur

• Manche Bereiche von f(x) ben¨otigen mehr Integrationspunkte als andere.

• Mehrgitter-Verfahren:

1. Starte bei einem groben Gitter (per Trapezintegration)

2. Verfeinere das Gitter durch Verfahren h¨oherer Ordnung (Trapez→Simpson)

3. Bestimme aus dem Genauigkeitsgewinn vom groberen zum feineren Gitter, wo feinere Intervalle n¨otigt sind.

4. Iteriere bis zur Konvergenz.

⇒ Die Gitterpunkte sammeln sich dort, wo

”viel passiert“.

• Probleme: Fehlerabsch¨atzung, Robustheit

• Beispiele von resultierenden Integrationsgittern:

0 20 40 60 80 100

0 0.2 0.4 0.6 0.8 1

y

x

Funktion Integrationspunkte

0 10 20 30 40 50 60 70 80 90 100

0 0.2 0.4 0.6 0.8 1

y

x

Funktion Integrationspunkte

(39)

4.6 Gauß-Integration

Bei allen obigen Newton-Cotes-Formeln werden nur Anzahl der äquidistanten Stützstellen und Werte der Koeffizienten variiert, aber für Erzielung noch höherer Genauigkeit kann man noch einen weiteren Parameter variieren: Positionierung von nicht-

¨aquidistanten St¨utzstellen. Man kann zeigen (z.B. Stoer/Bulirsch), daß

b

Z

a

f(x)dx =

N

X

i=1

w_if(x_i) (39)

exakt ist f¨ur f(x) = Polynom vom Grad ≤ (2N −1), wenn die x_i Nullstellen eines Orthogonalen Polynoms vom Grad N sind und die w_i geeignet gew¨ahlte Gewichte.

• Exaktheit f¨ur 2N −1 ist viel mehr als bei Newton-Cotes (N −1 bzw. N)

• Theorie dahinter nicht ganz trivial; hier unn¨otig; Anwendung auch ohne m¨oglich.

• {x_i, w_i} aus Tabellen (Abramowicz/Stegun) oder mit speziellen Programmen (Numerical Recipes) f¨ur beliebige N berechenbar (Hintergrund z.T. in Henrik Larssons Numerikskript WS15/16)

• bei bekannten {x_i, w_i} nicht langsamer als Newton-Cotes

• Fehlerformeln existieren (z.B. Stoer/Bulirsch), sind aber schwierig anwendbar (gehen nicht so einfach mit N^k oder h^k).

• Verfeinerung des Integrationsgitters/Mehrgitterverfahren schwieriger als bei Newton-Cotes

• bei Implementation zu beachten:

– wg. Definitionsbereich der Orthogonalen Polynome gelten Gauß-Integralformeln direkt nur im Intervall [−1,+1] ⇒ die St¨utzpunkte m¨ussen auf das allgemeine Intervall [a, b] umskaliert werden, und das Integralergebnis ebenso (mit 1/2×(b−a), analog zur Formel hfi(b−a)).

– Legendre-Polynome sind (un)gerade bzgl. Ursprung, daher sind in Tabellen oft nur halb so viele Punkte und Gewichte wie scheinbar n¨otig angegeben.

(40)

Einschub: Orthogonale Polynome

• Wichtige Klasse von Funktionen, welche ein vollst¨andiges Orthogonalsystem bilden: (klassische) Orthogonale Polynome

• Unterschiedliche Arten je nach Definitionsbereich und Art des Skalarprodukts:

hm|ni ≡ Z _b

a

p_m(x)p_n(x)W(x)dx = h_nδ_m,n = h_n×

(1, m = n,

0, m 6= n. (40)

• W(x) nennt man Wichtungsfunktion.

• Alle orthogonalen Polynome erf¨ullen die folgenden rekursiven Beziehungen:

p₋₁(x) = 0, (41)

p₀(x) = 1, (42)

p_j+1(x) = (x−a_j)p_j(x)−b_jp_j−1(x), j = 0,1,2, . . . (43)

• Die Polynome sind zudem L¨osungen von bestimmten Differentialgleichungen (DGL) und tauchen daher in der Physik sehr oft auf.

• Wichtige Eigenschaft: Ein orth. Polynom vom Grad n hat nverschiedene Nullstellen, welche jeweils zwischen den Nullstellen vom Polynom vom Grade n−1 zu finden sind.

(41)

Die wichtigsten orthogonalen Polynome:

Name [a, b] W(x) Rekursion DGL

Legendre [−1,1] 1 (j+ 1)p_j+1 = (2j + 1)xp_j −jp_j−1 (1−x²)y⁰⁰ −2xy +j(j + 1)y = 0 Laguerre [0,∞) exp(−x) (j+ 1)p_j+1 = (2j + 1−x)p_j −jp_j−1 xy⁰⁰+ (1−x)y⁰ +jy = 0 Hermite (−∞,∞) exp(−x²) p_j+1 = 2xp_j −2jp_j−1 y⁰⁰−2xy⁰ + 2jy = 0 Chebychev [−1,1] √

1−x² p_j+1 = 2xp_j −p_j−1 (1−x²)y⁰⁰ −xy⁰+ j²y = 0 Beispiel Legendrepolynome:

-1 -0.5 0 0.5 1

Legendre Polynome

x

P_n(x)

P₀(x) P₁(x) P₂(x) P₃(x) P₄(x) P₅(x)

p0(x) = 1 p₁(x) =x

p₂(x) = ¹₂(3x² −1) p₃(x) = ¹₂(5x³ −3x)

p₄(x) = ¹₈(35x⁴ −30x² + 3) p₅(x) = ¹₈(63x⁵ −70x³ + 15x)

p₆(x) = ₁₆¹ (231x⁶ −315x⁴ + 105x² −5) p_n(x) = ₂n¹n! × _dx^dⁿn(x² −1)ⁿ

https://de.wikipedia.org/wiki/Datei:Legendrepolynome6.svg

(42)

4.7 Schwierige Integranden

Unstetige Integranden Bei bekannten Sprungstellen: Aufteilung des Integrals, sodass die Sprungstellen jeweils die Integrati- onsgrenzen sind. Unbekannte Sprungstellen auffinden durch adaptive Quadratur.

Nadelimpulse/Singularit¨aten Aufteilung des Integrals wie bei Sprungstellen, dann: Variablentransformation (s.u.), oder: tanh-sinh-Quadratur:

• Trapezregel nach Variablensubstitution x = g(t) = tanh ^π₂ sinh(t)

, dadurch f¨allt der Integrand doppelt-exponentiell ab.

• R1

−1f(x)dx = R∞

−∞f g(t)

g⁰(t)dt ≈P∞

k=−∞ω_kf(x_k),

• x_k = tanh ^π₂ sinh(kh)

, mit Schrittweite h,

• ω_k =

1

2hπcosh(kh) cosh² ^π₂sinh(kh)

• Neben Singularit¨aten an den Grenzen auch gut geeignet f¨ur hochgenaue Integration.

Uneigentliche Integrale

• Benutzung von entsprechender Gauß-Quadratur (Laguerre, Hermite), sofern der Integrand daf¨ur geeignet ist (analytisch).

• |a| = ∞ oder |b| = ∞: Variablentransformation mit x = 1/t, dx = dt/t², also

b

Z

a

f(x)dx =

1/b

Z

1/a

1 t² f

1 t

dt (44)

und dann eine

”offene“ Formel anwenden (ohne Berechnung von f(a) bzw. f(b)), um Singularit¨at bei 1/t² zu vermeiden.

Achtung: Der Integrand muss schneller abfallen als x⁻²!

• Beide Grenzen ∞: Aufteilung in zwei Integrale.

• Oder: Integration ¨uber gen¨ugend großen Bereich falls bekannt ist, dass der Integrand schnell abklingt.

Stark oszillierende Integranden Gefahr von Rundungsfehlern (Subtraktion ¨ahnlich großer Zahlen), evtl. analytische Mittelung ¨uber Teilintervalle

(43)

4.8 Mehrdimensionale Integration

Bei einfachen Grenzen, schwach ver¨anderlichem Integranden und hohem Genauigkeitsanspruch: Wiederholte 1D-Integrationen:

Dabei sind Grenzen der 2., 3., . . . Integration durch Funktionen der vorhergehenden Variablen auszudr¨ucken, z.B. Integral ¨uber einen Kreis mit Radius 1 um den Ursprung:

1

Z

−1

√1−x²

Z

−√ 1−x²

f(x, y)dydx (45)

Problem dieses Zugangs: Anzahl der St¨utzstellen skaliert bei naivem Vorgehen exponentiell mit der Anzahl der Dimensionen:

x x x x

x x x x x

x x x

x x x x x

x x x

x x x x

x x x

x x x x

x x x

x x x x

x x x

(44)

Alternative: D¨unnbesetzte Gitter, z.B. Smolyak- oder Clenshaw-Curtis-Quadratur:

http://people.sc.fsu.edu/∼jburkardt/datasets/sparse_grid_cce/345.png

Hat der Integrand

”lokale Spitzen“ in unbekannten Stellen, ist hochdimensionale Integration hoffnungslos.

(45)

4.9 Einschub: Zufallszahlen

• “echte” Zufallszahlen: Detektion von externen Signalen (Mausbewegung, Spannungs-/Frequenzunterschiede, . . . );

dauert i.d.R. zu lange

• berechnete Pseudozufallszahlen: schell viele Zahlen erzeugbar; aber deterministisch und sich periodisch wiederholend.

• Starte von einer Zahl I₀ (random seed) als und berechne daraus iterativ neue Zufallszahlen. Einfaches Beispiel:

I_j+1 = mod(a×I_j + c, m), (46)

mit Parametern a, c (mod liefert den Rest einer Division zur¨uck). Liefert ramdom integers zwischen 0 und m.

Erzeugung von random reals im Intervall [0,1] ¨uber

r = I_j

m. (47)

• In der Praxis deutlich kompliziertere Formeln. De-facto-Standard: Mersenne-Twister,¹ wiederholt sich erst nach 2¹⁹⁹³⁷ −1≈ 4,3·10⁶⁰⁰¹ Berechnungen.

• Einfacher Test auf Regelm¨aßigkeiten: Graphische Auftragung einer Zahlensequenz als Punkte in einem 2D-plot:

Schlechter versus guter Zufallszahlengenerator:²

1https://en.wikipedia.org/wiki/Mersenne_twister

2http://www.reedbeta.com/blog/2013/01/12/quick-and-easy-gpu-random-numbers-in-d3d11/

(46)

Wichtige Praxistricks:

• f¨ur Fehlersuche: immer wieder denselben random seed verwenden.

• f¨ur Produktionsrechnungen quasi-zuf¨alliger

random seed: aktuelle Zeit, inkl. (Milli)sekunden

Zerst¨orung von Korrelationen und damit Verbesserung eines existierenden Zufallszahlengenerators:

• Initialisierung: generiere N Zufallszahlen und speichere sie in array r(i)

• bei allen folgenden Aufrufen:

– skaliere vorherige Zufallszahl auf integer-Wert j im Intervall [1, N]

– r(j) wird als Zufallszahl-output ausgegeben – generiere eine neue Zufallszahl und speichere sie

in r(j)