Computational Physics IProgrammiersprachen: EGAL !!! Jede Sprache hat seine Vor- und Nachteile.

Fortran: weit verbreitet, relativ einfach, schnell,komplexe Datenstrukturen nicht ganz einfach zu realisieren

C: weit verbreitet; man muss etwas mehr arbeiten, um es so schnell wie Fortran zu machen; alles ist möglich, aber nicht objektorientiert

C++ und wie C, aber objektorientiert, nicht ganz einfach zu lernenähnliches alles ist möglich

Java verwandt mit C++, Garbage Collection, langsam

Pascal erinnert an Schulzeiten

Python rapid prototyping, langsam; aber durch Bibliotheken mächtig (und dann schnell)

Ruby verwandt mit Python, nicht ganz so viele Mathe-Bibliotheken

Hier: Fortran, C++, Python

Was wollen wir lernen?

Wir wollen Probleme aus der Physik numerisch lösen.

Das macht aber doch auch die Numerische Mathematik!

Antwort: Ja, aber …

Beispiele physikalischer Probleme:

Newtonsche Bewegungsgleichungen : F = m a

Schrödingergln.:

Maxwell:

Wärmeleitungsgleichung: (parabolisch)

Poisson-Gleichung: (elliptisch)

@

@t✓ = �✓

Wellengleichung: (hyperbolisch)

Navier-Stokes:

hyperbolisch elliptisch parabolisch

Vlasov-Gleichung:

Statistik und QFT: Monte-Carlo Simulationen: Phasenübergang beim Ising Modell

Numerische Berechnung von Pfadintegralen

Instantonen

Warum wird es schnell sehr aufwändig?

Diskretisieren in Raum und Zeit

Beispiel: 1-Teilchen Schrödingergln.

diskretisiere Raum in alle 3 Raumrichtungen mit N Gitterpunkten (oder Moden)z.B. N=256, also insgesamt

kein Problem (Laptop, fast Handy)

jetzt 2-Teilchen Schrödingergln.

3-Teilchen Schrödingergln.

@

@t

u+ a

@

@x

u = 0

@

@tu+ u ·ru+rp = ⌫�u+ F

@

@tf(x,v, t) + v ·r

x

f +F

m·r

v

f = 0

2563 · 8 Byte = 134 MByte

2563 ⇤ 2563 · 8 Byte = 2251 TByte

2563 ⇤ 2563 ⇤ 2563 · 8 Byte = 37 ⇤ 1021 = 37 Zetta Byte

Inhaltsverzeichnis

1 Grundlagen

1.1 Zahlendarstellungen

Grundlegende Unterscheidungen:

• Ganzzahlen (Integer, int, etc.)

• Fliesskommazahlen (float, double, REAL)

• Spezielle Typen: Komplex, ...

Jede Zahl wird in einer bestimmten Anzahl von bits gespeichert (Typisch: 16, 32, 64,128 bit).

Anzahl Bits und Darstellung entscheiden uber den Zahlenvorrat eines Typs. Insbeson-dere:

• Dynamik: Verhaltnis von (betragsmassig) grosster zu kleinster Zahl.

• Auflosung: Di↵erenz zwischen benachbarten Zahlen.

Fliesskommazahlen sind rationale Zahlen. Halblogarithmische Darstellung auf 2er-Basis:

z = (�1)s · m · 2e

Typisch: “double” (FORTRAN: REAL*8), in 64 bit dargestellt nach IEEE-Standard (s.http://en.wikipedia.org/wiki/IEEE754):

• 1 bit Vorzeichen (s), 52 bit Mantisse (m), 11 bit Exponent (e).

• Dynamik: 11 bit Exponent mit Vorzeichen e = �210...210 = �1024...1024) Dynamik des Exponenten 2�1024...21024 ⇡ 10�308...10308

• relative Auflosung der Mantisse: 2�52 ⇡ 2.224 · 10�16

3

Matlab:

• “Default”-Typ ist “double” (bzw. komplex mit “double”-Darstellung fur Real- u.Imaginarteil)

