Numerik Version 0.6

Vorlesung

Numerische Methoden der Physik

(Modul VNUMP) WS 2011/12 Dozent: Prof. Dr. Marc Wagner
A Umsetzung ins LTEX Alexandros Bouras & Jonathan Enders
Korrektor Martin Stein

V ersion 0.6
(Anmerkungen, Ergnzungen und Fehler bitte melden.) [email protected]
INHALTSVERZEICHNIS
Inhaltsverzeichnis
1 Einleitung
1.1 Algebra 1.1.1 1.1.2 1.1.3
Numerik
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1 1 1 2
Ziel: Computer mit Verstand einsetzen (sonst oft sinnlos). . . . . . . Anwendungsgebiete . . . . . . . . . . . . . . . . . . . . . . . . . . . . Computeralgebrasysteme (Mathematica, Maple, ...) . . . . . . . . . .
Darstellung von Zahlen und Rundungsfehler

2.1 2.2 Ganze Zahlen (integer) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 2.2.1 2.2.2 2.3 2.3.1 2.3.2 Arithmetik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gleitkommazahlen (rationale Zahlen) . . . . . . . . . . . . . . . . . . . . . . Darstellung von Gleitkommazahlen
3
3 3 3 3 4 4 4 5
z = S M 2Ee
. . . . . . . . . . . . . . . . .
Zwei Standardtypen: oat (32 Bits) , double (64 Bits)
Rundungsfehler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Elementare Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerische Ableitung . . . . . . . . . . . . . . . . . . . . . . . . . .
Gewhnliche DGL s & Anfangswertprobleme

3.1 3.2 Physikalische Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Euler - Methode 3.2.1 3.2.2 3.3 3.3.1 3.3.2 3.3.3 3.3.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Zweckmig: Umschreiben auf ein System von Gleichungen 1. Ordnung Lsung durch Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fehelrabschtzung
8
8 8 8 8 9 14 16 17 22
Runge-Kutta-Verfahren
Dynamische Anpassung der Schrittweite . . . . . . . . . . . . . . . . Alternative Fehlerabschtzung / Schrittweitenanpassung: . . . . . . . Einschub: Dimensionslose Gren . . . . . . . . . . . . . . . . . . . .
Nullstellensuche, lsen nicht-linearer Gleichungen

4.1 4.2 4.3 4.4 4.5 Problemstellung, physikalische Motivation Bisektion (fr 4.2.1 . . . . . . . . . . . . . . . . . . .
24
24 24 24 25 26 27
N = 1)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Vor- und Nachteile: . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sekantenverfahren (fr
N = 1)
. . . . . . . . . . . . . . . . . . . . . . . . .
Newton-Raphson-Verfahren (fr Newton-Raphson-Verfahren (fr
Gewhnliche DGL s & Randwertprobleme

5.1 5.2 Physikalische Motivation & Problemstellung . . . . . . . . . . . . . . . . . .
N = 1) N > 1)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
28 28 29 32
Shooting Methode
5.2.1 5.2.2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Beispiel: Quantenmechanik, 1D unendlicher Potentialtopf
Beispiel: Quantenmechanik, 1D harmonischer Osziallator . . . . . . .
Lsung von linearen Gleichungssysteme

6.1 6.2 6.3 6.4 Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Gau-Jordan-Ellimination (direktes Verfahren) 6.2.1 Pivotisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
37 37 38 39 40 40
Gau-Elimination mit Rckwrtssubstitution (direktes Verfahren) . . . . . . LU-Zerlegung (direktes Verfahren) 6.4.1 Crouts Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . .
INHALTSVERZEICHNIS
ii
6.4.2 6.4.3 6.5 6.6 6.7
Lsen von
Ax = b
mittels LU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41 41 41 42 42 42 44 44
Berechnen von det(A) mittels LU . . . . . . . . . . . . . . . . . . . .
QR-Zerlegung (direktes Verfahren)
Iterative Verbesserung einer Lsung . . . . . . . . . . . . . . . . . . . . . . . Methode der konjugierten Gradienten (iteratives Verfahren) 6.7.1 6.7.2 6.7.3 Spezialfall: A sei symmetrisch und positiv denit Verallgemeinerung Konditionszahl einer Matrix, Vorkonditionierung . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Numerische Integration
7.1 Eindimensionale Integrale 7.1.1 7.1.2 7.2 7.2.1 7.2.2 7.2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Newton-Cotes-Formeln . . . . . . . . . . . . . . . . . . . . . . . . . . Gauss-Quadratur (hier Gauss-Legendre-Integration)
45
45 45 46 47 47 47 48
Mehrdimensionale Integration . . . . . . . . . . . . . . . . . . . . . . . . . . Geschachtelte eindimensionale Integrationen . . . . . . . . . . . . . . Monte-Carlo-Integration . . . . . . . . . . . . . . . . . . . . . . . . . Wann ist welches Verfahren geeignet? . . . . . . . . . . . . . . . . .
Eigenwertprobleme
8.1 8.2 8.3 8.4 8.5 Problemstellung, grundlegende Eigenschaften/Tatsachen . . . . . . . . . . . Prinzipielle Arbeitsweise numerischer Eigenwertverfahren . . . . . . . . . . . Jacobi-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Beispiel: fuer phys. Anwendung: kleine Schwingungen . . . . . . . . . . . . . Bibliotheken fuer Eigenwertprobleme . . . . . . . . . . . . . . . . . . . . . .
49
49 50 50 52 54
Interpolation, Extrapolation, Approximation Datenmodellierung

9.1 9.2 9.3 9.4 Polynominterpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Kubische Spline-Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . Methode der kleinsten Fehlerquadrate (least square) . . . . . . . . . . . . .
55
55 56 57 59
2 - Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10 Funktionsminimierung, Optimierung
10.1 Golden Section Search in 1D . . . . . . . . . . . . . . . . . . . . . . . . . . 10.2 Quadratische Interpolation in 1D . . . . . . . . . . . . . . . . . . . . . . . . 10.3 Minimumsuche mit Hilfe von Ableitungen in 1D . . . . . . . . . . . . . . . . 10.4 Simplex-Methode (D 10.5
62
63 64 64 64 66 67 67 69
> 1)
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
D > 1-Minimierung
durch wiederholte
1D
Minimierung
. . . . . . . . . . .
10.6 Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.6.1 Kombinatorische Minimierung . . . . . . . . . . . . . . . . . . . . . . 10.7 Kontinuierliche Minimierung . . . . . . . . . . . . . . . . . . . . . . . . . . .
11 Monte Carlo-Simulation statistischer Zustandssummen

11.1 Ising Modell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Grundlagen der Monte Carlo-Simulation . . . . . . . . . . . . . . . . . . . . 11.2.1 Metropolis-Algorithmen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2.2 Heatbath-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . 11.3 Monte Carlo-Simulation des Ising-Modells
69
69 70 71 71 72
ii
ABBILDUNGSVERZEICHNIS
iii
Abbildungsverzeichnis
2.1
oat:
mit
S = 1 Bit , E = 8 Bits und NM = 23 Bits ist seine Eigenschaften
folgendermaen determiniert.
Stackoverow: Online in Internet: URL: http://stackoverow.com/questions/2795629/cint-oat-double [Stand 15.11.2011]. . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2
double
Wikipedia Artikel Doppelte Genauigkeit: Online in Internet: URL: http://de.wikipedia.org/wiki/Doppelte_Genauigkeit [Stand 15.11.2011]
Euler Methode Oszilattors fr kleine
. . .
4 9 11 12 13 13 14 14 15 16 18 19 20 21 23
3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Runge-Kutta-Verfahren & die Euler-Methode am Beispiel des Harmonischen
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Runge-Kutta-Verfahren & die Euler-Methode am Beispiel des Harmonischen Oszilattors fr weit grere
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fehler der Euler-Methode bei wachsener Zeit & Frequenz . . . . . . . . . . . Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz
Abschtzung a Schrittweite . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.10 Nicht-Harmonischer Oszillator in 1D , V(x)=x , n =2,4,6,....
. . . . . . .
3.11 Runge-Kutta-Verfahrem zweiter Ordnung aufgetragen mit der Algebraischen
Lsung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 Runge-Kutta-Verfahrem dritter Ordnung aufgetragen mit der Algebraischen
Lsung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 Runge-Kutta-Verfahrem vierte Ordnung aufgetragen mit der Algebraischen
Lsung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.14 Harmonsicher Oszillator mit Bercksichtigung der Einheiten (oben) und ein-
mal im dimensionslosen Verhltnis

4.1
x=
x x0 , t =
(unten).
. . . . . . . . .
Sekanten-Verfahren
Quelle: Verndertes de:Bild:Sekantenverfahren Ani.gif, Ursprngliche Version von de:Benutzer:Ralf Pfeifer: URL:http://de.wikipedia.org/wiki/Benutzer:Ralf_Pfeifer
. . . . . . . . . . . .
25
4.2
Ausschnitte einer Animation des Newton-Verfahren 19.11.2011] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 29 30 32 33 34 35 35
Quelle: Ausfhrliche Animation. Online im Internet: URL: http://de.wikipedia.org/wiki/Datei:Newto 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 RK-Entwicklung bis
t = t2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Scannen der mglichen Eigenwerte
Illustration des Shooting Verfahrens fr die niedrigen drei angeregten Zustnde 31 Newton-Raphson Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . Grundzustand des harmonsichen Oszillators . . . . . . . . . . . . . . . . . . Shooting Verfahren: Angewendet fr den Harmonischen Oszillator; Randbedingungen ((x
1) = 0
) numerisch nicht realisierbar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Grobes Scannen der mglichen Energieeigenwerte Grobes Scannen der mglichen Energieeigenwerte
iii
ABBILDUNGSVERZEICHNIS
iv
5.9
(Nicht-nomierte) Wellenfunktionen der niedrigsten vier Energieeigenzustnde. Jeder Zustand entspricht einer Mode der Wellenfunktion. Der Grundzustand als eine viertel Periodenschwingung, der erste angeregte Zustand als eine halbe Periodenschwingung, der zweite angeregte Zustand als eine dreiviertel Periodenschwingung und zuletzt eine ganze Periodenschwingung des dritten angeregten Zustandes dargestellt. . . . . . . . . . . . . . . . . . . . .
36 43 53 56 56 57 59 60 60
6.1 8.1 9.1 9.2 9.3 9.4 9.5 9.6
CG- Verfahren fr
N =2
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
N=10 Massenpunkte (m=Masse) mit Federn (Federkonstante k) verbunden Interpolation von 4 Datenpunkten mit einen Grad-3-Polynom . . . . . . . . Interpolation von 10 Datenpunkten mit einen Grad-9-Polynom . . . . . . . Interpolation von 10 Datenpunkten mit einem kubischen natuerlichen Spline Least squares Approximation von 10 Datenpunkten mit Polynomen von diedrigem Grad. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Least squares Approximation von 4 Datenpunkte mit konstantem Ansatz, also f (x, a) = a . . . . . . . . 2 -Fitting von 4 Datenpunkte . . . . . . . . . . . . . . . . . . . . . . . . . . mit konstantem Ansatz, also
f (x, a) = a
. . .
iv
EINLEITUNG
1 Einleitung
Denition aus Wikipedia:
Die Numerische Mathematik, kurz Numerik genannt, beschftigt sich als Teilgebiet der Mathematik mit der Konstruktion und Analyse von Algorithmen fr kontinuierliche mathematische Probleme. Hauptanwendung ist dabei die approximative Berechnung von Lsungen mit Hilfe von Computern.
1.1
Algebra
Numerik
Heutzutage beschftigen sich Physiker in der Regel nicht mehr rein analytisch mit der Physik, sondern sind auf numerische Anstze und Algorithmen angewiesen. Heutzutage werden Computeralgebrasystemen (CAS) Methoden zum Lsen oder Vereinfachen von algebraischer Aufgaben oder Integrale genutzt. Solche Programme knnen nicht nur mit Zahlen rechnen, sondern auch mit symbolischen Ausdrcken wie Variablen, Polynomen und Matrizen. (Wikipedia )
1.1.1
Ziel: Computer mit Verstand einsetzen (sonst oft sinnlos).

Schreiben von (Simulations-)Programmen z.B.: in fr wissenschaftliche (auch Industrie wie Autobranche) Anwendungen geeigneten Computersprachen Fortran, C, C++,... Dieses wird speziell auf das zept) Numerik impliziert Gleitkommazahlen was Rundungsfehler mit sich bringt. Diese mssen verstanden bzw. in Kenntnis genommen werden, um die Korrektheit der berechneten Ergebnisse auf eine gewisse Stellenanzahl zu gewhrleisten. Solche numerische Verfahren sind oft Rechenintensiv je genauer man ein Problem numerisch bestimmen will und knnen deshalb je nach Komplexitt Stunden, Tage oder sogar Monate andauern. Oft werden Hochleistungrechner bzw. Systeme eingesetzt, um die Rechenzeit zu verringern welche mehrere Prozesse parallel und schneller bearbeiten knnen. Auerdem ist die Wahl der Algorithmen von entscheidender Bedeutung und nicht zuletzt ist auch eine Codeoptimierung zweckmig. In der Regel ist die Herangehensweise eine Mischung aus selbstgeschriebenen Programmteilen und existierenden numerischen Bibliotheken wie GSL, LAPACK, ARPACK, etc.
vorliegende Problem zugeschnitten (kein Universalre-
1.1.2
Anwendungsgebiete
Oft nden sich zu einigen Problemstellungen keine explizite (analytische) Lsung (Dreikrperproblem) oder, im Falle das eine Lsung existiert, der Rechenaufwand zu gro ausfllt. Desweiteren ndet man Anwendungen in
Linerare Gleichungssysteme
EINLEITUNG

1.1.3
Eigenwertprobleme (Mehrdimensionale) Integration Dierentialgleichungen Nullstellensuche, Optimierung.
und vieles mehr.
Computeralgebrasysteme (Mathematica, Maple, ...)
Computeralgebrasysteme (Mathematica, Maple, ...) ergnzen die Numerik
Ein Computeralgebrasystem (CAS) ist ein Computerprogramm, das Methoden der Computeralgebra nutzt. Konkreter kann es Rechenaufgaben aus verschiedenen Bereichen der Mathematik lsen und dabei nicht nur (wie ein Taschenrechner) mit Zahlen, sondern auch mit symbolischen Ausdrcken (Variablen, Funktionen, Polynomen und Matrizen) umgehen. (Wikipedia)
Diese liefern exakte Ergebnisse, welches das Lsen des Gesamtproblems beschleunigen sei es durch das Lsen von Standardintegralen, Vereinfachung von langen Ausdrcken, Koordinatentransformationen (z.B: von kartesisch in Kugelkoordinaten), oder hnliches.
DARSTELLUNG VON ZAHLEN UND RUNDUNGSFEHLER
2 Darstellung von Zahlen und Rundungsfehler

2.1 Ganze Zahlen (integer)
arbeiten (d.h. speichern und rechnen ber gekoppelte Flip Flops bzw. ReBits (als Bezeichnung fr eine Binrzier die blicherweise aus 0 und 1
1
Computer gister) mit besteht).
Integer-Zahl
z = bN 1 ....b3 b2 b1 b0 .
mit
bj
1; 0 z
durch 0 und 1 beschrie-
Diese binre Darstellung, von der positiven Ganze Zahlen ben, nden sie so im Speicher des Computers.
N =1
z =
j=0
bj z j
Typisch ist die Wahl der Bits N = 32 jedoch ndet man auch N=8,16,64,128 je nach Problemstellung und verwendeter Hardware Fr
N = 32
Zahlen (+1 Bit fr Sign).
0zz
2
32
1 = 4.294.967.295
. Analog geht das auch fr die Negativen
2.1.1
Arithmetik
ist in der Regel exakt; Ausnahmen:
Die Arithmetik
Wertebereich wird berschritten Division, Wurzel,... liefern Ganze Zahl d.h. eine Komma Zahl wird abgeschnitten (
7 3
= 2)
2.2
2.2.1
Gleitkommazahlen (rationale Zahlen)

Darstellung von Gleitkommazahlen
gibt das Vorzeichen an
1 Die Zahl 10 wird im Binrsystem z = 1010 dargestellt. 2 Die Arithmetik umfasst vor allem das Rechnen mit den Zahlen, also die Grundrechenarten Addition
(Zusammenzhlen), Subtraktion (Abziehen), Multiplikation (Vervielfachen), Division (Teilen) sowie die zugehrigen Rechengesetze. Zur Arithmetik gehren auch die Gesetze der Teilbarkeit der ganzen Zahlen sowie die Division mit Rest. Weiter zu erwhnen ist das Rechnen mit Brchen. [Quelle: Wikipedia: Aritkel Arithmetik: Online im Internet: URL: http://de.wikipedia.org/wiki/Arithmetik[Stand:01.12.2011]]
S = 1 M=
3
z = S M 2Ee
NM j=0
mj ( 1 )j , m0 = 1 2
mj {0; 1}
Exponent: E ist Integer & e eine Integerkonstante.
3 In der Wissenschaft arbeitet man immer mit normierten Zahlen. Falls eine Binre Zahl normiert ist,
muss die Floating Zahl mit 1,??? .Die Auszeichnende Bit Zier ist also eine 1. In diesem Fall nehmen wir das zur Kenntnis und brauchen diese nicht explizit auszuschreiben was uns auch ein Bit erspart. Dies bezeichnet man auch als Phantombit.
2.2.2
Zwei Standardtypen: oat (32 Bits) , double (64 Bits)
Abbildung 2.1:
oat: mit S = 1 Bit , E = 8 Bits und NM = 23 Bits ist seine Eigenschaften
folgendermaen determiniert.
Stackoverow: Online in Internet: URL: http://stackoverow.com/questions/2795629/c-int-oatdouble [Stand 15.11.2011].
Wertebereich: M=1,1+ ,1+2 ,..,2-2 ,2wobei die
beste relative Genauigkeit charakterisiert und damit
1, 19
107 e = 127
wird.
= ( 1 )23 2
Die kleinste und grte Zahl durch
1038 ...1038 determiniert
2
E(8Bits) e
= 22
8 127
= 2126 ....2+126
Abbildung 2.2:
Wikipedia Artikel Doppelte Genauigkeit: Online in Internet: URL: http://de.wikipedia.org/wiki/Doppelte_Genauigkeit [Stand 15.11.2011]
2. Double
double
=Doppelte Genauigkeit.
Besitzt statt 4 Byte jetzt 8 Byte (64Bit) also dop-
pelte Bitzahl gegenber dem Float. Fr die kleinste bzw. grte Zahl ergibt sich analog zur Rechnung von oben:
2Ee 10308 bis 10308

2.3 Rundungsfehler
Reele Zahl im Computer wird im Allgemeinen wegen endlicher Bitzahl nicht exakt angegeben. von M:
Z = S M 2El
Genauigkeit
107 fr oat4
und
1016
fr
double5
2.3.1
Elementare Beispiele
Extremfall : 1+ = falls
1
38
| |<
4 Gre 4Bytes 5 Gre 8Bytes
fr den Wertebreich : 10 fr den Wertebreich : 10
stimmt das Ergebnis nicht ganz berein (unnormalized numbers).
308
bis 10 & die Genauigkeit liegt bei 10 bis 10 & die Genauigkeit liegt bei 10
38 308
6 15
Subtraktion zweier hnlich groer Zahlen (
Z1 Z2 = Z3 = O(Z1 = O(Z2 1 10 2 10 = 3 10 1.??? 1.??? = O(10n ) O(Z1 )
10n )
10n )
Ergebnis kleine Zahl )
Weil groe Zahlen weniger Nachkommastellen besitzen (Exponent bringt Kommaverschiebung), wird die kleine Zahl nicht mehr genau dargestellt. Stellen lschen sich aus Resultierende Mantisse
6
107n (oat) & 1016n (double) Stellen z.B: 6 abweichen ist nur auf 1 Stelle exakt. um 10
hat nur eine Genauigkeit von
Die ersten n
Dierenz zweier oats, die relativ nur
Beispiel: die Dierenz zweier oats, die relativ nur um
106 abweichen, ist
nur auf 1 Stelle exakt!! Numerische Ableitung

Annahme: f(x) soll bekannt bzw. auswertbar sein, whrend dessen die Ableitung
2.3.2
f(x) = f
(x) es nicht ist (siehe ausfhrlicher in Kapitel 5.2 Shooting Methode).
In diesem Falle knnte man f(x) numerisch durch die Dierenz annhern:
f (x)
Taylorentwicklung
1 f (x) + f (x) h + f 2
(x) h + O(h )
3
(asymmetrisch / schlechtes Verfahren)
f (x)
(symmetrisch / gutes Verfahren)
f (x + h) f (x) + O(h) h
f (x)
Problem: wenn einerseits h klein ist
f (x + h) f (x h) + O(h ) 2h
&
die Dierenz groer Zahlen eine kleine Zahl geliefert ( - andererseits wenn h gro ist, werden die Fehler Also suchen wir ein
f (x + h) f (x)
f (x + h) f (x h)
O(h), O(h )
wird durch
2.3.1). gro!
h = hopt
im mittleren Bereich
Abschtzung
hopt (fr
den asymmetrischen Fall)
Relativer Fehler verursacht durch
O(h):
f (x) f (x)
1f O(h) = 2 f (x) f
(x)h f f
(x)
(x)h (x)
Relativer Fehler verursacht durch
f (x + h) f (x) : f (x)
f (x+h)f (x) h
f (x) f (x)
f (x) f (x)h
6 Als Mantisse bezeichnet man die Ziernstellen einer Gleitkommazahl. Der Teil vor dem e nennt man
Martisse whrend der folgende Teil nach dem e als Exponent bezeichnet wird. - 3.2exp-(5)
wobei und
=107n (oat) bzw. = 1016n (double) betrgt , f(x) ein Zahlenwert ist der Nenner f (x + h) f (x) die Dierenz ist und der Ausdruck dann die rel.
Genauigkeit multipliziert mit der max. Genauigkeit beschreibt. Das Optimum bei Gleichheit (fr den asymmetrischen Fall)
(x) h f (x)
opt
opt
f (x) f (x) hopt
h
f (x) f (x)
f(x) | f (x)
Abschtzung
opt
(x)h f (x)
3
opt
hopt (fr
den symmetrischen Falle)
hopt
,
f (x) |opt f (x)

2 3

h) durchfhren.
In der Praxis: Nicht nur Abschtzung, sondern auch Test der Stabilitt bzgl. der numeri-
schen Parameter(indem Fall hier
Zusammenfassung des Kapitels

Integer Zahl:
z=
N =1 bj z j mit j=0
bj
1; 0
Gleitkommazahl:
NM j=0
z = S M 2Ee wobei S = 1= Vorzeichen ; M = 1 j mj ( 2 ) , m0 = 1 ; mj {0; 1} Exponent: E ist ein Integer & e eine S = 1 Bit , E = 8 = 107 S=1 = 1016
Bit , Bits und
Integerkonstante zum Vorzeichenwechsel im Exponenten. oat: besitzt 4 Byte mit einer Genauigkeit von
NM = 23
Bits ; mit
double: besitzt 8 Byte mit einer Genauigkeit von
E = 11
Bits und
NM = 56
Bits ; mit
Aufgrund von Rundungsfehler sind Manahmen zur optimalen Genauigkeit zu treen - fr den Asymmetrischen Fall:
h h
opt
f (x)hopt f (x) f (x) |opt f (x) opt
f (x) f (x)
- fr den Symmetrischen Fall:
3
2 3
f (x) f (x) |opt
GEWHNLICHE DGLS & ANFANGSWERTPROBLEME
3 Gewhnliche DGLs & Anfangswertprobleme

3.1
Physikalische Motivation
Gekoppelte Newtonsche Bewegungsgleichungen fr N Massenpunkte
mj
mit
j =
(1, ..., N )
und Anfangsbedingungen
rj (t = 0) = rj,0 ; rj (t = 0) = vj,0
mj
rj = Fj (r1 , ..., rn , r1 , ..., rn , t)
Analytisch allgemein unlsbar, z.B: primitives Dreikrperproblem (Sonne, Erde, Mond) Dagegen befassen wir uns erst mit Randwertprobleme (z.B: QM, Schrdingergleichung,
(x1 ) = 0, (x2 ) = 0)
3.2
3.2.1
in Kapitel 5
Euler - Methode
Zweckmig: Umschreiben auf ein System von Gleichungen 1. Ordnung
Um ein Sysem 2ter Ordnung zu lsen, kann man es auf ein Problem der 1ten Ordnung reduzieren mit dem Preis, dass sich die Anzahl der Gleichungen verdoppeln.
y(t + ) =
f (y(t + ), t + )
Newtonsche Bewegungsgleichung quivalent zu
rj = vj
v =
1 Fj (r1 , ..., rn , v1 .., vn , t) mj y

die aus mehreren Funktionen bzw. Vektoren
Deniere eine Funktion/Zeilenvektor besteht:
y(t) = (r1 , ..., rn , v1 ....., vn , t) 1 1 F1 (r1 , ..., rn , r1 , ..., rn , t), ..., FN (r1 , ..., rn , r1 , ..., rn , t)) m1 mN
muss fr beliebig gegebene (y, t) auswertbar sein!
f (y, t) = (v1 , ..., vN ,
3.2.2
Damit gilt:
y(t) = f (y(t), t) ; f
Lsung durch Iteration

Taylorentwicklung
y( + t)
y(t) + y(t) + O( ) = y(t) + f (y(t), t) + O( )
Abbildung 3.1: Euler Methode
Iteration auch mit negativen
Schritten mglich.
Problem: Ungeschickte Diskretisierung hat groen Fehler zur Folge

3.3
Fehler
O( )
pro Schritt
Zeitentwicklung zu t=T
Gesamtfehler
T/ Schritte O( T ) = O( ) schlecht
Nicht bercksichtigt ist eine qualitativ andere Entwicklung nach geringer Zeit ... z.B: chaotisches System (fast gleiche Starbedingungen, aber unterschiedliche Trajektorien)
Runge-Kutta-Verfahren
(Nebenbemerkung: wenn keine Argumente angegeben sind dann sind
y = y(t) , f (y, t) =
f (y(t), t)
gemeint;)
Immernoch betrachten wir unsere Newton sche Problemstellung:
mj
rj = Fj (r1 , ..., rn , r1 , ..., rn , t)
mit r(t = 0) = r0 ; r(t = 0) = r0

Iterative Lsungen wie das Eulerverfahren:y(t + )
O( )
wir uns ein anderes Verfahren an:
T aylorentwicklung
y(t) + f (y(t), t) +
hat sich (siehe Abb.3.1, 3.2, 3.4) auf langer Zeitskala nicht bewhrt. Also schauen
Grundidee des Runge-Kutta-Verfahren ist wie Euler, aber geschicktere Diskretisierung sodass Fehler strker in Second order Runge-Kutta
unterdrckt ist.
k1 = f (y(t), t) =1 Eulerschritt
9
(3.1)
10
1 1 1 k2 = f (y(t) + k1 , t + ) = Eulerschritt in 2 2 2 2 y(t + ) = y(t) + k2 + O( 3 ) 2

3.2 sich mit einem halben Euler-Schritt (t
(3.2)
(3.3)
sodass man bei Gleichung 3.1 sich mit einemganzen Euler-Schritt und bei Gleichung
+ 1 ) 2
dem Ergebnis annhert.
Beweis das aus Gleichung3.2 ,Gleichung3.3 folgt:
1 1 k2 = f (y(t) + 1 k1 t + 1 ) = f + ( f 1 f + f 2 ) + O( ) = f + 2 f + O( ) 2 2 y 2 t 1 y(t + ) = y + f + 2 y + O( ) mit y = f ; y = f
Das RK-Verfahren 2. Ordnung ist nicht eindeutig. Es gibt viele Vorschriften, die zu einem Fehler von
O( )
fhren.
Die Diskretisierung kann bis zu beliebiger Ordnung verbessert werden:
Third order RK
k1 = f (y(t), t) 1 k2 = f (y(t) + 2 k1 t + 1 ) 2 1 k3 = f (y(t) + 1 (k 1 + k2 )t + 2 ) 4 y(t + ) = y(t) + 1 (k1 + k2 + 4k3 ) + O( 4 ) 6

Fourth order RK
7
k1 = f (y(t), t) k2 = f (y(t) + 1 k1 , t + 1 ) 2 2 k3 = f (y(t) + 1 k2 , t + 1 ) 2 2 k 4 = f (y(t) + k3 , t + ) y(t + ) = y(t) + 1 (k1 + 2k 2 + 2k3 + k4 ) + O( 5 ) 6
Wie man dieses Verfahren in C nun implementiert steht dem Leser im Anhang oder unter der
http://th.physik.uni-frankfurt.de/~mwagner/teaching/numerik/RK_.C zur Verfgung.

Hier reicht uns das Resultat, was uns der Algorithmus liefert, graphisch zu betrachten:
7 In der Praxis wird hug das RK- Verfahren 4. Ordnung angewendet. Hohe Ordnungen ermglichen
groe Schritte , jedoch sind durchaus mehrere Operationen pro Schritt ntig.
10
11
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 2.5 2 1.5 1 x/x0 0.5 0 -0.5 -1 -1.5 0 2 4 t 6 8 10 Euler-Methode 2nd order Runge-Kutta 3rd order Runge-Kutta 4th order Runge-Kutta analytische Loesung
Abbildung 3.2: Runge-Kutta-Verfahren & die Euler-Methode am Beispiel des Harmoni-
schen Oszilattors fr kleine
Auf dem Plot sieht man die verschiedenen vorgestellten Verfahren (Euler-Methode und 2nd,3rd & 4th order Runge-Kutta-Verfahren) ber die analytische Lsung des Harmonsichen Oszillators mit der Bewegungsgleichung
x = x
und der Anfangsbedingung
x(t = 0) = x0
aufgetragen. Schon fr geringe
weicht die Euler-Methode gegenber
der analytischen Lsung deutlich ab, dagegen stimmen die anderen Verfahren mit bloem Auge nahezu berein.
11
12
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 2.5 2 1.5 1 x/x0 0.5 0 -0.5 -1 -1.5 90 92 94 t 96 98 100 Euler-Methode 2nd order Runge-Kutta 3rd order Runge-Kutta 4th order Runge-Kutta analytische Loesung
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 2.5 2 1.5 1 x/x0 0.5 0 -0.5 -1 -1.5 990 Euler-Methode 2nd order Runge-Kutta 3rd order Runge-Kutta 4th order Runge-Kutta analytische Loesung
992
994 t
996
998
1000
Abbildung 3.3: Runge-Kutta-Verfahren & die Euler-Methode am Beispiel des Harmoni-
schen Oszilattors fr weit grere

Fr grere
sieht man auch die Abweichung des Runge-Kutta-Verfahrens zweiter
Ordnung, whrend die hheren Ordnungen zur Beschreibung der Oszillation immernoch dem analytischen Ergebnis folgen.
12
13
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 150 Fehler Euler-Methode 100 (xnumerisch - xanalytisch)/x0
50
-50
-100
-150 0 20 40 t 60 80 100
Abbildung 3.4: Fehler der Euler-Methode bei wachsener Zeit & Frequenz
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 0.2 Fehler 2nd order Runge-Kutta 0.15 0.1 (xnumerisch - xanalytisch)/x0 0.05 0 -0.05 -0.1 -0.15 -0.2 0 20 40 t 60 80 100
Abbildung 3.5: Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz Hier werden nun nicht mehr die Amplitude, sondern deren Abweichung bzw. Fehler gegenber der analytischen Lsung fr wachsende
t,
dargestellt.
13
14
HO: V(x) = m2x2/2, ABs x(t=0) = x0, v(t=0) = 0, Schrittweite = 0.1 Fehler 3rd order Runge-Kutta 0.004
(xnumerisch - xanalytisch)/x0
0.002
-0.002
-0.004
20
40 t
60
80
100
Abbildung 3.6: Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz
Abbildung 3.7: Fehler der Runge-Kutta-Verfahren bei wachsender Zeit & Frequenz
3.3.1
Fehelrabschtzung
Fr einen RK Schritt ist es hilfreich und wnschenswert die Fehlerabschtzung zu wissen.
14
15
RK-Schritt mit
y(t) y(t
+ )
2:
dies dient eher zur Fehlerabschtzung als zum Arbeiten. 2RK-Schritte mit

y(t) /2 /2 y(t + ) n+1 mit: n = Ordnung von RK - Verfahren c 2 c ( )n+1 2 2 2
Abbildung 3.8:
Abschtzung a
Damit geschtzter absoluter Fehler
abs = Abschtzung = a
Alternative: relativer Fehler,
| y2 y | 2 2n 1
rel =
Da das Verhltnis kann man
abs | y(t) 2 | 2
1 : 2n 1 von y(t) 2 (t + )& y(t) (t+ ) 2 versuchen y(t + ) durch Extrapolation zu bestimmen.
Fehlerabschtzung ermglicht lokale Extrapolation Korrigiere:
y(t) 2 y(t) 2 + 2 2
y(t) 2 y (t) 2 2n 1
Jedoch ist dies Fragwrdig , weil diesbezglich keine Fehlerabschtzung existiert.
Fehlerabschtzung ermglicht bei gegebener maximaler Fehlerschranke abs,max Abschtzung

8
des grtmgichsten Schrittweite
max :
1 abs, max n+1 ) abs
8 funktioniert auch fr relative Fehler
n+1 max n+1
abs,max abs

15
max
= (
16
3.3.2
Dynamische Anpassung der Schrittweite
1. kleine Schrittweite 2. groe Schrittweite

langsam ungenau
Es muss je nach Problemstellung ein Kompromiss zwischen der Genauigkeit und der Schnelligkeit der Auswertung Individuell getroen werden.
Ideal: Je nach Verlauf der Funktion und dessen Ableitung kann man die Schrittweite
varrieren. Also wenn die Funktion glatt verluft kann man grere Schritte whlen, dagegen bei radikalen Vernderungen , also Schrittgre
f(x)gro, kleine Schritte zu whlen sind: also die
immer neu berechnen lassen.
Abbildung 3.9:
Schrittweite
Eine mgliche Realisierung bietet folgender Algorithmus: 1. Input: Maximal erlaubter Fehler Initale (grobe) Schrittweite Anfangsbedingungen
abs,max
y(t = 0)
wie gehabt
2.
t=0
3. RK-Schritte:
y(t) y (t + ) y(t) =
4. Schtze nun den Fehler ab:
9

2 2
y2 (t + ) 2
abs = Abschtzung = a
5. Passe Schrittweite an:
10
| y2 (t + ) y (t) | 2 2n 1
neu = 0, 9 (
1 abs,max n+1 ) abs
9 Fr welche Norm | y2 (t + ) y (t) | man sich endgltig entscheiden muss oder am besten sich dafr
2
eignet, kann man nicht allgemein festlegen. Ob Euklidische Norm, Maximum-Norm oder hnliches eine geeignete Wahl sei, hngt von der Problemstellung ab und erfordert ein wenig Erfahrung oder Fingerspitzengefhl um sich fr eine zu entscheiden. holungen zu verhindern.
10 Der Wert 0,9 ist ein empirisch bestimmter typischer Wert der dazu dienen soll huge Schrittwieder16
17
6. Falls mit
y(t + ) 2 ( z.B. Ausgabe in Datei) 2 t = t + wird als nchster Schritt ausgefhrt, gehe also zu 3. RK-Schritt mit = neu also mit aktuell geschtzter optimaler Schrittweite. falls jedoch abs abs,max
akzeptiert der Algorithmus unser mssen wir unsere Schrittweite verringern, da wir diese Ungenauigkeit nicht akzep-
abs abs,max
=
3.3.3
2 ).
tieren wollen.
neu (stelle also neue und kleiner Schrittweite ein) Gehe dann wieder zu 3. RK-Schritte
Alternative Fehlerabschtzung / Schrittweitenanpassung:
Verwende RK-Schritte n-ter und n+1-ter Ordnung (statt verschiedene Schrittlngen
und
Hier wre noch zu bemerken das diese Verfahren seit Kapitel 2 nur Grundlegende Techniken
sind und das
bessere und ausgefeiltere
Verfahren existieren.
17
18
AHO: V(x) = xn, n = 2, = 0.5, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 Euler-Methode analytische Loesung 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(x0 /m) t

n-2
10
AHO: V(x) = xn, n = 20, = 1.0, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 Euler-Methode 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(xn-2/m) t 0 6 8 10
Abbildung 3.10: Nicht-Harmonischer Oszillator in 1D , V(x)=x , n =2,4,6,....
Aufgetragen sind die Euler-Methode auf die Analytische Lsung und auf der unteren Graphik die Schrittweite e
aufgetragen. Diese verdeutlicht die Vernderung der Schrittgr-
in Abhngigkeit des Gradienten.
18
19
AHO: V(x) = xn, n = 2, = 0.5, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 2nd order Runge-Kutta analytische Loesung 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(x0 /m) t

n-2
10
AHO: V(x) = xn, n = 20, = 1.0, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 2nd order Runge-Kutta 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(xn-2/m) t 0 6 8 10
Abbildung 3.11: Runge-Kutta-Verfahrem zweiter Ordnung aufgetragen mit der Algebrai-
schen Lsung
Aufgetragen sind das 2nd Order RK-Verfahren auf die Analytische Lsung und auf der unteren Graphik die Schrittweite der Schrittgre
aufgetragen. Diese verdeutlicht die Vernderung
in Abhngigkeit des Gradienten
19
20
AHO: V(x) = xn, n = 2, = 0.5, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 3rd order Runge-Kutta analytische Loesung 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(x0 /m) t

n-2
10
AHO: V(x) = xn, n = 20, = 1.0, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 3rd order Runge-Kutta 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(xn-2/m) t 0 6 8 10
Abbildung 3.12: Runge-Kutta-Verfahrem dritter Ordnung aufgetragen mit der Algebrai-
schen Lsung
Aufgetragen sind das 3rd Order RK-Verfahren auf die Analytische Lsung und auf der unteren Graphik die Schrittweite Schrittgre
aufgetragen. Diese verdeutlicht die Vernderung der
20
21
AHO: V(x) = xn, n = 2, = 0.5, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 4th order Runge-Kutta analytische Loesung 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(x0 /m) t

n-2
10
AHO: V(x) = xn, n = 20, = 1.0, ABs x(t=0) = x0, v(t=0) = 0, abs,max = 0.001, initial = 1.0 2 4th order Runge-Kutta 1.5
0.5 x/x0 0 -0.5 -1 -1.5 0 2 4 sqrt(xn-2/m) t 0 6 8 10
Abbildung 3.13: Runge-Kutta-Verfahrem vierte Ordnung aufgetragen mit der Algebraischen
Lsung
Aufgetragen sind das 4th Order RK-Verfahren auf die Analytische Lsung und auf der unteren Graphik die Schrittweite Schrittgre
aufgetragen. Diese verdeutlicht die Vernderung der
Man erkennt deutlich das fr grere Ordnungen des Runge-Kutta-Verfahren die Genauigkeit der Punkte ansteigt inkl. sich schon bei geringen Vernderungen des Gradienten sich eine Schrittgre
gewhlt wird.
21
22
3.3.4
Einschub: Dimensionslose Gren
Die meisten physikalischen Gren tragen Einheiten (Lngen,Zeiten,Energien,...), jedoch besitzt der Compu ter nicht die Fhigkeit oder das Verstndnis in Einheiten zu denken, sondern rechnet nur mit Zahlen. Welche Mglichkeiten bieten sich an, trotz Einschrnkung zu arbeiten?
Methode 1:
Deniere Einheiten fr den Computer, z.B: alle Lngen in Meter oder der Wert 3,77 im Computer fr Lnge als 3,77m als Nutzer verstehen. Dies zieht selbstverstndlich auch alle Lngen-Werte in dieser Einheit einzutippen mit sich. Alles andere wre Unsinn. Whle typische Einheiten fr die jeweilige Problemstellung, um zu groe bzw. zu klei-
ne Zahlen zu vermeiden (Lichtjahre bzw. parsec in der Astrophysik wogegen Femtometer oder Atomare Einheit (a.u.) in der Teilchenphysik sich anbieten).
Methode 2:
Betrachte stets geeignete dimensionslose Verhltnisse die die Physikalischen Abhngigkeiten & Prozesse charakterisieren wie z.B: die dimensionslose Kennzahlen Mach-Zahl) Aber auch einfache und klassische Beispiele geben den Vorteil der dimensionlosen Verhltnisse preis:
11
(Reynoldszahl,
L=
Klassischer harmonischer Oszillator

1 2 mx
fr die Bewegungsgleichung des H.O.) Messe die Zeit in Einheiten von
1 2 m
x x = x
(Die Masse
m ist also eine bedeutungslose Gre = x
: t = t
d x dt
Auerdem sind Anfangsbedingungen ntig, die z.B:
x(t = 0) = x0 ; x(t = 0) = 0
Messe x in Einheiten von der Amplitude Nun hat unsere Bewegungsgleichung aus
x0 := x = L =
1 2
in Einheiten durch gezielte Wahl von dimensionslose Verhltnisse auf die Form d x x = einerseits bersichtlich reduziert, andererseits oentsichtlich die Mgdt lichkeit geschaen die Physik in den Computer zu bertragen bzw. fr die Numerik relevante Form zum Auswerten der physikalischen Problemstellung gebracht.
mx
x x0
d x x = dt 1 2 m x
x = x
11 Der Vorteil der dimensionslosen Kennzahlen liegt in der Mglichkeit, durch wenige, beispielhafte Messungen im Modellversuch die Lsung fr beliebige andere Flle zu ermitteln, bei denen die dimensionslosen Kennzahlen gleich gro sind wie im Modellversuch.(Quelle: Wikipedia: Dimensionslose Kennzahlen: Online im Internet: URL: http://de.wikipedia.org/wiki/Dimensionslose_Kennzahl [Stand:20.11.2011]
22
23
Abbildung 3.14: Harmonsicher Oszillator mit Bercksichtigung der Einheiten (oben) und
einmal im dimensionslosen Verhltnis
x=
x x0 , t =
(unten).
Zustzlich ist eine solche numerische Rechnung unmittelbar fr unendliche viele Physikalische Parameter (in diesem Fall die freie Wahl der Amplitude auch der (Kreis-)Frequenz
x0
wie
s & x0 s ist T = 2 und ist deshalb allgemein T = 0, 165s oder der gleichen welches sich nur auf ein
verwendbar: jede Periode des Sinus fr beliebige
gltig
nicht wie sonst mit
Beispiel bezieht.
23
NULLSTELLENSUCHE, LSEN NICHT-LINEARER GLEICHUNGEN
24
4 Nullstellensuche, lsen nicht-linearer Gleichungen

4.1 Problemstellung, physikalische Motivation
Gegeben: N nicht-lineare Gleichungen mit N Unbekannten (Dies knnen auch lineare Gleichungen sein, jedoch werden hierfr spter weitaus ezientere Methoden vorgestellt.)
fj (x1 , ..., xN ) = 0 ,
bzw. kompakt
j = 1, ..., N
f (x) = 0
Gesucht:
sind Lsungen
die Nullstellen der Funktion
f (x)
sind.
Hierbei handelt es sich um ein Standardproblem, fr das sich zahlreiche Anwendungen bei physikalischen Problemstellungen & Rechnungen nden. Man bentigt sie z.B. bei Lsung der Schrdingergleichung mit Hilfe von Shooting Methoden (siehe ncshtes Kapitel Randwertproblem).
4.2
Bisektion (fr
N = 1)
x1 , x2 mit f (x1 ) < 0 und f (x2 ) > 0 kann z.B. durch das Plotten der Funktion f (x) gefunden werden. Bisektion: ndet eine Nullstelle von der Funktion f , die zwischen x1 & x2 liegt. Algorithmus: Wir denieren den Mittelwert x = 2 (x1 + x2 ) und betrachten folgende 1
2 Flle: 1. Falls 2. Falls sonst
Startpunkt:
f (x1 ) f () < 0 x2 = x x
sonst
x
=x
x1 x2 hinreichend klein ist, bricht man den Algorithmus12 ab, weil x1 und/oder x2
die Eigenschaft als approximative Nullstelle hinreichend erfllt. nchste Iteration.
Konvergenz: Nun ist es relevant auch den Fehler dieser approximativen Nullstelle
x1 bzw. x2
zu bestimmen und damit ein Ausdruck der Genauigkeit des Verfahrens erhalten.
Der Fehler
der approximativen Nullstelle
x1
wird ber
f (x1 + ) = 0
deniert und
nach n Schritten durch:
1 n | x1 x2 | 2
nach oben beschrnkt.
Aufgrund des Faktors
1n 2 fllt der Fehler exponentiell ab. Alle 3 bis 4 Schritte gewinnt man 1 2 n (n+1 linear in n )
also eine Dezimalstelle an Genauigkeit.
Linear konvergent: bezeichnet man Verfahren fr die z.B. n+1 4.2.1 Vor- und Nachteile:
konvergiert immer (wenn die Funktion glatt verluft)
Bisektion
Lineare Konvergenz stellt sich als schlechtes Verfahren heraus aufgrund der langen Rechenzeit & Auswertung der Funktion (z.B. in der Gittereichtheorie dauert die Auswertung u.U. auf einem/mehreren HPC-Systemen Wochen lang).
12 Ein Algorithmus ist eine aus endlich vielen Schritten bestehende eindeutige Handlungsvorschrift zur
Lsung eines Problems oder einer Klasse von Problemen. [Quelle: Wikipedia: Algorithmus]
24
25
4.3
Sekantenverfahren (fr
N = 1)
Abbildung 4.1: Sekanten-Verfahren

Quelle: Verndertes de:Bild:Sekantenverfahren Ani.gif, Ursprngliche Version von de:Benutzer:Ralf Pfeifer: URL:http://de.wikipedia.org/wiki/Benutzer:Ralf_Pfeifer
Startpunkt: Beliebiger
schen
x1 , x2 mit | f (x2 ) |<| f (x1 ) |. f , die nicht notwendig zwi-
Sekantenverfahren ndet mglicherweise eine Nullstelle von
x1 und x2
liegt.
Grundprinzip13 : Iteration hnelt dem Newton-Verfahren. Man nhert sich approximativ

der Nullstelle indem man den Schnittpunkt der Sekante von als neuen Wert nimmt.
f (x1 ) & f (x2 ) mit der x-Achse
xn+1 xn x = f (xn ) f (xn )
xn xn1 f (xn ) f (xn1 )

zu verweisen, da die Dierenz
Auch hier ist wieder auf Kapitel sehr ungenaues Ergebnis liefert!
2.3.1 Rundungsfehler
fast gleicher Zahlen nahe der Konvergenz zu einem Rundungsfehler fhrt, welches evtl. ein
Algorithmus:
1. 2.
x = f (xxn xn1 ) f (xn ) n )f (xn1 xn+1 = xn + x |x|

hinreichend klein ist
3. Falls
sonst wird die nchste Iteration ausgefhrt und beginnt erneut mit Schritt 1. (siehe Animation: Verweis auf Funote 13)
Abbruch und
xn+1 = approximative N ullstelle
Konvergenz:
13 Ausfhrliche Animation der Iteration

Quelle: Verndertes de::Bild:Sekantenverfahren Ani.gif, Ursprngliche Version von de:Benutzer:Ralf Pfeifer: URL:http://de.wikipedia.org/wiki/Benutzer:Ralf_Pfeifer
1,618... 14 n+1 n ,
also besser als linear, d.h. Bisektion
14 kann man zeigen...
25
26
Vor- und Nachteile:
Gute Konvergenz Konvergiert nicht immer...
4.4
Newton-Raphson-Verfahren (fr
N = 1)
Startpunkt: beliebiges
leitung
15
x1 P1 (xn1 ; f (xn1 )) und P2 (xn ; f (xn ))
Grundprinzip: hnlich wie das Sekantenverfahren, man verwendet aber hier die Ab-
f (xn )
statt die Sekante zwischen
Abbildung 4.2: Ausschnitte einer Animation des Newton-Verfahren Quelle: Ausfhrliche Animation. Online im Internet: URL: http://de.wikipedia.org/wiki/Datei:NewtonIteration_Ani.gif[Stand: 19.11.2011]
Algorithmus:
1. 2.
x = f (xn )
1 f (xn )
xn+1 = xn + x Falls |x| hinreichend

proximative Nullstelle.
klein ist, bricht man den Algorithmus ab;
xn+1
ist dann ap-
Ansonsten fhrt man die nchste Iteration aus und beginnt also erneut mit Punkt 1.
Konvergenz: n+1
f () x 2f () x
2 n
16
man beachte die Quadratische proportionalitt. Das Verfahren konvergiert schneller als das obrige Sekantenverfahren!
Vor- und Nachteile:
Sehr gute Konvergenz Konvergiert auch hier nicht immer
f'
muss bekannt sein oder einfach (= billig) auswertbar sein.
15 Die Ableitung f (x) muss bekannt sein bzw. billig auswertbar sein 16 Eignet sich als Hausbung dies zu berprfen & zu zeigen.
26
27
4.5
Newton-Raphson-Verfahren (fr
N > 1)
Nullstellensuche fr
N =2 N >2
0,
N >1
sehr schwierig! Dazu folgendes Beispiel:
f1 (x1 , x2 ) = 0 ,
gefunden werden.
f2 (x1 , x2 ) = 0 f1 (x1 , x2 ) = 0
und
Es mssen alle Schnittpunkte der Isolinien von [Einfgen der Graphik]
f2 (x1 , x2 ) = 0
Es mssen alle Schnittpunkte der
N 1 Dimensionalen Hyperisochen fj (x1 , ..., xN ) =
j = 1, ..., N
17
gefunden werden. Dies zu berechnen ist nahe zu unmglich.
Dies funktioniert in der Regel ganz gut, sofern man eine ungefhre Vorstellung bzw. guten Schtzwert fr die gesuchte Nullstelle
x1
hat.
Startpunkt: Schtzwert Grundprinzip:
x1 fr
die Nullstelle.
0 = fj (xn + )
mit
Taylor Entwicklung
fj (xn ) + (
fj (xn ))k + O( 2 ) xk
18
(4.1)
0 = f (x ) + (
j n
j = 1, ..., N ; und ( xk fj (xn )) = Jjk (xn ) : 2 Vernachlssige O( )

exakte Korrektur
Jacobi-Matrix
approximative Korrektur
x x = J(xn )1 fj (xn )
xk fj (xn ))k
= fj (xn ) + J(xn )x
bzw.
N = 1-Fall
natrlich enthalten:
x =
Algorithmus:
1. 2.
1 1 f (xn ) = f (xn ) . f (xn ) J(xn )
x = (J(xn ))1 f (xn ) xn+1 = xn + x |x|

hinreichend klein ist, bricht man ab,
3. Falls
xn+1
als approximative Nullstelle.
17 z.B. gewisse physikalische Argumente und hnliches. 18 Die Jacobi-Matrix kann fr N 3 einfach von Hand
als lineares Gleichungssystem gelst werden.
Sonst sind numerische Methoden empfehlenswert, die spter erlutert werden.
27
GEWHNLICHE DGLS & RANDWERTPROBLEME
28
5 Gewhnliche DGLs & Randwertprobleme

5.1
Physikalische Motivation & Problemstellung

Betrachten wir Newton sche Bewegungsgleichungen der Form:
mj rj = Fj (rk , rk , t)
Zeitpunkt
mit Randbedingungen
rj (t1 ) = rj,1 ; rj (t2 ) = rj,2
dies entspricht der Frage, welche Bahn beschreibt ein Massenpunkt, der sich zum
t1
hier und
t2
dort bendet?
Betrachten in der QM die Schrdinger-Gleichung in einer Dimension: d 2m dx (x) + V (x)(x) = E(x) mit Randbedingungen der Form 0 = (x = a) wegen dem Potentialtopf V (|x| a) = Die Energie
(x = a) =
ist im Allgemeinen unbekannt, also ein Eigenwertproblem, und somit
keine Dierentialgleichung. Das Umschreiben auf ein DGL - System ist jedoch einfach: Fasse die Energie
als Funktion (E
= E(x))
auf, jedoch mit der Sondereigenschaft
E (x) = 0
DGL-System fr
(x) , E(x)
(y
(dass die Eigenschaft
E(x) = const
enthlt)
Whle Standardform wie in Kapitel 3:
d dx y(x)
= f (y(x), x)
,f
haben
Komponenten)
Die Randbedingungen lauten:
gj (y(x1 )) = 0 ,
5.2
j = 1, ..., n ;
hj (y(x2 )) = 0 ,
j = 1, .., N n
Shooting
Methode
Grundidee: Whle oder rate die Anfangsbedingungen
gj (y(x1 )) exakt erfllen und x = x2 (siehe Methoden aus

erfllen.
y(x1 ), welche die Randbedingungen nach der Runge-Kutta-Entwicklung der Funktion y nach Kapitel 3) Randbedingungen hj (y(x2 )) = 0 nherungsweise
Finde iterativ mit Methoden aus Kapitel 4 (z.B: Newton-Raphson-Verfahren) immer bessere Anfangsbedingungen numerisch exakt erfllen.
y(x1 ),
welche die Randbedingungen mit
hj (y(x2 )) = 0
am Ende
Beispiel: Mechanik,
m = F (x) x
x(t1 ) = a ,
x(t2 ) = b.
y = (x, v) ,
wre:
(x) f = (v, Fm )
wie in Kapitel 3.
Whle die Anfangsbedingungen
y(t1 ) = (a, )
:
sodass
x(t1 ) = a ,
x(t2 ) b
(ideal
x(t2 ) = b)
erfllt;
RK-Entwicklung von
t = t1
bis
t = t2
28
29
Abbildung 5.1: RK-Entwicklung bis
t = t2
Je nachdem wie gut man die Anfangsbedingungen rt, kommt es zu einer Abweichung der Gre
h(y(t2 )) = x(t2 ) b = 0
Optimiere rens:
z.B. mittels des aus Kapitel 4.5 eingefhrte Newton-Raphson Verfah-
Fasse
h(y(t2 )) = x(t2 ) b
als Funktion von
auf, da
x(t2 )
gerade davon
abhngt. Bestimme numerisch, anhand Kapitel 2.32, die Ableitung
d d h(y(t2 )), welche
fr das Newton-Raphson Verfahren gem Gleichung (4.1) erforderlich ist. Fhre Newton-Raphson Schritt aus, um besseres
zu erhalten:
h(y(t2 )). d d h(y(t2 )). h(y(t2 )) = 0

numerisch exakt erfllt.
Wiederhole Runge-Kutta Entwicklung und Newton-Raphsen Schritt so lange bis gewhltes
die Randbedingungen
5.2.1
Beispiel: Quantenmechanik, 1D unendlicher Potentialtopf
Betrachten wir die Schrdingergleichung mit folgenden Randbedingungen: d 2m dx (x) = E(x) mit den Randbedingungen (x = 0) = 0, (x Jetzt formen wir unsere Gleichung dimensionslos um und erhalten:
= L) = 0.
x = x/L
0) = 0,
d () x d2 x
1 x = 2mL2 E h2 () = E() mit den Randbedingungen ( = x x
( = 1) = 0 x
Zum spteren Vergleich mit der numerischen Lsung fhren wir die analytische Lsung vor:
() = sin(n), n = 1, 2, 3, ... x x E = 2 n2
Numerische Lsung: Forme Eigenwertgleichung (Schrdingergleichung) in ein System 1. Ordnung DGLs um. Somit sind die DGL s mit dem Runge-Kutta-Verfahren behandelbar:
29
30
() = () x x x)() () = E( x x x x E () = 0 damit E() = const.
Shooting mit den Anfangsbedingungen
19
:
20
( = 0) = 0, x
-
( = 0) = # = 0, x
x E( = 0) =?
Grobes Scannen der mglichen Eigenwerte ...
Potentialtopf: Grobes Scannen der moeglichen Energieeigenwerte 1
0.8
0.6 (x=1)
0.4
0.2
-0.2 0 20 40 E 60 80 100
Abbildung 5.2: Scannen der mglichen Eigenwerte
Im Bereich von
0 E 100
liegen 3 Eigenwerte,
E0 10 , E1 40 , E2 90
wie es Abb.
5.2 graphisch als Nullstellen darstellt.

die Startwerte bzw. Anfangsbedungen wurden gut (aus)gewhlt, sodass die numerische

19 Irgendeine Zahl = 0. 20 Wird mit Newton-Raphson Verfahren so eingestellt, dass die Randbedingungen ( = 1) = 0 erfllt x
ist.
Lsung sich sehr Nahe am analytischen Ergebnis bendet .
30
31
Potentialtopf: Iteratives Finden der Wellenfunktion des 2. angeregten Zustands 0.15 Vor dem 1. Newton-Raphson-Schritt Nach dem 1. Newton-Raphson-Schritt
0.1
0.05
-0.05
-0.1 0 0.2 0.4 x 0.6 0.8 1
Potentialtopf: Iteratives Finden der Wellenfunktion des 2. angeregten Zustands 0.05 Vor dem 1. Newton-Raphson-Schritt Nach dem 1. Newton-Raphson-Schritt 0.04
0.03
0.02
0.01
-0.01 0.95
0.96
0.97 x
0.98
0.99
Potentialtopf: (Nicht-normierte) Wellenfunktionen der niedrigsten Energieeigenzustaende 0.4 Grundzustand 1. angeregter Zustand 2. angeregter Zustand
0.3
0.2
0.1
-0.1
-0.2 0 0.2 0.4 x 0.6 0.8 1
Abbildung 5.3: Illustration des Shooting Verfahrens fr die niedrigen drei angeregten Zu-
stnde
31
32
Abbildung 5.4: Newton-Raphson Iteration

Abbildung 5.4 zeigt, dass nach bereits 3 Iterationen die Shooting Methode auf ber 7 Stellen exakt ist.
5.2.2
Beispiel: Quantenmechanik, 1D harmonischer Osziallator
d m 2 2 (x) + x (x) = E(x) 2m dx 2
Dimensionslose Gleichung: Forme Lnge aus
, m, : [] =
1 s
[ ]=
kgm2 s ,
[m] = kg,
d dx 1 d a d x 2
x = x/a = () + x () = x x
d d2 x
Die Randbedingungen lauten:
Lnge
a=
2E x ()
x = E() ( +) = 0 x
21
( ) = 0 , x
Einschub: Paritt
Paritt
P:
Raumspiegelung,
P x = x,
P = (x)
gerade ungerade
Eigenwerte und Eigenwertfunktionen:
P (x) = (x) P P (x) = 2 (x) 2 = 1 = 1 = + bzw.

1
[H, P ] = 0, falls V (x) = V (x),

existiert. dass sie entweder -
P = + : P (x) = (x) = +(x) P = : P (x) = (x) = (x)

also fr spiegelsymmetrisches Potential. Der Kom-
mutator drckt aus, dass ein gemeinsamer Satz von Eigenfunktionen zu
und
quivalent: Wellenfunktion der Schrdingergleichung knnen so gewhlt werden,
P = + : P (x) = (x) = +(x) gerade (0) = 0 P = : P (x) = (x) = (x) ungerade (0) = 0
ist in der Numerik problematisch.
21
erfllen.

32
33
zurck zu den Randbedingungen: Ersetze die Randbedingung
Versuch:
(0) = 0
( ) = 0 x
durch eine dieser beiden Bedingungen
oder
(0) = 0
Ersetze
im klassisch
( +) = 0 durch x = L = 0 mit L hineichend weit x a verbotenen Bereich E V (L) (dort fllt exponentiell ab. siehe
Abbildung 5.5 )
Abbildung 5.5: Grundzustand des harmonsichen Oszillators im Computer:
() = () x x 2 E())() x () = ( x x x () = 0 E x
Shooting mit Anfangsbedingungen
( = 0) = 1 (beliebig, Hauptsache = 0) x ( = 0) = 0 (damit P = + ausgezeichnet) x x E( = 0) =? 22 (wird mit Newton-Raphson gung ( = x L ) = 0 erfllt ist). a
Verfahren so eingestellt, dass Randbedin-
1 muss bereits in der Nhe des gesuchten Zustands liegen; Grundzustand E = 2 E = 1 (bei analytisch unlsbaren Problemen, was die Regel ist, wird der Schtzwert durch grobes scannen ermittelt),
22
also
E = ( = 0) 1 x
whlen, damit shooting diesen Zustand liefert.
33
34
HO: RB (x>>1) = 0 numerisch nicht realisierbar 2 E = 1 106 E = 1 + 106
1.5
0.5
-0.5
-1 0 1 2 3 4 x 5 6 7 8 9
Abbildung 5.6: Shooting Verfahren: Angewendet fr den Harmonischen Oszillator; Rand-
bedingungen ((x
1) = 0
) numerisch nicht realisierbar.
In Abbildung 5.6 sieht man das die Randbedingung steigenden Lsung wird fr groe
x=
L a
= 0
numerisch nur
schlecht realisierbar ist. Selbst eine winzige Beimischung der verbotenen exponentiell an-
Zweckmigeres Vorgehen:
oder oder...
dominieren. Nutze diese Eigenschaft aus, indem man weit im klas-
sischen verbotenen Bereich mit beliebigen Anfangsbedingungen startet. z.B.:
x= x= E x= x= x=
L a L a L a
=1 =0 =? =0 =1
L a L a
Verwende nun das Shooting-Verfahren, um stellen von
( = 0) = 0 bzw. ( = 0) = 0 x x
durch Ein-
zu realisieren.
34
35
HO: Grobes Scannen der moeglichen Energieeigenwerte 1.2e+22 1e+22 8e+21 RB (siehe Legende) 6e+21 4e+21 2e+21 0 -2e+21 -4e+21 -6e+21 -8e+21 0 2 4 6 E 8 10 12 (x=0) (x=0)
Abbildung 5.7: Grobes Scannen der mglichen Energieeigenwerte
HO: Grobes Scannen der moeglichen Energieeigenwerte log(|(x=0)|+1) log(|(x=0)|+1)
50
40 RB (siehe Legende)
30
20
10
0 0 2 4 6 E 8 10 12
Abbildung 5.8: Grobes Scannen der mglichen Energieeigenwerte
35
36
HO: (Nicht-normierte) Wellenfunktionen der niedrigsten Energieeigenzustaende 1 Grundzustand 1. angeregter Zustand 2. angeregter Zustand 3. angeregter Zustand
0.5 / maxy|(y)|
-0.5
-1 0 1 2 x 3 4 5
Abbildung 5.9: (Nicht-nomierte) Wellenfunktionen der niedrigsten vier Energieeigenzustn-
de. Jeder Zustand entspricht einer Mode der Wellenfunktion. Der Grundzustand als eine viertel Periodenschwingung, der erste angeregte Zustand als eine halbe Periodenschwingung, der zweite angeregte Zustand als eine dreiviertel Periodenschwingung und zuletzt eine ganze Periodenschwingung des dritten angeregten Zustandes dargestellt.
36
LSUNG VON LINEAREN GLEICHUNGSSYSTEME
37
6 Lsung von linearen Gleichungssysteme

Numerische Verfahren fr partielle DGL s fhren meistens auf lineare Gleichungssysteme z.B. in der Numerische Elektrodynamik, Strmungsmechanik, Kontinuumsmechanik.
6.1
Problemstellung
Sei A eine quadratische gige Gleichungen.)
N N
Matrix mit
det(A) = 0
(A beschreibt also linear unabhnwre es ratsam in der Lite-
(Fr nicht quadratische Matrizen oder solche mit ratur nachzuschlagen.)
det(A) = 0
bj ,
j = 0, ..., M 1 = Axj = bj A1
N-komponentige Vektoren
Typische Problemstellungen: Lse (eine oder mehrere rechte Seiten).
Berechne
(Niemals zur Lsung von
Axj = bj
gem
xj = A1 bj
einsetzen. Rundungsfehler
werden i.d.R. sehr gro! Nheres dazu wird in Kapitel 7 erklrt).
Zwei Klassen von Verfahren, die im Laufe des Kapitels vorgestellt werden:
Direkte Verfahren:
Fhren in endlich vielen Schritten zur Lsung. Rundungsfehler
Iterative Verfahren:
deren Positionen).
(bei groen N) werden problematisch. Iterative Annherung an die Lsung. Unempndlich gegen
Rundungsfehler, jedoch sehr aufwndig bei dichtbesetzten Matrizen. Sinnvolle Anwendung fr dnn besetze Matrizen. ( Speichert nicht alle Eintrge z.B. Eintrge mit 0-en werden nicht abgespeichert, sondern nur deren Eintrge ungleich 0 und
Welche Grenordnung fr N handhabar? (hngt von dem(n) Rechner(n) und ihren Leistungen ab)
Dicht besetzt: Dnn besetzt:
N = O(1000) N O(106 )
Numerische Probleme :
Rundungsfehler fhren zu linear abhngige Gleichungen (wird vom Algorithmus erkannt). Rundungsfehler knnen grer werden als das gesuchte Ergebnis (hug dann, wenn die Gleichungen nahezu linear abhngig sind. Funktioniert bei linear unabhngigen Gleichungen). Um zu berprfen, ob dies aufgetreten ist, berechnet man hug das Ergebnis mit dem gegebenen
Axj
und vergleicht
bj
6.2
Gau-Jordan-Ellimination (direktes Verfahren)

(Enthlt Problemstellung Berechne
Ziel: Lse
Axj = bj (x0 , x1 , ...xN 1 ))
A1 :
Whle
bj = ej
A
Grundidee: Ziehe Gleichung so lange geeignet voneinander ab bis Lsung vorliegt.
37
38
Schrebweise:
Axj = bj
a00 a01 a02 ... a10 a11 a12 ... a20 a21 a22 ... ....
b00 b01 ... b0M 1 b10 b11 ... b1M 1 b20 b21 ... b2M 1 ....
(xj wird bei expandierter Schreibweise weggelassen, da von Umformungen unbeeinusst).
1.Schritt: Ausrumen der 0. Spalte:
= 0...N 1 (1) (damit a00 = 1) (1) = 0...M 1 b0j = (1) (1) ajk = ajk aj0 a0k , j = 1...N 1, k = 0...N 1 (1) (1) bjk = bjk aj0 b0k , j = 1...N 1, k = 0...M 1 (1) (damit aj0 = 0, j = 1...N 1) (1) (1) (1) (1) b00 b01 ... 1 a01 a02 ... (1) (1) (1) (1) b10 b11 ... 0 a11 a12 ... Axj = bj (1) (1) (1) (1) 0 a21 a22 ... b20 b21 ... a0j = ....
Spalte:
(1)
a0j a00 , j b0j a00 , j
b0M 1 (1) b1M 1 (1) b2M 1
(1)
....
2. Schritt bis Nter (induziert durch n-) Schritt: Ausrumen der 1. bis N-1
0
Die nchsten 4 Zeilen bringen n-1te a-Spalte auf die Form
0 1 .
0 : 0
an1j = an1j /an1n1 , j = n 1...N 1 bn1j = bn1j /an1n1 , j = 0...M 1 ajk = ajk bjk = bjk
23
(n)
(n1)
(n1)
(n)
(n1)
(n1)
(n)
(n1)
ajn1 an1k , j = n q, k = n 1, ...N 1 ajn1 bn1k , j = 1, k = 0, ...M 1

(0) (0) (n1) (n)
(n1) (n)
(n)
(n1)
Gleichgungen enthalten 1.Schritt wenn
ajk := ajk , bjk := bjk b00 b01 ... b0M 1 (N ) (N ) (N ) b10 b11 ... b1M 1 (N ) (N ) (N ) b20 b21 ... b2N 1 ... . . ... x0 ... ... xM 1 x0 , ..., xM 1
(N ) (N ) (N )
1 0 0 ...
0 1 0 ...
0 0 1 ...
... ... ... ...
Ix
= bj
(N )
Damit entsprechen b-Spalten den gesuchten Lsungen
Vor und Nachteile:

Alle
bj
mssen gleichzeitig behandelt werden, sonst inezient
relativ langsames Verfahren, wenn Problemstellung
Axj = bj
wenn
6.2.1
O.K. fr Berechnung von
A1
Pivotisierung
Probleme bei Gau-Jordan:
Voraussetzungen
an1 n1 = 0
(n1)
eventuell nicht erfllt.
23 a(n1)
j,n1n1
=0
38
39
Falls
an1 n1 0
(n1)
ist, dann numerisch hochgradig instabil aufgrund Division durch
kleiner Zahl. (Rundungsfehler sind besonders gro). Anordnung der linearen Gleichungen in der Matrix Zeile in spezieller Weise benutzt, dann 1.Zeile,...?!
A, bj
willkrlich; dennoch wird 0.
Lsung: Ordne Gleichung d.h. Zeilen von A und bj zweckmig um. Teilpivotisierung: Vertausche vor dem n-ten Schritt n-1-te und j-te Zeile (wobie
j=
n 1, .., N 1), wobei sich die j-te Zeile durch maximales |

wie oben.
(n1) aj n1 |auszeichnet; alles andere

immer erfllt (weil
Nach
diesen Zeilenaustausch haben wir
an1 n1 = 0
(n1)
det(A) = 0
[Grundvoraussetzung])
Rundungsfehler
im Allgemeinen stark reduziert.
skalierte Teilpivotisierung:
Gewisse Willkr auch bei Teilpivotisierung enthalten - z.B: sei vor 1.Schritt
| a00 |sehr
(0)
klein (1.Zeile wrde mit anderer vertauscht)
- Multipliziere 1. Zeile mit sehr groer Zahl (verndert Lsung und Numerik)... dann kein Austausch der 0. Zeile aber es wrde i.d.R. zu numerischen Instabilitten fhren bzw. verursachen. Beseitige Willkr durch Ersetzung im Austauschkriterium von dieser Teilpivotisierung:
(n1) an1 n1
| aj,n1 | maximumk=0,...,N 1 | ajk |

(0)
(n1)
(behalte Zeile, die relativ gesehen ein groes Element vorne stehen hat, rume beim Rest die aktuelle Spalte aus).

Pivotisierung ist Notwendig! Es reicht das eine 0 in der Matrix enthalten ist welches die
Rechnung ruiniert. Es existieren weitere Pivotstrategien, z.B: zustzliche Umordnung der Spalten (Vollpi-
votisierung). Aufwand und Gewinn ist je nach Problemstellung abzuwiegen.
6.3
Gau-Elimination mit Rckwrtssubstitution (direktes Verfahren)
hnlich zu Gau-Jordan:
1. bs N-te Schritt ist nahe zu identisch wie beim Gau-Jordan, aber verndere im n-ten Schritt nur Zeilen unterhalb der n-1-ten Zeile, um dort 0-len in der n-1-ten Spalte zu erzeugen (hier nur N-1 Schritte ntig).
a00 a01 (1) 0 a11 (2) 0 0 a22 0 0 0 ... Dreieck von 0 len
(0)
b00 ... ... b0,M 1 (1) b10 ... ... b1,M 1 ... ... ...
(0)
39
40
Dann Rckwrtssubstitution (Berechnen von - zuerst: - dann:
xn ):
xN 1 n =
xN 2 n =
- Allgemein:
xjn =
(N 1) bN 1 n (N 1) aN 1 N 1 (N 2) (N 2) bN 2 n aN 2 n1 xN 1 n (N 2) aN 2 N 2 (j) bjn N 1 ajk xkn k=j+1 wobei (j) ajj
in Reihenfolge j=N-1,N-2,...,1,0 ausge-
werten werden muss!
Vor- und Nachteile:
Alle
bj mssen
gleichzeitig behandelt werden, sonst inezient etwa ein Faktor 1,5 schneller als Gau-Jordan (nicht besser
Frwenige
bj (M<<N)
als LU-Zerlegung was in der Regel eher verwendet wird.)
6.4
LU-Zerlegung (direktes Verfahren)
LU-Zerlegung von Matrix A:
A=LU 1 0 0 10 1 0 L = = 20 21 1 30 31 32 0 0 =lower triangular matrix 0 1
0 = 0 0 =upper triangle matrix 0 0 0 Ax = bj (damit

auch
Ermglicht eziente Lsung von
A1 )
sowie von
det(A)
(s.u.).
6.4.1
Crouts Algorithmus
Zu lsen sind
N2
Gleichungen die Komponenten von der Matrix A:
aij = ik kj bzgl.
ij , ij
Knnen trivial gelst werden, wenn diese in der richtigen Reihenfolge behandelt werden: Fr 1.
j = 0, 1, ..., N 1
(fr alle Spalten) fhre folgende beiden Schritte aus:
ij = aij
k=0
2.
ik kj , i = j + 1, ..., N 1
(6.1)
ij =
1 (ij ji
ik kj ), i = j + 1, ..., N 1 mitij =Pivotisierung ij = 0

sein kann!
(6.2)
Pivotisierung auch hier essentiell, da
Wie in Kapitel 6.2.1, also z.B: Teilpivotisierung oder skalierte Teilpivotisierung
40
41

6.4.2
i=j 1, ..., N 1 .24

Whle bei 6.1 6.2 fr
ideale Zeile, also vertausche j-te Zeile mit k-ter Zeile,
k = j, j +
Eziente Bestimmung der idealen Zeile problemlos, da die Ausdrcke in Gleichung
ij
identisch.
Lsen von
Ax = b mittels LU
Weiterhin betrachten wir Gleichgunen der Art
Ax = b = L U x =: y
1. Schritt: Berechne
y , x
ber Vorwrtssubstitution:
Lse
Ly = b
yj =
bj
j1 k=0 jk yk jj =1
j = 0, 1, ..., N 1.
ber Rckwrtssubstitution
25
2. Schritt: Berechne
xj =
Vorteile:
yj
N 1 k=j+1
Lse
Ux =
jk xk
jj
j = N 1, N 2, .., 1, 0
Nach LU-Zerlegung knnen mehrere
bj
seperat bhandelt werden (nicht so bei Gau-
Jordan,Gau).
Nicht langsamer wie det(A), mglich.

6.4.3
Gau-Jordan, Gau, sowoohl fr
Axj = bj
als auch fr
A1 .
Berechnen von det(A) mittels LU

N 1
det(A) = deta(LU ) = det(L) det(U ) = =1 det(A) = (1)v

6.5
N 1 26 j=0 jj j=0
jj
Bei Pivotisierung muss das Vorzeichen unter Umstnden verndert werden:
v =#der
Zeilenvertauschung
QR-Zerlegung (direktes Verfahren)

QR-Zerlegung von
A = QR27
Q kann trivial, R einfach invertiert werden. Anwednungen sind also hnlich wie bei dem LU-Verfahren bzw. LU-Zerlegung. Doppelete
28
Rechenzeit wie bei LU, dafr numerisch mglicher Weise stabiler.
Details: spter bei Eigenwertproblemen.
41
42
6.6
Iterative Verbesserung einer Lsung
Numerisch bestimmte Lsung Rundungsfehlern. Verbessere Lsung:
eines linearen Gleichungssystems i.A. nicht exakt wegen
z.B: mit LU

Sei
Ax = b = b A(x +
Dann
x ges. Korrektur
)=b
bzw.
Ax = b Ax = b
Lse dieses System zur Verbesserung
x x + x

Verbesserung kann mehrmals angewendet werden (wobei einmal i.d.R. ausreicht).
Sehr zu empfehlen: Kostet bei Verwendung von LU kaum zustzliche Rechenzeit, bringt unter Umstnden groen Genauigkeitsgewinn.
6.7
Methode der konjugierten Gradienten (iteratives Verfahren)
Hier behandeln wir lediglich nur ein Beispiel des konjugierten Gradienten (CG- Verfahren) - Problem:
N N M atrizen mit N 10000 sehr speicherintensiv, (10000)2 8 Byte 1GB = Grenordnung eines Hauptspeichers
double
- dnn besetzte Matrizen (fast alle Eintrge 0) dagegen problemlos speicherbar, speichere nur von 0 verschiedene Eintrge. z.B:
10000 3 8Byte 1M B
-
N = 10000 , Tridiagonalmatrix
vernachlssigbarer Speicherplatz.
- Das bietet jedoch keinen Vorteil bei Verwendung direkter Verfahren (formen dnn besetzte Matrix i.d.R. in dicht besetzte Matrizen um).
Wnschenswert: Verfahren die zur Lsung von
Ax = b
lediglich orginales
A und 6.7.1
(und evtl. einige wenige Hilfsverfahren) bentigen.
Spezialfall: A sei symmetrisch und positiv denit

-
Ziel: Lse
Ax = b
(eine rechte Seite, kein
A1 ,
kein det(A))
- Zu Grunde liegende Idee: Minimiere

29
1 f (x) = 2 xAx bx
f (x) = Ax b = 0
Vorgehen:
24 LU-Zerlegung gehrt somit zu Zeilenvertauschter Matrix A 25 siehe auch Abeschnitt 6.3 26 Vorsicht: in der Regel liefert fr groe N N 1 jj winzige oder riesige Zahlen aus. 27 wobei Q= eine orthogonale Matrix und R=j=0 obere Dreiecksmatrix. eine 28 Daher wird das LU-Verfahren blicher fr Lsungen lineare Gleichungssysteme genutzt
j=0
Rate Lsung:
x0 = # R
(kann beliebig verkehrt sein z.B:
x0 = 0)
Whle eine geeignete Suchrichtung
Pj (siehe
weiter unten).
42
43
Minimiere
f (xj + j Pj )
bzgl.
j .
j+1
= x j + j Pj | b Axj |
hinreichend klein,
Falls
j =j+1
0
Abbruch.
x = x
j+1 = gesuchte Lsung.
Gehe Zum Punkt Whle eine geeignete Suchrichtung Fr Implementation notwendige Details:
Pj
zurck.
r = b Ax P =r
0 0
- j
0 (rj ist zu
xj
gehriges Residuum).
Whrend jeder Iteration
30
rj rj Pj APj
rj+1 = rj j APj j =
rj+1 rj+1 rj rj
Pj+1 = rj+1 + j Pj Pj ), dass xj+1 Minimum von f nicht P0 , P1 , P2 , ..., Pj ist. Lsung von
Verfahren so konstruiert (i. B. Wahl der Suchrichtung nur bzgl.
Pj ,sondern
aller bisherigen Suchrichtungen
Ax = b
nach sptestens N Schritten erreicht.
Rundungsfehler: Kein exaktes Ergebnis nach N Schritten. gewnschte Genauigkeit erreicht ist. CG- Verfahren fr
Iteriere so lange, bis die
N = 2:
Abbildung 6.1: CG- Verfahren fr
N =2
Eine Geometrische Interpretation fr
N > 2 ndet man z.B: im Skript von Professor
30 siehe Numerical Recipies, Kapitel 2.7.6
Hans Grabmller unter der Online URL: http://fauams5.am.uni-erlangen.de/am1/de/scripts_grabmue.
43
44
6.7.2
Verallgemeinerung
Falls A nicht symmetrisch und/oder positiv denit ist:
Biconjugate gradient method Minimum residual method Generalized minimum residual method ...
Hier empehlt sich an weitere Literatur zuwenden, da dieses speziale Thema den Rhamen dieser Vorlesung verlsst.
6.7.3
Jede .
Konditionszahl einer Matrix, Vorkonditionierung
N N M atrix A kann verlegt werden gem A = U diag(0 , 1 , , ..., N 1 )V T
U,V: Orthogonale Matrizen.
j 0
(Singulrwerte).
Konditionszahl:
cond(A) :=
maxj (j ) minj (j ) .
Groe Konditionszahl eignet sich nicht gut fr CG: Einerseits sind viele Iterationen ntig und andererseits ist die Numerische Genauigkeit stark limitiert.
A symmetrische positiv denit
j sind Eigenwerte von A. 1 f (x) = 2 xAx bx = const 1 j
Halbachsen der Ellipsoide
Numerisch problematisch, wenn die Funktion lang und spitz ist (und damit die Konitionszahl ein groen Wert annimmt). Dagegen ist eine sphrische Funktion ideal.
Warnung
Matrix A ist nicht immer symmetrisch, positiv denit. Lse
AT A
x = AT b
statt
Ax = b
z.B: mit CG.
Niemals machen !!!

6.7.2.)
symm., pos. def init Ax = b(siehe Abschnitt
In diesem Falle verwendet man verallgemeinere Methoden fr
cond(AT A) = (cond(A))2
verhlt sich also numerisch viel schlechter.
Hug ist eine Vorkonditionierung hilfreich:
A-1 Ax = A-1 b statt Ax = b , wobei = b leicht lsbar ist (z.B: analytisch) - Ax -1 -1 - A A 1 und damit cond(AA) 1, also A A A.
Lse
numerisch leichter handzuhaben als
44
NUMERISCHE INTEGRATION
45
7 Numerische Integration
7.1
Eindimensionale Integrale
Ziel: Berechne bestimmtes Integral
I=
b
a
dxf (x)
7.1.1
Newton-Cotes-Formeln
Grundprinzip: Naehere Naehere
f (x)
durch Polynom, integriere Polynom dann exakt
f (x)
durch Lagrange-Interpolation
Waehle
x+1
(zugehoerige
xj j = 0, 1, ..., n Funktionswerte fj = f (xj ))

Stuetzstellen
Deniere: lj (x)
:=
xxk k=j xj xk (damit lj (xi )
= ji )
(lj gibt Grad des Polynoms)
Approximiere Integrand: Es gilt wenn

pic
f (x)
g(x) =
n j=0 fj lj (x)
f (xj ) = g(xj )
ist
n=2
g(x)
eine Parabel
Integriere approximierten Integranden
I=
b
a
dxf (x)
b
a
dx
n j=0 fj lj (x)
b n j=0 fj ( a dx lj (x))
n j=0 fj j
j :
Gewichte der Stuetzstellen
Fehlerabschaetzung: (Ableitung siehe Grabmueller-Skript)
*n *n
gerade:
I = |
b
a
dx(f (x) g(x))| =
g Kn (n+2) () (x+2)! f
ungerade:
I =
u Kn (n+1) (), (n+1)! f
[a, b]
mit
g Kn = u Kn =
b ab
a
x x
n j=0 (x xj )dx n j=0 (x xj )dx
(b a)n+2
exakt integriert!
Polynome vom Grad Trapezregel (n
(b a)n+1 n + 1 bzw. n werden
= 1):
x0 = a, x1 = b 1 I = f (x)dx ( 1 f0 + 2 f1 )h 2 b 3 I = 1 a (x a)(x b)dxf () = h f () = O(h3 ) 2 12 = 2):

(aequidistante Stuetzstellen)
Simpson-Regel (n
x0 = a, x1 = a + h, x2 = a + 2h = b b 4 I = a f (x)dx ( 1 f0 + 3 f1 + 1 f2 )h, 3 3
(h
ba 2 )
I = h f (4) () = O(h5 ) -sehr 90 siehe ausgeteilte Kopien!
viel besser fuer
klein!
Wiederhole Trapezregel
45
46
Zerlege Intervall
I=
1 1 ba a f (x)dx ( 2 f0 + f1 + f2 + ... + fN 1 + 2 fN )h = TN h, h = N (Faktoren 1 1 1 2 fehlen da 2 fn + 2 fn = fn bei einem und dem darauolgenden Intervall) 1 I = N O(h3 ) = O( N 2 )
[a, b]
in
gleiche Teile und wende jeweils die Trapezregel an
Verfeinere Teilung schrittweise um Faktor 2, d.h.
N 2N ,
bis hinreichende
Genauigkeit erreicht ist, Fehler wird pro Schritt um Faktor 4 reduziert.
Wiederholte Simpsonregel
Naehere Integral gemaess:
I=
b
a
4 dxf (x) 3 T2N 1 TN ( 4 1 )I = I 3 3 3
Naive Erwartung an den Fehler:

4 3 T2N
1 I = O( N 2 )
Tatsaechlicher Fehler viel kleiner:
1 TN = ( 3 =(
2 3 hf0 1 3 hf0 1 3 f0
+ 4 hf1 + 4 hf2 + 4 hf3 + ... + 4 hf2N 2 + 4 hf2N 1 + 2 hf2N 3 3 3 3 3 3 2 hf2 ... 2 hf2N 2 1 hf2N ) 3 3 3 2 + 4 f1 + 3 f2 + 4 f3 + ... + 2 f2N 2 + 4 f2N 1 + 1 f2N ) h 3 3 3 3 3
7.1.2
Das ist gerade die wiederholte Simpson Regel
Gauss-Quadratur (hier Gauss-Legendre-Integration)

Grundprinzip: Waehle Stuetzstellen nicht aequidistant; Nutze Freiheit bei Stuetzstellenwahl, um den Fehler stark zu unterdruecken Legendre Polynome: maess
pj (x)
Entstehen durch Orthogonalisierung von
x0 , x1 , x2 , ...
ge-
1 dx
pj (x) pk (x) =
2 2j+1 jk ,
p0 = 1, p1 = 1, p2 =
als Nullstellen von
3x2 1 2
waehle Stuetzstellen OBdA den).
x0 , ..., xn1
pn (x),
also
pn (xj ) = 0
a = 1, b = 1
(kann durch lineare Koordinatentransformation erreicht wer-
Approximiere Integral wie in Kapitel 7.1 .
I=
1 f (x)dx
n1 j=0 fj wj
I = O(h2n+1 ) fuer ein Teilstueck I = O(1/N 2n ) bei wiederholter Anwendung 2n 1

werden exakt integriert Beweis: sei dann
Polynome bis zum Grad
* f (x) q(x) r(x) * f (x)dx =
f (x)
Polynom vom Grad
2n 1
f (x) = q(x)pn (x) + r(x) (Restpolynom) darstellbar ( , Polynome vom Grad n 1) 1 1 1 n1 j=0 r(xj )j 1 1 dx(q(x)pn (x) + r(x)) = 1 dx r(x) = q(x)pn (x) faellt weg, da q(x) durch p0 , ..., pn1 darstellbar ist und diese orthogonal zu pn (x) sind! 1 1 f (x)dx = n1 f (xj )j q.e.d. j=0
gemaess
Berechnung der Nullstellen von Kapitel 4
pn
(also der Stuetzstellen) z.B. mit Verfahren aus
Fehlerabschaetzungen nur hilfreich, wenn Ableitungen
f (n)
beschraenkt; sonst wach-
sen die Vorfaktoren des Fehlers staerker als dessen Unterdrueckung mit
hn
oder
1 Nn
46
47
7.2
Mehrdimensionale Integration
Zusaetzliche Schwierigkeiten: Hohe Anzahl an Stuetztellen, wenn im 1d-Fall benoetigt werden, dann fuer ein vergleichbares Integral in
Dimensionen
N Stueck ND
Integrationsbereich kann "komplizierte" Form haben:

pic
7.2.1
Geschachtelte eindimensionale Integrationen

Im Folgenden
D=3
(Integration in z.B.
x, y , z )
und
Bestimme minimales und maximales Bestimme minimales und maximales analog: Damit:
x: x1 y
x2 y1 , y2 )
als Funktion von x (
z1 (x, y), z2 (x, y) z y x I = V dxdydzf (x, y, z) = x12 dx y12 dy z12 dzf (x, y, z) (x, y)
fuer das z-Integral dann y Integral fuer alle
Verwende 1D Methoden fuer alle
x mit 1D Methoden und dann das x-Integral mit 1D-Methoden
soll heissen: berechne Funktionswert/Integralwert und gebe weiter ...
Diese Aufspaltung ist nur fuer einfache Integrationsbereiche. Evtl. Aufspaltung in mehrere Integrale noetig! pic
7.2.2
Monte-Carlo-Integration
statistische Abschaetzung des Integrals, d.h. Ergebnis kommt mit Fehlerbalken (analog zu experimentellen Ergebnissen) waehle zufaellig gleichverteilt N Punkte im Integrationsgebiet
I=
=
mit
V : x1 , x2 , ..., xN
dD xf (x) = V < f > V ( <f

Schaetzwert fuer
2 ><f >2
N 1
)1/2
V < f >:
I zu 68% liegt der tatsaech< f >
2 2 ( <f ><f > )1/2 : Fehlerbalken, "1 -Abweichung" N 1
liche Wet von
<f >
im bestimmten Intervall
< f >=
1 N
1 N
N j=1 f (xj ) N 2 j=1 (f (xj ))
< f 2 >=
Nachteil: konvergiert sehr langsam, da Fehler Vorteil:
1 N
sehr grosse D moeglich (allerdings nur ezient, wenn Integrand "gutartig" (=

nicht stark gepeakt))
komplizierte Integrationsbereiche problemlos, lediglich Funktion g erforderlich,

die entscheidet, ob
g(x) =
waehle
1 0 xj
falls
x in V xV
liegt:
sonst
folgendermassen:
47
48
* * *
7.2.3
Waehle einfaches Volumen Waehle gleichverteilt Akzeptiere pic7
welches
umgibt
als
xj ,
xV falls g(x) = 1,
sonst goto Punkt 2
Wann ist welches Verfahren geeignet?

Wenn hohe Genauigkeit erforderlich:
D 1D (D-fach
geschachtelte
1D
Integration)
scheidet das Monte-Carlo-Verfahren(MC), wegen seiner langsamen Konvergenz
D 1D
pic 8
1 , aus. N
Komplizierter Integrationsbereich einfach fuer MC, schlecht fuer erfordert glatte Integranden (beschraenkte Ableitungen)
D 1D
sonst Fehlerabschaetzungen wertlos; MC erfordert Integranden, die nicht in kleinen Gebieten stark gepeakt sind, sonst kommt es zu riesigen statistischen Fehlern.
Im Falle eines komplizierten Integrationsbereich und bei einem Integrand der nicht stark gepeakt ist,fhrt es zu geringe Genauigkeit. Einfacher Integrationsbereich, Integrand glatt,
MC geeignet.
nicht riesig
D 1D
stark oszillierdende Integranden/ unstetige Integranden
MC
stark gepeakten Integranden: Teile Integrationsgebiet in mehrere Bereiche, in denen jeweils Integrand gleiche Groessenordnung hat
(erfordert Kenntnis der Positionen
der Peaks, wenn nicht bekannt, dann sehr sehr schwer!).
48
EIGENWERTPROBLEME
49
8 Eigenwertprobleme
8.1
Problemstellung, grundlegende Eigenschaften/Tatsachen

Eigenwertproblem:
Avj = j vj (A j I)vj = 0
quadratisch,
Vorgehen:
A,
N N vj = 0), j = 1...N
j , vj :
Eigenwerte (EW) und -vektoren (EV) (wichtig:
Bestimmungsgleichung fuer
j N
(eventuell komplexe) Nullstellen
det(A I) = 0
Polynom vom Grad
N,
hat
EWs und EVs (nicht notwendig verschieden, linear unabhaengig)
Schiften von EWs: Ersetze d.h.
A A = A + I
Dann gilt
A vj = (j + )vj ,
hat gleiche EVs aber um
geschiftete EWs
EWs 0 bzw. sehr kleine/grosse EWs vom mathematischen Standpunkt nicht speshifting haeug angewandt, um Numerik zu verbessern
ziell
Eigenschaften von EWs und EVs:

A A A
reell, symmetrisch (A
= A):
EWs und EVs reell EWs reell
komplex, hermitesch (A
= A):
nicht symm/ hermitesch: EWs und EVs i.A. komplex
wenn
* *
normal (A
= AA ):
hermitesch (z.B. QM) bedingt normal
falls EWs paarweise verschieden falls EWs entartet
EVs orthogonal
EVs koennen orthogonalisiert werden
(z.B. Gram-Schmidt auf "Gruppen entarteter EVs")
Generalisiertes EW-Problem:
Avj = j Bvj
kann ordinaeres EW-Problem umgeformt werden
B 1 Avj = j vj
wenn
* *
A, B
symmetisch,
positiv denit, dann bessere Methode:
Fuehre Chaelesky-Zerlegung aus: lich Wurzelziehen
B = LLT
Literatur; aehnlich LU, aeh-
L linke Dreiecksmatrix LT rechte Dreiecksmatrix T Damit Avj = j LL vj 1 A(LT )1 LT v T L j = j L vj = j vj A vj

(ordinaeres EW-Problem) Vorteil:
ist symmetrisch (im Gegensatz zu
B 1 A)
einfach, weil Dreiecksmatrix
... und dafuer gibt es bessere Verfahren Wenn EVs von
gesucht: loese:
vj = LT vj
49
EIGENWERTPROBLEME
50
8.2
Prinzipielle Arbeitsweise numerischer Eigenwertverfahren

schrittweise Veraenderung von
mittels Aehnlichkeitstransformationen:
1 1 1 A P1 AP1 P2 P1 AP1 P2 ... 1 1 1 ... Pn ...P2 P1 AP1 P2 ...Pn 1 1 1 (Pn ...P2 P1 = Q1 , P1 P2 ...Pn = Q)
Waehle
Pj
moeglich;
Q A
Q1 AQ diagonal ist (in der Regel nicht in endlich vielen Schritten multipliziere so lang mit Pj , Pj1 , bis Q1 AQ nahezu diagonal)
so, dass hat EWs
1 AQ
= diag(1 , 2 , ..., N )
EVs
vj = ej (ej :
Einheitsvektor)
hat ebenfalls EWs
det(A I) = det(Q1 )det(A I)det(Q) = det(Q1 AQ Q1 IQ)

also
und
Q1 AQ Qej ,
haben gleiches "charakteristisches Polynom" und damit gleiche
EWs hat EVs also Spalten von
sind EVs von
A:
Q1 AQvj
= j vj (vj = ej )
AQvj = j Qvj
Zusammenfassung:
Bringe
gemaess
Q1 AQ
auf eine Diagonalform
Lies EWs von
aus Diagonalen von
Q1 AQ
ab berechnet
Falls EVs benoetigt, koennen sie mittels

werden, EVs von
Q = P1 P2 ...Pn
sind Spalten von Q
8.3
Jacobi-Verfahren
Voraussetzung:
reell und symmetrisch
dann auch
normal
reell,
vj
reell und orthogonalisierbar.
Q = v1 v2 ... vN Q1 = QT = diag(1 , 2 , ..., N )

.
v1 v2 ... vN
A, QT AQ =
Eigenvektoren bilden orthogonale Matrix; diese diagonalisieren
Vorteile: sehr einfach, insbesondere fuer kleine EWs sehr genau. Nachteil: Langsamer als z.B. die QR-Methode (um konstanten Faktor aber bedeutender Faktor).
50
EIGENWERTPROBLEME
51
Wdh. 8. Eigenwertprobleme
1 1 1 A P1 AP1 P2 P1 A Q1
bis
P1 P2 ...
Q
(Spalten sind EVs von A) vonA)
diag(1 ,2 ,...n )(EWs
A reel, symmetrisch j reell vj orthogonalisierbar,
koennen reell gewaehlt werden
Verwende Drehmatrizen fuer
Pj
in pq-Ebene, um nach
Verwende fuer
PjT APj
Apq
Pj Rotation (sog. Jacobi-Rotation) = Aqp = 0 zu erzielen: ... . p ... ... q ... .
AA =
1 | | ... c s Pj = ... s c ... ... ... c = cos , s = sin

Akl = (PjT )km Amn (Pj )nl

(Pj )mk

k, l = p, q : Akl = Akl App = (Pj )mp Amn (Pj )np = = c2 App + s2 Aqq 2csApq c wenn n = p s wenn n = q

Aqq = c2 Aqq + s2 App + 2csApq Apq = Aqp = (c2 s2 )Apq + sc(App Aqq ) Akp = Apk = cAkp sAkq (k = p, q ) Akq = Aqk = cAkq sAkp (k = p, q ) Apq = Aqp == 0 t = tan
1 1 2( t ! c2 s2 2sc Aqq App 2Apq
=:
1
sgn()
t)
||+(2 +1) 2
1
t2 + 2t 1 = 0 t = (2 + 1) 2 =
erfahrungsgemaess stabilere Numerik) falls
(waehle kleineres
t;
gibt
sehr gross naehere
t=
1 2
Fuer Numerik ist Umschreibung zweckmaessig
Apq = Aqp = 0 App = App tApq Aqq = Aqq + tApq Akp = Apk = Akp s(Akq + Akp ) (k = p, q ) Akq = Aqk = Akq + s(Akp Akq ) (k = p, q )
Bei Jacobi-Rotation mit
Pj
so, dass
Apq = Aqp = 0 vergroessern sich unter UmstaenP1 P2 ...

ueberhaupt gegen eine Diago-
den andere o-diagonal Elemente Konvergiert eine Folge von Jacobi-Rotationen nalmatrix?
51
EIGENWERTPROBLEME
52
Deniere Abweichung von Diagonalform: Man kann zeigen:
S=
k=l
A2 kl S
gegen 0, wenn "grosse
S = S 2A2 ( kl
damit geht
Elemente" wegrotiert werden) Wie waehlt man
p, q ? Apq
wegrotiert hat (
Jacobi 1846: Immer das groesste o-diagonal Element

viel zu teuer fuer Computer / suche .......)
Heute: zyklische Jabobi-Methode
p.q = 1.2, 1.3, ... ... 1.N 2.3, 2.4 ... 2.N 3.4 ... N 1.N
oberes Dreieck
also in fester Reihenfolge durch z.B.
EVs (falls benoetigt) ueber Spalten von
Q = P1 P2 ...Pn
Initialisiere Dann bei

.......
Q=I
auch
A A = PjT APj
Q Q = QPj
8.4
Beispiel: fuer phys. Anwendung: kleine Schwingungen

Betrachte Federkette: pic (in 1D)
L=
N 1 2 j=1 2 mxj
N 1 1 j=1 2 k(xj
xj+1 )2
(quadratisch in
xj , xj )
BGL linear
m1 = +k(x2 x1 ) x m2 = +k(x3 x2 ) + k(x1 x2 ) x

koennen daher in Mtrixform geschrieben werden:
M x = kx, x = (x1 , x2 , ..., xN ) m m M : Massenmatrix, M = 0 k k k 2k k K :Kraftmatrix, K =
0 m ...
k 2k k k 2k ... k
Dimensionslose BGLs:
d2 x dt2
= k x (x =
x L, t =
k mt
L:
charakteristische Laenge, z.B. Auslenken des 1. Massenpunkte gemaess ABs
52
EIGENWERTPROBLEME
53
1 1 1 2 1 1 2 1 K= 1 2 ... 1

x = vei t
Loesung mit Exponentialansatz:
2 vei t = Kvei t Kvj = j vj , 2

also reduziert auf EW-Problem mit RWs
j 2
und EVs
vj
Da
symmetrisch, Jacobi-Verfahren anwendbar
N = 10 Massenpunkte (Masse m) mit Federn (Federkonstante k) verbunden 1
0.5
x/L
-0.5 Massenpunkt 1 Massenpunkt 2 Massenpunkt 3 Massenpunkt 10 0 2 4 6 (k/m)1/2 t 8 10 12
-1
Abbildung 8.1: N=10 Massenpunkte (m=Masse) mit Federn (Federkonstante k) verbunden
Betrachte
N = 10
und spezielle ABs:
x1 (t = 0) = L, xj (t = 0) = 0, j = 2...10, cos( j t) + Bj sin( j t)) ( Normalschwin
xj (t = 0) = 0, j = 1...10
Allgemeine Loesung:
gung)
x=
N j=1 vj (Aj
Realisiere ABs:
! x(t = 0) = N vj (Bj ) = 0 i=1 Bj = 0 (weil vj orthogonal, also linear unabhaengig) ! x(t = 0) = N vj Aj = (1, 0, 0, ..., 0) j=1
Multipliziere mit
x=
vk , nutze Orthogonalitaet Ak = vk,1 ( 1. Eintrag des vk ) N j=1 vj vj+1 cos( j t) ( Modell z.B. fuer 1D-Kristall, SchallgrschwinN
Massenpunkten, welches
digkeit im Medium ablesbar) Problemlos verallgemeinerbar fr "Jedes" System von der Form WW)
kleine Bewegungen um die stabile Gleichgewichtslage vollfuehrt, wird durch BGLs
M x = kx
beschrieben. ( Taylor-Entwicklung der evtl. komplizierten
"kleine Schwingungen"
53
EIGENWERTPROBLEME
54
8.5
Bibliotheken fuer Eigenwertprobleme

Numerische Eigenwertverfahren i.A. sehr kompliziert
Verwendung existierender Bibliotheken empfehlenswert
LAPACK, ARPACK, (GSL) gute Bibliotheken enthalten eine Reihe spezieller Funktionen fuer
alle EWs/ wenige spezielle EWs ( z.B. die kleinsten) mit/ohne EV-Berechnung dicht/duenn besetzte Matrizen reelle/komplexe Matrizen symmetrische bzw. hermitische Matrizen/Matrizen ohne diese Eigenschaft
54
INTERPOLATION, EXTRAPOLATION, APPROXIMATION DATENMODELLIERUNG
55
9 Interpolation, Extrapolation, Approximation Datenmodellierung
Problemstellung:
Gegeben: Funktionswerte f0 , f1 , ..., fN einer unbekannten Funktionf (x) bei x0
<
x1 < ... < xN

Gesucht:
bzw.
f (x) in einem Bereich xmin x xmax f (y) fuer festes y = x0 , x1 , ..., xN

Fehlerabschaetzung
(kann natuerlich nur naeherungsweise erreicht werden noetig).
Interpolation, wenn
lation.
x0 xmin , mmax xN
bzw. wenn
x0 y xN ,
Extrapo-
Approximation, wenn gefundenes f nur naeherungsweise f (x0 )
= f0 , ..., f (xN ) =
fN
erfuellt, Interpolation .
Vorgehen:
Modelliere Funktion
mit plausiblem Ansatz, z.B. Polynome, physikalisch mo-
tivierte Funktion, etc.
Viele Theoreme in mathematischer Literatur, die Aussagen, welche Funktionen

gut mit welchen Ansaetzen beschrieben werden koennen... in aller Regel nutzlos!
Bezug zur Physik:

f0 , ..., fN
experimentelle Messwerte
f0 , ..., fN entstammen aufwendigen numerischen Rechnungen/Simulationen analytische Funktion f (x) wuenschenswert (zumindest approximativ)
z.B. Potentialwerte, Streuquerschnitte
9.1
Polynominterpolation
verwende Polynom vom Grad lieren Eindeutige Loesung (unschwer beweisbar) Einfache Konstruktion mit Hilfe von Lagrange-Polynomen
N,
um
N +1
Funktionswerte
f0 , f1 , ..., fN
zu interpo-
lj (x) =
N xxk k=0,k=j xj xk
oensichtlich Polynome von Grad
lj (xk ) = jk
Damit
f (x) =
N j=0 fj lj (x)
Weitere Algorithmen zur Polynominterpolation existieren (eventuell numerisch stabiler, weniger aufwaendig, ...) Polynominterpolation fuer grosse
N (N > 5)
in der Regel nicht zu Empfehlen:
Polynom oszilliert stark, aehnelt nicht der natuerlichen Kurve
55
56
Interpolation von 4 Datenpunkten mit einem Grad-3-Polynom 1.4 1.2 1 0.8 0.6 0.4 0.2 Interpolationspolynom Datenpunkte 0 -1 0 1 2 3 4
Abbildung 9.1: Interpolation von 4 Datenpunkten mit einen Grad-3-Polynom
Interpolation von 10 Datenpunkten mit einem Grad-9-Polynom 1.4 1.2 1 0.8 0.6 0.4 0.2 Interpolationspolynom Datenpunkte 0 -1 0 1 2 3 4
Abbildung 9.2: Interpolation von 10 Datenpunkten mit einen Grad-9-Polynom
9.2
Kubische Spline-Interpolation
Setze
Grad-3-Polynome stueckweise zusammen, so dass
die Funktionswerte
f0 , f1 , ..., fN
an den Uebergangsstellen interpoliert,
die 1. und 2. Ableitung an den Uebergangsstellen stetig sind

Vorteil: Polynomgrad ist unabhaengig von lich starken Oszillationen
auch bei grossen
keine unnatuer-
56
57
Interpolation von 10 Datenpunkten mit einem kubischen natuerlichen Spline 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -1 0 1 2 3 4 kubischer natuerlicher Spline Interpolationspolynom vom Grad 9 Datenpunkte
Abbildung 9.3: Interpolation von 10 Datenpunkten mit einem kubischen natuerlichen Spline
Konstruktion von Splines
Interpolation der
fj
kann bereits mit Grad-1-Polynomen erreicht werden:
yj (x) = Afj + Bfj+1 , A =

9.1)
xj+1 x xj+1 xj ,
B=
xxj xj+1 xj (Lagrange-Polynome, Sektion
Annahme: Neben
f0 , ..., fN
waeren auch 2. Ableitungen
f0 , ..., fN
vorgegeben
Eindeutige Grad-3-Polynome yj (x) = Afj + Bfj+1 + Cfj + Dfj+1

1 C = 1 (A3 A)(xj+1 xj )2 , D = 6 (B 3 B)(xj+1 xj )2 6
Es fehlt noch Stetigkeit der 1. Ableitungen; diese legt
f1 , f2 , ..., fN 1
fj fj1 xj xj1
fest.
yj1 (xj ) = yj (xj ), j = 1, ..., N 1 f0

xj xj1 fj1 6
und
xj+1 xj1 fj 3
xj+1 xj fj+1 6
fj+1 fj xj+1 xj
lineares
Gleichungssystem fuer Bestimmung von
fj
(siehe Kapitel 6 ... sogar tridiagonal).
fN
koennen beliebig vorgegeben werden (z.B.
f0 = fN = 0,
sogenannte
natuerliche Spline)
9.3
Splines und Aehnliches bildet einen eigenen "Wissenschaftszweig" CAGD (computer aided geometric design)= mathematische Beschreibung von Kurven und Flaechen fuer wissenschaftliche/medizinische Visualisierung, Tricklme, Computerspiele.
Methode der kleinsten Fehlerquadrate (least square)

j unterliegen haeug statistischen Schwankungen (experimentelle Daten, Resultate
von Monte-Carlo-Integration oder Simulation), gesuchtes
f (x) soll diese Schwankun-
gen nicht wiederspiegeln
Approximation oft geeigneter als Interpolation
waehle Ansatz
f (x; a0 , a1 , ...)
(Parameter des Ansatzes)
57
58
z.B. wieder ein Polynom von niedrigem Grad, oder oder ...
f (x, a) = a0 + a1 x + a2 x2
f (x, a) = f (x, a) =
a0 x fuer ein Coulomb-artiges Potential a0 a1 x fuer Potential mit beschraenkter Reichweite xe
Bestimme Parameter
durch Minimierung von quadratische Abweichung von
G(a) =
punkt bezueglich Parameter
N j=0 (f (xj ; a)
fj )2
(a) G(a)
vom
j -ten
Daten-
a, a
loese also linear in
=0
f
z.B.

f (x; a) = G(a) =
M k=0 ak fk (x) M N j=0 ( k=0
= xk
wenn Ansatz Polynom
M -ten
Grad
ak fk (xj ) fj )(
M l=0 al fl (xj )
fj )
fj kann in Matrixform geschrieben werden f0 (x0 ) f1 (x0 ) ... fM (x0 ) f0 f0 (x1 ) f1 A= . , b = . . . . . f0 (xN ) fM (xN ) fN (Matrix "hoeher" als "breit", d.h. N > M )
M k=0 ak fk (xj )
G(a) =
N j=0 (
M k=0 Ajk ak
bj )(
M l=0 Ajl al
bj )
(*) in Komponenten
( am G(a) =)2
N j=0 (
M k=0 Ajk ak
bj ) Ajm = 0
AT mj
AT Aa = AT b
lineares Gleichungssystem, welches die Parameter
bestimmt;
Loesung einfach, siehe Kapitel 6 Parameter
nicht-linear in
(*) ist nicht-lineares Gleichungssystem
Loesung schwierig, siehe Kapitel 4
Startwerte der Parameter a sollten bereits in der Naehe des gemachten Ergebnisses
liegen
58
59
Least-squares Approximation von 10 Datenpunkten mit Polynomen von niedrigem Grad 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -1 0 1 2 3 4
approximierendes Polynom, Grad 0 approximierendes Polynom, Grad 1 approximierendes Polynom, Grad 2 kubischer natuerlicher Spline Datenpunkte
Abbildung 9.4: Least squares Approximation von 10 Datenpunkten mit Polynomen von diedrigem Grad.
9.4
2 -
Fitting
kleinen Fehlern
Fehlerbalken der Datenpunkte wurden bisher ignoriert. Fr realistischere Ergebnisse sollten Datenpunkte mit
f stark und solche mit grossen Fehlern
kaum beeinussen
Ersetze Fehlerfunktional
durch Bei kleinem Fehler
2 = Gj :
f (xj ;a)fj 2 N ) j=0 ( Gj 2 gross sonst wird

Fehler des
Gj
muss
f (xj )
nahe bei
fj
liegen ...
j -ten
Datenpunkts, also bei
xj
Wert
fj Gj
Resultierendes/minimales
liefert gleichzeitig ein Guetekriterium fuer den Fit:
Falls
guter Fit (also Physik wird durch Ansatz korrekt beschrieben)
jeder Summand in reduziertes
sollte
O(1) =
sein
2 =
2 N M
d.o.f.
degrees of freedom
2 Falls d.o.f.
>> 1 << 1
schlechter Ansatz Fehlerbalken
Falls d.o.f.
Gj
unrealistisch bzw. Datenpunkte korreliert
theoretische Grundlagen dazu in Literatur ueber Datenanalyse
59
60
Least-squares Approximation von 4 Datenpunkten mit einer Konstante 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -1 0 1 2 3 4
Abbildung 9.5: Least squares Approximation von 4 Datenpunkte mit konstantem Ansatz, also
f (x, a) = a
2-Fit einer Konstante an 4 Datenpunkte (2/d.o.f. = 0.97) 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -1 0 1 2 3 4
Abbildung 9.6:
2 -Fitting
von 4 Datenpunkte mit konstantem Ansatz, also
f (x, a) = a
In den obigen Abbildungen wird die Motivation des
- Verfahren deutlich gegen-
ber dem least-square-Verfahren. Whrend das least-square-Verfahren ohne Bercksichtigung der Gewichtung einzelner Fehler die Konstante approximiert mit deutlicher Abweichung, Gewichtet das
- Verfahren anhand der Fehler die Messwerte und approximiert
ein vernnftiges Ergebnis.
2 = 0=
afj 2 l j=0 ( Gj ) N afj d 2 j=0 G2 da = 2 j
60
61
a=
1/G2 j
j N 2 k=0 1/Gk j (Gewicht
Datenpunkts gross, wenn des
fj
Gj klein)
Fehlerabschaetzung der Fit-Parameter Wenn sich
a: fj , k = 1...K
(k)
zusammensetzt (i.d.R.
1 fj )2 ) 2 und diese noch vorliegen
(fj , Gj )
aus "Einzelmessungen"
fj =
1 K
(k) K k=1 fj , Gj
1 ( K(K1)
(k) K k=1 (fj (k)
Fuehre K
Fits mit Datenpunkten
Formeln zur Mettelwert und
(fj , Gj ), k = 1...K Fehlerbestimmung fuer a
aus, verwende dann analoge
Wenn nur
(fj , Gj )
vorliegt
Resampling: Nimm statisch unabhaengige Normalverteilung fuer
fj
an, ge-
(k) neriere zufaellig fj und verfahre wie oben. 1 1 faktor , nicht K K(K1)!
Fehlerbestimmung nur mit Vor-
61
10
FUNKTIONSMINIMIERUNG, OPTIMIERUNG
62
10 Funktionsminimierung, Optimierung
Problemstellung
Gegeben: Funktion Gesucht: Wert von
f (x) x,
der
lokal oder global minimiert
Aequivalent zur Funktionsmaximierung: Maximum von Ziele:
f (x) ist Minimum von f (x)
Schnell, d.h. mit moeglichst wenig Funktionsaufrufen Eventuell globales Minimum

typische Strategien:
sehr schwierig, i.B. wenn Dimensionswahl D>1:
* *
Wiederholte Minimumssuche mit vielen verschiedenen Startwerten Stoere gefundene (eventuell lokalen) Minima und suche von dort weiter
Terminologie/Klassikation von Punkten picMinMax
A,B,C: lokale Maxima D: globales Maximum (da am Rand des Denitionsbereichs kann
gelten)
f (xD ) = 0
E,F: lokale Minima G: globales Minimum (f

lich)
(xG ) = 0
weil im Inneren des Denitionsbereichs)
a,b,c: lokalisieren Minimum (bei
xG ), siehe Abschnitt 10.1 (nur in D = 1 moeg-
Nullstellensuche vs. Funktionsminimierung
Auf ersten Blick aehnlich
f (x) = 0 f1 , f2 , ...unabhaengige Funktionen nicht klar in welche RichAber: tung in x-Raum zu gehen ist, um gleichzeitig alle Komponenten auf 0 zu setzen
f (x) = 0
Komponenten
f verbunden
( f )1 , ( f )2 , ueber infach entgegen des

einfach in hoehe-
Gradienten nach unten fahren"pic
schwer
in hoeheren Di-
vergleichsweise
ren Dimensionen
mensionen
(schlechte) Idee: Fuehre Nullstellensuche auf Funktionsminimierung zurueck, suche Minimum von
F (x) =
2 j (fj (x))
Nicht gut ... Funktionsminimierung ndet in der Regel lokales Minimum mit das keiner Nullstelle entspeicht. pic
F > 0,
62
10
63
10.1
Golden Section Search in 1D

O.B.d.A. Minimumsuche Ausgangspunkt: Minimum bereits in einem Intervall lokalisiert.
Vergleich mit Nullstellensuche ... Nullstelle zwischen
und
lokalisiert, wenn
f (a) < 0
und
f (b) > 0 a
(und
stetig)
Fuer Extrema sind drei Punkte noetig

Minimum zwischen und
lokalisiert, wenn
f (a) > f (b)
und
f (b) < f (c),
a < b, < c
picEx Nicht klar, ob Minimum links oder rechts von
liegt
Verkleinere Intervall
[a, c]
durch Auswertung von
bei
y, b < y < c
(bzw.
a<y<b
* *
...analog)
Fall 1:
f (y) > f (b) aneu = a, bneu = b, cneu = y f (y) < f (b) aneu = b, bneu = y, cneu = c
picFall1 Fall 2:
picFall2
Wie waehlt man
y,
um das Intervall
[a, c]
moeglichst schnell zu verkleinern?
Deniere relative Intervallaengen:

picIntervallaengen
Relative neue Intervalllaenge von
* *
[aneu , cneu]
Fall 1: Fall 2:
w+z 1w
(Dierenz zwischen
Minimiere schlechtesten Fall durch Gleichsetzen:
w + z = 1 w z = 1 2w = (1 w) w graphisch pic
[a, b]
und
[b, c])
Konvergenzgeschwindigkeit:
Wie muss w urspruenglich gewaehlt werden, damit es in jedem Fall gleich bleibt?
picKonv mit
w=
3 5 2 3 5 2
z = 1 2w w2 3w + 1 = 0 2.618... = 0.381... = 0.38... w=

3 5 bewegt sich die Minimumsuche gegen 2
ist stabiler Fixpunkt (ohne Beweis), d.h. auch bei Start
w =
abseits von
Damit cneu
aneu
3 5 = (1 w)(c a) = (1 )(c a) 2
0.618
w=
3 5 2
also lineare Konvergenz, aber schlechter als z.B. Nullstellensuche mit Bisektion, siehe Kapitel 4.2 Goldener Schnitt
1 5 2
Minimale Genauigkeit bei der Minimumsuche:
63
10
64
In der Gegend des Minimums
xmin : f (x) fmin +
fmin 2 (x
xmin )2
Unterscheidung eines Funktionswertes

f /2(xxmin )2 mehr moeglich, wenn min fmin
double, Kajuete?)
f (x) bei x vom Minimum fmin nicht < ( rel. Genauigkeit, 107 oat, 1016 2fmin 1 )2 fmin
O(1)...?
moeglich!
Minimale Intervalllaenge
x xmin ( O( )
In der Regel nur Genauigkeit von
10.2
Quadratische Interpolation in 1D
Approximiere
mit Parabel duch (a, f (a)), (b, f (b)), (c, f (c)), (a,
b, c
z.B. Punkte
die Minimum lokalisieren) picQuad Minimum der Parabel
y = b
1 (ba2 (f (b)f (c))(bc)2 (f (b)f (a)) 2 (ba)(f (b)f (c))(bc)(f (b)f (a))
nahe dem Minimum von
f,
wenn
pa-

10.3
rabelaehnlich, d.h. relativ glatt bezogen auf durch
a, b, c
deniertes Intervall
Problem: Es gibt immer Fuktionen, die die Parabelaehnlichkeit nicht erfuellen
kombiniere schnelle quadratische Interpolation mit langsamer sicherer Methode
(z.B. Golden-Section-Search) Ausfuehrlicher Beispielalgorithmus: Numerical Recipes, 10.3.
Minimumsuche mit Hilfe von Ableitungen in 1D

viele Moeglichkeiten, Ableitungsinformation zu nutzen Beispiel:
Ausgangspunkt: via Berechne
a, b, c f (a) f
lokalisiertes Minimum
f (b)
und
(falls
f (a) < f (b), f (c) f)
sonst)
schaetze Nullstelle von

picMiniDer
mit einem Schritt des Sekantenverfahrens (Abschnitt
4.3) (Nullstelle entspricht Minimum von
Sekantenverfahren konvergiert besser als linear:

1.618 n+1 = cn
(Abschnitt 4.3); damit auch Minimumsverfahren
Problem auch hier. Konvergenz nicht gesichert; z.B. kann neuer Punkt
serhalb von
aus-
[a, c] liegen kombiniere auch hier mit sicherer Methode (z.B. Golden-Section-Search, Abschnitt 10.1) (Beispiel: Numerical Recipes, 10.4)
10.4
Simplex-Methode (D
> 1)
D > 1,
nicht unbedingt sehr ezient
Einfaches Optimierungsverfahren fuer
geeignet, wenn man schnell (menschliche Zeit) ein Ergebnis will (und das Problem
nicht sehr viel Computerzeit benoetigt)
64
10
65
Simplex
Wird in
deniert
Dimensionen durch
D+1
Punkte
pj
und deren Verbindungslinien
D = 2: D = 3:
Dreieck Tetraeder
Wir betrachten nur nicht-degenerierte Simplices
p1 p0 , p2 p0 , ..., pD p0
Raum auf
linear unabhaengig, spannen den
D-dimensionalen
Lokalisieren eines Minimums (wie in picSimplex
D=1
mit 3 Punkten) in
D>1
nicht moeglich
Grundidee: Simplex kriecht schrittweise den Berg hinunter ... deformiert sich bei der durch
f (x)
denierten Landschaft ... erreicht dadurch ein lokales Minimum
Intiabr Simplex:

p0 =Schaetzwert
geben)
fuer Minimum typische Laengenskala in
pj = p0 + j ej (j :
j -Richtung;
von Anwender vorzu-
Simplexdeformationsschritte:
Schritt 1: Spiegelung des maximalen Punkts

pic1Spiegelung
Schritt 2: Expansion mit Faktor 2 (nur falls

pic2Exp
f (p2,neu ) < f (p0 )) p2,neu

noch immer maximaler
1 Schritt 3: Kontraktion mit Faktor 2 (nur falls

Punkt) pic3Kontraktion
Schritt 4: Mehrfache Kontraktion (nur falls

pic4MultKontraktion wieder von vorne
f (p2,neu ) > f (p2,neu ) mit Faktor
1 2 ):
Abbruchkriterium:
Im Gegensatz zu z.B. wenn

|pD p0 | 1 |p +p0 | 2 D
D=1
sind Abbruchkriterium in
D>1
nicht oensichtlich
<
(vom Anwender vorzugeben)
(p0 ist minimaler,
pD
maximaler Punkt)
dann abbrechen, Minimum
=p0
oder Abbruch, wenn Fortbewegung aller Simplex-Punkte marginal... In jedem Fall zur Sicherheit Simplex-Methode noch einmal starten, am gefundenen Minimum (mit grossen initialen Simplex, wie oben beschrieben) (sollte bei Erfolg das Minimum schnell wiedergeben)
65
10
66
10.5
D > 1-Minimierung
Vorgehen: 1. Starte bei Punkt 2. Gib
durch wiederholte
1D
Minimierung
p. u
vor.
1D
Minimierungsrichtung
3. Minimiere
entlang
u,
also nde
so, dass
f (p + u)
minimal (verwende
1D-
Methode, Abschnitt 10.1 bis 10.3) 4. Ersetze
p p + u. D
Dimensionen) naeherungsweise erreicht (siehe
5. Breche ab, falls Minimum (in Abschnitt 10.4), sonst gehe zu 2)
Schwierigkeit liegt in einer geschickten Wahl der Minimierungsrichtung terscheiden sich Methoden dieses Typs
u;
darin un-
Typisches Problem bei ungeschickter Richtungswahl, verdeutlicht an Beispiel.

D=2 f:
Langgezogenes Paraboloid, Ausrichtung schraeg zu Koordinatenachsen
Minimierungsrichtungen: e1 , e2 , e1 , ...
Folie
viele viele Schritte noetig sind
Bei Wahl von Minimierungsrichtungen

...
u1 =
1 2
1 1
u2 =
1 2
1 1
u1 , u2 ,
wenige Schritte fuehren ins Minimum
Beliebter Anfaenger fehler: Waehle Minimierungsrichtungen stets entlang des Gradienten, d.h. die Richtung des steilsten Abstiegs
leicht besser als e1 , e2 , e1 , ..., also Minimierung entlang von Koordinatenachsen Dennoch i.A. sehr viele Minimierungsschritte erforderlich
Zwei moegliche Strategien zur geschickten Wahl: 1. Waehle neue Richtung so, dass Minimierung bezueglich vorheriger Richtungen erhalten bleibt (sprechweise: neue und alte Richtungen sind konjugiert)
Minimierung wird in spaetestens
Schritten erreicht
Bedingung fuer konjugierte Richtungen
* *
u, v : u-Richtung;
dann
p das Minimum gilt u f (p) = 0

sei entlang der
von
entlang einer Geraden in
Waehle konjugierte Richtung
v so, dass der Gradient in u-Richtung v -Minimierungsrichtung verschwindet (praeziser: auf einer Graden durch p in Richtung v ); also u f (p+v) = 0 fuer alle . nicht fuer beliebige Funktionen f erfuellbar moeglich, wenn f quadratisch 1 f (x) = c bx + 2 xAx f (x) = b + Ax
66
10
67
Ausgangspunkt: Forderung: ( beliebig) wenn Falls
u f (p) = 0 u(b + Ap) = 0 u f (pv) = 0 u(b + A(p + v)) = uAv = 0
v das erfuellt ist es zu u konjugiert f naeherungsweise quadratisch:

1 2 j,k
Taylor-Entwicklung: f (x) f (o) + j j f |o xj +

=c =bj
j k f |o xj xk
=:Ajk
(auch andere Wahl als Entwicklungspunkt f(.) moeglich, z.B. scheinlich guenstiger) Damit Minimum nach
wahrnaehe-
Schritten (=D
=1
Minimierung in
rungsweise konjudierten Richtungen) nur naeherungsweise erreicht.
Aber: Schnelle, d.h. quadratische Konvergenz (kann man zeigen!) Falls reich
unaehnlich zu Paraboloid
konjugierte Richtungen wenig hilf-
2. Richtungen entlang von Taelern
Wenn man fast am Minimum ist, ist

verwenden
stets naeherungsweise quadratisch
... am Ende fuer schnelle hohe Genauigkeit immer konjugierte Richtungen
Detaillierte Beschreibung ezienter Methoden: Numerical Recipes, Kapitel

10.7/10.8
10.6
Simulated Annealing
Bisher nur das Finden lokaler Minima diskutiert Wie ndet man das globale Minimum einer Funktion
f?
Wiederholte Minimierung mit verschiedenen Startwerten ... wenn immer das

gleiche Minimum/wenige identische Minima gefunden werden mag das globale dabei sein
Haeug aber riesige Menge von lokalen Minima in unueberschaubaren hochdimensionalen Funktionen
f (x)
... was tun?
10.6.1
Simulated Annealing
Kombinatorische Minimierung
statt kontinuierlichen Wertebereich beschrieben durch Zustnde bzw. Kongurationen
xgewaltige
Menge diskreter
s S = {S1 , S2 , S3 , ..., SN }
Ziel: Finde
31
und mit der Bewertungsfunktion
f (s).
Sj ,
dass f minimiert.
Denn Simulated-Annealing zu Grunde liegende Analogie: eine sich abkuehlende kristallisierende Fluessigkeit, bzw. ausgluehender (Englisch: to anneal) Stahl:
31 N = 10100 z.B., so dass es unmglich ist alle Sj durchzuprobieren.
67
10
68
schnelle Abkuehlung:
Nicht genug Zeit fuer Atome, um sich auszurichten und stabilen Kristall zu
bilden ... sehr zerbrechlich
langsame Abkuehlung:
Atome bilden einheitlichen Kristall (das energetische Minimum der vielen
moeglichen Atomanordnungen) ... sehr stabil Uebertragung auf Minimierung von
f (s): s0 s1 s2 ...
wobei
Bewegung durch Kongurationsraum, che Kongurationen sind
sj
und
sj+1
aehnli-
Akzeptiere bewusst hin und wieder eine Verschlechterung, d.h.

berwunden werden)
f (sj+1 ) > f (sj )
(oft kann ein globales Minimum nur erreicht werden, wenn einige Berge ue-
Akzeptiere Verschlechterung zu Begin haeug, dann mit immer kleinerer Wahrscheinlichkeit (Stahl-Analogie: Wenn heiss, ordnen sich Atome bereitwillig um, ueberwinden Energierbarrieren, wenn kalt, frieren sie kristallisiert in einem Energieminimum fest)
Simulated-Annealing-Algorithmus
1. starte mit beliebigem
sS
und geeigneter Temperator
T
(stochastische Minimie-
(a) Waehle zufaellig aehnliches rung) (b) Falls
sneu S ,
d.h. sneu
s
f (sneu ) f (s) s = sneu

sonst
s = sneu
(c) Verringere (d) Gehe zu 2.
mit Wahrscheinlichkeit
e
E
f (sneu )f (s) T
(Analogie zu Boltzmann-Faktor
e kT
langsam (in Praxis Verringerung immer erst nach vielen Schritten)
Wenn
hinreichend klein friert Algorithmus in einem Minimum ein; wenn Parameter
geschickt gewaehlt (richtige
T -Verringerung, vernuenftige Zustandsaenderungen s
sneu
lange Laufzeit) haeug in globalem Minimum)
Beispiel: Traveling Salesperson problem Gegeben:
Staedte auf einer 2D-Landkarte also
(xj , yj ), j = 1...N
Gesucht: kuerzester Rundweg (alle Staedte muessen einmal besucht werden)
S enthaelt alle Permutationen von (1, 2, ..., N ), N = 100 eine unueberschaubare Anzahl f=
100 j=1 ((xj
1
hat also
N!
Elemente bei z.B.
xj1 )2 + (yj yj1 )2 ) 2 (x0 , y0 x100 , y100 )
Anwendung von Simulated Annealing:
Zustandsaenderungen
s snew :
68
11
MONTE CARLO-SIMULATION STATISTISCHER ZUSTANDSSUMMEN
69
* * * * *
Teilstueck ausschneiden und umdrehen: ...-17-3-9-6-14-...
...-17-6-9-3-14-...
Teilstueck ausschneiden und an anderer Stelle einfuegen ...-5-3-9-1-16-14-21-...
...-5-16-14-3-9-1-21-...
Temperaturveraenderung:
waehle initiales
groesser als typisches
|f (sneu ) f (s)|
Schritten oder
(generiere einige
zufaellige Konguartionen, um das abzuschaetzen) Verringere
um 10% nach entweder
100N
10N
erfolgreichen
Updates (was auch immer eher erfuellt ist) Abbruch, wenn Algorithmus festfriert Folie Modikation:
f f+
(j j1 )2 , j=+1
falls Stadt links vom Fluss,
j = 1 sonst > 0 Angst den Fluss zu ueberqueren < 0 Schmuggler will dauernd Fluss
10.7
ueberqueren
Kontinuierliche Minimierung
Simulated Annealing ohne grosse Modikationen auf Minimierung kontinuierlicher Funktionen
f (x)
anwendbar
Kongurationen
Punkt
im Denitionsbereich von
Zufaellige eziente Wahl von

ren
xnew
schwierig
Haeuges Problem: Es werden zu selten Richtungen gewaehlt, die bergab fueh-
siehe Numerical Recipes, Kapitel 10.12.2
11 Monte Carlo-Simulation statistischer Zustandssummen

11.1
Ising Modell
Haeug verwendetes primitives Modell fuer Magnetismus Auf regelmaessiges kubisches Gitter angeordnete Spins
sj = 1
in
Dimensionen
picIsing Eine spezielle Konguration / Zustand:
C = (s1 , s2 , s3 , ...)
2N Zustaende bei
2)
Spins
Beispiel: Winziger 3D Magnet mit 100 Spins in jeder Raumrichtung
2(100
= 21000000 10300000 H=
moegliche Zustaende
Hamilton-Operator:
> 0fuer
Ferromagneten
sj sk
<j,k>
ueber naechste Nachbarn
j sj
aeusseres Magnetfeld
69
11
70

11.2
Zustandssumme:
Z=
eH(C) ( =
1 kT )
statistischer Erwartungswert Observablen
O(s1 , s2 , ...)
=C
< O >=
O(C)eH(C) Z
Problem: Summen koennen im Prinzip vom Computer berechnet werden, in der Praxis aber unmoeglich, da Summe z.B. aus
10300000
Termen besteht
Verwende sogenannte Monte Carlo-Simulationen
Grundlagen der Monte Carlo-Simulation

MC-Simulationen basieren auf Markov-Ketten
In diesem Kontext ist eine Markov-Kette ein Zufallsgenerator fuer Zustaende
c1 , c2 , c3 , ..., cN
Damit
c eH(J ) Z
< O >
1 N
N c j=1 O(j ) (Fehlerrechnung noetig)
N 10
keiten
<<
Anzahl der moeglichen Zustaende
M 10300000
Markov-Kette generiert in der Regel wahrscheinliche / repraesentative Zustaende Markov-Kette besteht aus Zustaenden
c1 , c2 , ..., cM
(*)
und Uebergangswahrscheinlich-
Pjk , j, k 1, ..., M
k
picMarkovChain Bedingung
Pjk = 1
springe von Zustand zu Zustand gemaess der Uebergangswahrscheinlichkeiten konstruire
Pjk
Pjk
so, dass nach hinreichend langer Laufzeit
P (cj ) =
(**)
eH(cj ) gilt. Z
(lange Laufzeit: Entwicklung muss von Startzustand unabhaengig sein) Dafuer muss fuer
Pjk
gelten
j
e
P (cj ) Pjk = P (ck )

H(cj ) Z eH(ck ) Z
(notwendig aber nicht hinreichend) Haeug konstruiert man
Pjk
(***)
so, dass detailed balance erfuellt ist:
P (cj )Pjk = p(ck )Pkj
(staerkere Bedingung als (**))
Beweis:
j (
)
j
P (cj )Pjk =
P (ck )Pkj Pjk

gemaess (***) waehlen
im Folgenden zwei typische MC-Algorithmen, die
70
11
71
11.2.1
Metropolis-Algorithmen
Schlage zufaellig in
cl = cj
einen neuen Zustand
ck
wjk
vor
(hierbei muss Eigenschaft
wjk = wkj
(einfach zu erfuellen) erfuellt sein)
Akzeptiere neuen Zustand mit Wahrscheinlichkeit
min(1, e(H(ck )H(cj )) ),

Beweis von (***):
also
cl+1 = ck ,
somit behalte alten Zustand, also
cl+1 = cj
P (cj )
e H(cj ) Z
Pjk
wjk min(1,e(H(ck )H(cj )) )
= P (ck )
eH(ck ) Z
Pkj
wkj min(e(H(cj )H(ck )) ,1)
1=1
Damit
c1 , c2 , ...
moeglichst schnell eine repraesentative Zustandsmenge bilden, muss
Foldendes gelten:
Der gemaess
wjk
vorgeschlagene Zustand
ck
sollte sich moeglichst stark von
cj unterscheiden schnelle Bewegung

durch Zustandsraum
ck ck
darf nicht zu oft abgelehnt werden, sonst langsame Bewegung durch den
Zustandsraum (Akzeptanzrate aehnlich zu
50%
i.d.R. vernuenftig) (impliziert oft, dass
cj
ist)
Vor-/Nachteile:
sehr einfach zu realisieren sehr allgemein anwendbar langsame Bewegung durch Zustandsraum Ezienz haengt stark von
wjk
ab; i.d.R. Rumprobieren fuer
wjk -Optimierung
Spins
erforderlich Moegliche Anwendung auf Ising-Modell: Waehle zufaellig aus, stelle sie zufaellig (50% +1, 50 % -1) (damit
wjk
n << Spinanzahl = wkj ) ein
11.2.2
Heatbath-Algorithmus
stelle einen (oder wenn mathematisch moeglich auch mehrere) zufaellig gewaehlten Freiheitsgrad (beim Ising-Modell einen Spin) entsprechend der Wahrscheinlichkeiten
eH(c) ein Z
(alle anderen FHGs werden dabei festgehalten) ein FHG wird bei festgehaltenem anderen ins thermische Gleichgewicht gebracht ( Heatbath)
Pjk =
P (ck ) l P (cl )
P (ck )
Summe ueber Zustande, die durch Veraenderung des ausgewaehlten FHGs erreicht
werden koennen Nachweis (***):
P (cj )
Pjk
= P (ck ) Pkj 1 = 1
P (ck ) 1 P (cj ) #
# der FHGs
1
71
11
72
Vor-/Nachteile
i.d.R. sehr schnelle Bewegung durch den Zustandsraum Uebergangswahrscheinlichkeiten
Pjk
haeug nur schwer/aufwaendig/gar nicht
berechenbar (z.B. kontinuierliche Zustandsraeume)
11.3
ein
Heatbath nur fuer wenige einfache System anwendbar
Anwendung auf Ising-Modell: Waehle zufaellig einen Spin, stelle ihn gemaess
Pjk
Monte Carlo-Simulation des Ising-Modells

Im Folgenden:
2 Dimensionen, 20x20 Spins
B=0 T
(besser
Verschiedene
T /,
dimensionslose Temperatur)
Heatbath-Algorithmus
Thermalisierung durch unabhaengige Starts brechen (heiss (sj zufaellig), kalt (sj
+1))
Folie 1 Approximative Bestimmung des Phasenuebergangs (= Curie-Temperatur) mit Hilfe der Observablen
< (sj < sk >)2 >

lich
(<
sk >:
Mittelwert auf aktueller Konguration)
Fuer praezieseres Resultat Extrapolation zu unendlich vielen FHGs/ Spins erforder-
kein Phasenuebergang im endlichen System
starker Anstieg bei
2.0...T /...2.5
Folie 2 Bestaetigung des analytischen Ergebnisses
Tcrit
2 ln(1+ 2)
Einzelne Spinkongurationen: Folie 3
Metastabiler Zustand (Weisssche Bezirke) (bei heissem Start und Megnetisierung fuer
T < Tc )
T < Tc T > Tc
keine Magnetisierung fuer
72

Numerik Version 0.6

Hochgeladen von

Copyright:

Verfügbare Formate

Numerik Version 0.6

Hochgeladen von

Dokumentinformationen

Copyright

Verfügbare Formate

Dieses Dokument teilen

Dokument teilen oder einbetten

Freigabeoptionen

Stufen Sie dieses Dokument als nützlich ein?

Sind diese Inhalte unangemessen?

Copyright:

Verfügbare Formate

Numerik Version 0.6

Hochgeladen von

Copyright:

Verfügbare Formate

Vorlesung

Numerische Methoden der Physik

Korrektor Martin Stein

Darstellung von Zahlen und Rundungsfehler

Zwei Standardtypen: oat (32 Bits) , double (64 Bits)

Rundungsfehler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Elementare Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerische Ableitung . . . . . . . . . . . . . . . . . . . . . . . . . .

Gewhnliche DGL s & Anfangswertprobleme

Nullstellensuche, lsen nicht-linearer Gleichungen

Vor- und Nachteile: . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Newton-Raphson-Verfahren (fr Newton-Raphson-Verfahren (fr

Gewhnliche DGL s & Randwertprobleme

Beispiel: Quantenmechanik, 1D unendlicher Potentialtopf

Beispiel: Quantenmechanik, 1D harmonischer Osziallator . . . . . . .

Lsung von linearen Gleichungssysteme

6.4.2 6.4.3 6.5 6.6 6.7

Berechnen von det(A) mittels LU . . . . . . . . . . . . . . . . . . . .

QR-Zerlegung (direktes Verfahren)

Interpolation, Extrapolation, Approximation Datenmodellierung

10.6 Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10.6.1 Kombinatorische Minimierung . . . . . . . . . . . . . . . . . . . . . . 10.7 Kontinuierliche Minimierung . . . . . . . . . . . . . . . . . . . . . . . . . . .

11 Monte Carlo-Simulation statistischer Zustandssummen

S = 1 Bit , E = 8 Bits und NM = 23 Bits ist seine Eigenschaften

Stackoverow: Online in Internet: URL: http://stackoverow.com/questions/2795629/cint-oat-double [Stand 15.11.2011]. . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9

Runge-Kutta-Verfahren & die Euler-Methode am Beispiel des Harmonischen

3.10 Nicht-Harmonischer Oszillator in 1D , V(x)=x , n =2,4,6,....

3.11 Runge-Kutta-Verfahrem zweiter Ordnung aufgetragen mit der Algebraischen

mal im dimensionslosen Verhltnis

Ausschnitte einer Animation des Newton-Verfahren 19.11.2011] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 29 30 32 33 34 35 35

Scannen der mglichen Eigenwerte

) numerisch nicht realisierbar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.1 8.1 9.1 9.2 9.3 9.4 9.5 9.6

Ziel: Computer mit Verstand einsetzen (sonst oft sinnlos).

vorliegende Problem zugeschnitten (kein Universalre-

Eigenwertprobleme (Mehrdimensionale) Integration Dierentialgleichungen Nullstellensuche, Optimierung.

und vieles mehr.

Computeralgebrasysteme (Mathematica, Maple, ...)

Computeralgebrasysteme (Mathematica, Maple, ...) ergnzen die Numerik

DARSTELLUNG VON ZAHLEN UND RUNDUNGSFEHLER

2 Darstellung von Zahlen und Rundungsfehler

Computer gister) mit besteht).

Zahlen (+1 Bit fr Sign).

. Analog geht das auch fr die Negativen

Gleitkommazahlen (rationale Zahlen)

Exponent: E ist Integer & e eine Integerkonstante.

DARSTELLUNG VON ZAHLEN UND RUNDUNGSFEHLER

Zwei Standardtypen: oat (32 Bits) , double (64 Bits)

oat: mit S = 1 Bit , E = 8 Bits und NM = 23 Bits ist seine Eigenschaften

Stackoverow: Online in Internet: URL: http://stackoverow.com/questions/2795629/c-int-oatdouble [Stand 15.11.2011].

Wertebereich: M=1,1+ ,1+2 ,..,2-2 ,2wobei die

beste relative Genauigkeit charakterisiert und damit

Die kleinste und grte Zahl durch

1038 ...1038 determiniert

Besitzt statt 4 Byte jetzt 8 Byte (64Bit) also dop-

2Ee 10308 bis 10308

4 Gre 4Bytes 5 Gre 8Bytes

fr den Wertebreich : 10 fr den Wertebreich : 10

stimmt das Ergebnis nicht ganz berein (unnormalized numbers).

DARSTELLUNG VON ZAHLEN UND RUNDUNGSFEHLER

Subtraktion zweier hnlich groer Zahlen (

Z1 Z2 = Z3 = O(Z1 = O(Z2 1 10 2 10 = 3 10 1.??? 1.??? = O(10n ) O(Z1 )