• Ganzzahlen mussen explizit generiert werden, z.B. “int32(2)”.

• Ganzzahlarithmetik rundet auf oder ab (untypisch)! Vgl. Sprache “C” etc.: NurAbrundung (Abschneiden).

1.2 Diskretisierung

Diskretisierung von Funktionen

Funktion, z.B. f(x) wird durch endliche Menge von Zahlen, z.B. fi

(Zahlindex i), dar-gestellt.

f(x)) fi

, i = 1, 2, 3, ...

Beispiel:

• fi

ist Wert von f(x) an bestimmter Stelle xi

• fi

ist Mittelwert von f(x) in einem bestimmten Intervall x 2 Ii

• fi

ist ein bestimmter Fourerkoe�zient von f(x)

• fi

ist sonstiger Entwicklungskoe�zient in einem anderen Funktionensystem (z.B.Polynome etc.)

Entsprechend: Diskretisierung von Operatoren

Operator, z.B. L[f ]) Rechenvorschrift, z.B. Ld

[fi

]

• Index d: Diskreter Operator hangt von der Art der Diskretisierung von f ab.

• Ld

[fi

] ist i.A. eine Naherung fur L(f)

• Bei gleicher Diskretisierung von f ist die Diskretisierung von L i.A. nicht eindeutig.

4

1.3 Di↵erentiation

Einfachstes Beispiel: Funktion f(x) in einer unabhangigen Variablen x, di↵erenzierbar.

Gesucht: Diskrete Form der Ableitung f 0(x) =: D[f ] (D: “Di↵erential-Operator”) alsRechenvorschrift fur diskrete Funktionswerte f

i

.

Aus Definition der Ableitung (rechtsseitig)

df

dx(x) = f 0(x) := lim

h!0

f(x + h)� f(x)h

Approximation fur die Ableitung f 0 an der Stelle x mit endlichem h:

f 0(x) ⇡ f(x + h)� f(x)h

Dies ist die Vorwarts-Di↵erenz.

1.3.1 Konvergenztest

Beispiel Vorwarts-Di↵erenz:

• bekanntes f(x), f 0(x)

• wahle Stelle x, “Stutzstellenabstand” h > 0

• berechne diskreten Operator Dh

[f ] := f(x+h)�f(x)h

fur kleiner werdendes h (“Li-mes”)

• berechne Diskretisierungsfehler Eh

:= |f 0(x)�Dh

[f ]|

Ergebnis fur alle (di↵baren) f(x) und x:

1. Vorwarts-Di↵erenz ist konvergent: limh!0 D

h

[f ] konvergiert bis zum Rundungs/Darstellungs/Abschneidefehler

2. Vorwarts-Di↵erenz ist konsistent: Dh

[f ] konvergiert gegen den “richtigen” Wertf 0(x)

3. Konvergenz ist 1. Ordnung in h: Eh

geht mindestens wie h gegen null

5

Alternative: Linksseitige (Ruckwarts-)Di↵erenz:

f 0(x) ⇡ f(x)� f(x� h)h

) Ruckwarts-Di↵erenz ebenfalls konsistentkonvergiert ebenfalls mit 1. Ordnung in h

Weitere Alternative: Zentrale Di↵erenz:

f 0(x) ⇡ f(x + h)� f(x� h)2h

) Konvergiert mit 2. Ordnung in hDiskretisierungsfehler E

h

/ h2 im Limes h! 0

Bestimmung der Konvergenz-Ordnung aus Taylor-Entwicklung:

f(x + h) = f(x) + hf 0(x) +h2f 00(x)

2+

h3f 000(x)6

+ O(h4)

f(x� h) = f(x)� hf 0(x) +h2f 00(x)

2� h3f 000(x)

6+ O(h4)

Daraus:

f(x + h)� f(x)h

= f 0(x) +hf 00(x)

2+

h2f 000(x)6

+ O(h3) = f 0(x) + O(h)

f(x)� f(x� h)h

= f 0(x)� hf 00(x)2

+h2f 000(x)

6+ O(h3) = f 0(x) + O(h)

aber:f(x + h)� f(x� h)

2h= f 0(x) +

h2f 000(x)6

+ O(h3) = f 0(x) + O(h2)

1.3.2 “Stencil”-Notation

(“Schablone”), deutscher Term: “Stempel” oder “Stern”. Oder Stencil.

6

−12 −10 −8 −6 −4 −2 0−25

−20

−15

−10

−5

0

log h

log

|Feh

ler|

Finite−Differenzen Konvergenz

Abbildung 1: Konvergenz von Vorwarts- (rot), Ruckwarts- (grun) und zentraler(blau) Di↵erenz aus Test an Beispielfunktion.

Stencil bezuglich Stutzstellenabstand h:

1h

⇥0 �1 1

⇤:=

f(x + h)� f(x)h

= Standard Vorwarts-Di↵erenz

1h

⇥�1 1 0

⇤:=

f(x)� f(x� h)h

= Standard Ruckwarts-Di↵erenz

12h

⇥�1 0 1

⇤:=

f(x + h)� f(x� h)2h

= Standard zentrale-Di↵erenz

Bem:12h

⇥�1 0 1

⇤=

12

✓1h

⇥0 �1 1

⇤+

1h

⇥�1 1 0

⇤◆

Elimination des Fehlerterms / h

Zentrale Di↵erenz mit breitem Stencil, ersetze h! 2h:

14h

⇥�1 0 0 0 1

⇤=

f(x + 2h)� f(x� 2h)4h

= f 0(x)+4h2f 000(x)

6+O(h3) = f 0(x)+O(h2)

7

Elimination des h2-Fehlerterms:

✓4 ⇤ 1

2h

⇥�1 0 1

⇤� 1

4h

⇥�1 0 0 0 1

⇤◆= 3 ⇤ f 0(x) + O(h3)

Damit:

f 0(x) =1

12h

⇥1 �8 0 8 �1

⇤+ O(h3)

1.3.3 Hohere Ableitungen

Analoges Vorgehen. Beispiel: 2. Ableitung. Starte wieder mit Taylor:

f(x + h) = f(x) + hf 0(x) +h2f 00(x)

2+

h3f 000(x)6

+ O(h4)

f(x� h) = f(x)� hf 0(x) +h2f 00(x)

2� h3f 000(x)

6+ O(h4)

Daraus Standard 2.-Ableitung mit 3-Punkt Stencil

f 00(x) =f(x + h) + f(x� h)� 2f(x)

h2+ O(h2) =

1h2

⇥1 �2 1

⇤+ O(h2)

1.3.4 Spektralverhalten eines diskreten Operators

Bisher: Feste Testfunktion f(x), h varriert ) Konvergenz

Jetzt: Betrachte verschiedene Funktionen, speziell Fourier-Moden

f(x) = eikx f 0(x) = ikeikx

Beachte: Operator ist linear, damit enthalt Wirkung auf Fourier-Moden praktisch kom-plette Information!

Bestimme f 0 aus zentraler Di↵erenz:

8

f(x + h)� f(x� h)2h

=eik(x+h) � eik(x�h)

2h

=eikx

2h

⇣eikh � e�ikh

⌘=

ieikx

hsin(kh)

= f 0(x)| {z }exakterWert

sin(kh)kh

Spektraler Diskretisierungsfehler: sin(kh)kh

� 1

Amplitudenfaktor der zentralen Di↵erenz:

0.5

x6543210

1.0

0.75

0.25

0.0

sin(x)/x

• langwellige (niederfrequente) Moden kh ⌧ 1: sin(kh)kh

⇡ 1 + O(kh)2, Diskretisie-rungsfehler 2. Ordnung

• Zentrale Di↵. wird null bei kh = ⇡ bzw. Wellenlange � := 2⇡

k

= 2h. Danachfalsches Vorzeichen!

• Betrag der zentralen Di↵. von eikx ist immer 1h

• Spektrum der zentralen Di↵. ist rein reell, kein Phasenfehler.

9